source: palm/trunk/SOURCE/flow_statistics.f90 @ 674

Last change on this file since 674 was 674, checked in by suehring, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 51.5 KB
Line 
1 SUBROUTINE flow_statistics
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7! Former revisions:
8! -----------------
9! $Id: flow_statistics.f90 674 2011-01-18 16:33:31Z suehring $
10!
11! 673 2011-01-18 16:19:48Z suehring
12! Top bc for the horizontal velocity variances added for ocean runs.
13! Setting the corresponding bottom bc moved to advec_ws.
14!
15! 667 2010-12-23 12:06:00Z suehring/gryschka
16! When advection is computed with ws-scheme, turbulent fluxes are already
17! computed in the respective advection routines and buffered in arrays
18! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the
19! numerics and to avoid unphysical kinks near the surface.
20! So some if requests has to be done to dicern between fluxes from ws-scheme
21! other advection schemes.
22! Furthermore the computation of z_i is only done if the heat flux exceeds a
23! minimum value. This affects only simulations of a neutral boundary layer and
24! is due to reasons of computations in the advection scheme.
25!
26! 624 2010-12-10 11:46:30Z heinze
27! Calculation of q*2 added
28!
29! 622 2010-12-10 08:08:13Z raasch
30! optional barriers included in order to speed up collective operations
31!
32! 388 2009-09-23 09:40:33Z raasch
33! Vertical profiles of potential density and hydrostatic pressure are
34! calculated.
35! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
36! end.
37! Temperature gradient criterion for estimating the boundary layer height
38! replaced by the gradient criterion of Sullivan et al. (1998).
39! Output of messages replaced by message handling routine.
40!
41! 197 2008-09-16 15:29:03Z raasch
42! Spline timeseries splptx etc. removed, timeseries w'u', w'v', w'q' (k=0)
43! added,
44! bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr)
45! (like other scalars)
46!
47! 133 2007-11-20 10:10:53Z letzel
48! Vertical profiles now based on nzb_s_inner; they are divided by
49! ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered
50! velocity components and their products, procucts of scalars and velocity
51! components), respectively.
52!
53! 106 2007-08-16 14:30:26Z raasch
54! Prescribed momentum fluxes at the top surface are used,
55! profiles for w*p* and w"e are calculated
56!
57! 97 2007-06-21 08:23:15Z raasch
58! Statistics for ocean version (salinity, density) added,
59! calculation of z_i and Deardorff velocity scale adjusted to be used with
60! the ocean version
61!
62! 87 2007-05-22 15:46:47Z raasch
63! Two more arguments added to user_statistics, which is now also called for
64! user-defined profiles,
65! var_hom and var_sum renamed pr_palm
66!
67! 82 2007-04-16 15:40:52Z raasch
68! Cpp-directive lcmuk changed to intel_openmp_bug
69!
70! 75 2007-03-22 09:54:05Z raasch
71! Collection of time series quantities moved from routine flow_statistics to
72! here, routine user_statistics is called for each statistic region,
73! moisture renamed humidity
74!
75! 19 2007-02-23 04:53:48Z raasch
76! fluxes at top modified (tswst, qswst)
77!
78! RCS Log replace by Id keyword, revision history cleaned up
79!
80! Revision 1.41  2006/08/04 14:37:50  raasch
81! Error removed in non-parallel part (sums_l)
82!
83! Revision 1.1  1997/08/11 06:15:17  raasch
84! Initial revision
85!
86!
87! Description:
88! ------------
89! Compute average profiles and further average flow quantities for the different
90! user-defined (sub-)regions. The region indexed 0 is the total model domain.
91!
92! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
93! ----  lower vertical index for k-loops for all variables, although strictly
94! speaking the k-loops would have to be split up according to the staggered
95! grid. However, this implies no error since staggered velocity components are
96! zero at the walls and inside buildings.
97!------------------------------------------------------------------------------!
98
99    USE arrays_3d
100    USE cloud_parameters
101    USE cpulog
102    USE grid_variables
103    USE indices
104    USE interfaces
105    USE pegrid
106    USE statistics
107    USE control_parameters
108
109    IMPLICIT NONE
110
111    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
112    LOGICAL ::  first
113    REAL    ::  dptdz_threshold, height, pts, sums_l_eper, sums_l_etot, ust, &
114                ust2, u2, vst, vst2, v2, w2, z_i(2)
115    REAL    ::  dptdz(nzb+1:nzt+1)
116    REAL    ::  sums_ll(nzb:nzt+1,2)
117
118    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
119
120!
121!-- To be on the safe side, check whether flow_statistics has already been
122!-- called once after the current time step
123    IF ( flow_statistics_called )  THEN
124
125       message_string = 'flow_statistics is called two times within one ' // &
126                        'timestep'
127       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
128     
129    ENDIF
130
131!
132!-- Compute statistics for each (sub-)region
133    DO  sr = 0, statistic_regions
134
135!
136!--    Initialize (local) summation array
137       sums_l = 0.0
138
139!
140!--    Store sums that have been computed in other subroutines in summation
141!--    array
142       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
143!--    WARNING: next line still has to be adjusted for OpenMP
144       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
145       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
146       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
147
148!
149!--    Copy the turbulent quantities, evaluated in the advection routines to
150!--    the local array sums_l() for further computations
151       IF ( ws_scheme_mom )  THEN
152!       
153!--       According to the Neumann bc for the horizontal velocity components,
154!--       the corresponding fluxes has to satisfiy the same bc.
155          IF ( ocean )  THEN
156             sums_us2_ws_l(nzt+1,sr) = sums_us2_ws_l(nzt,sr)
157             sums_vs2_ws_l(nzt+1,sr) = sums_vs2_ws_l(nzt,sr)   
158          ENDIF
159!         
160!--       Swap the turbulent quantities evaluated in advec_ws.
161          sums_l(:,13,0) = sums_wsus_ws_l(:,sr)       ! w*u*
162          sums_l(:,15,0) = sums_wsvs_ws_l(:,sr)       ! w*v*
163          sums_l(:,30,0) = sums_us2_ws_l(:,sr)        ! u*2
164          sums_l(:,31,0) = sums_vs2_ws_l(:,sr)        ! v*2
165          sums_l(:,32,0) = sums_ws2_ws_l(:,sr)        ! w*2
166          sums_l(:,34,0) = sums_l(:,34,0) + 0.5 *                             & 
167                (sums_us2_ws_l(:,sr) + sums_vs2_ws_l(:,sr)                    &
168                + sums_ws2_ws_l(:,sr))                      ! e*
169          DO  k = nzb, nzt
170             sums_l(nzb+5,pr_palm,0) = sums_l(nzb+5,pr_palm,0) + 0.5 * (      &
171                sums_us2_ws_l(k,sr) + sums_vs2_ws_l(k,sr) +                   &
172                sums_ws2_ws_l(k,sr))
173          ENDDO
174       ENDIF
175       IF ( ws_scheme_sca )  THEN
176          sums_l(:,17,0) = sums_wspts_ws_l(:,sr)      ! w*pt* from advec_s_ws
177          IF ( ocean ) sums_l(:,66,0) = sums_wssas_ws_l(:,sr) ! w*sa*
178          IF ( humidity  .OR.  passive_scalar ) sums_l(:,49,0) =              &
179                                                   sums_wsqs_ws_l(:,sr) !w*q*
180       ENDIF
181!
182!--    Horizontally averaged profiles of horizontal velocities and temperature.
183!--    They must have been computed before, because they are already required
184!--    for other horizontal averages.
185       tn = 0
186
187       !$OMP PARALLEL PRIVATE( i, j, k, tn )
188#if defined( __intel_openmp_bug )
189       tn = omp_get_thread_num()
190#else
191!$     tn = omp_get_thread_num()
192#endif
193
194       !$OMP DO
195       DO  i = nxl, nxr
196          DO  j =  nys, nyn
197             DO  k = nzb_s_inner(j,i), nzt+1
198                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
199                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
200                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
201             ENDDO
202          ENDDO
203       ENDDO
204
205!
206!--    Horizontally averaged profile of salinity
207       IF ( ocean )  THEN
208          !$OMP DO
209          DO  i = nxl, nxr
210             DO  j =  nys, nyn
211                DO  k = nzb_s_inner(j,i), nzt+1
212                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
213                                      sa(k,j,i) * rmask(j,i,sr)
214                ENDDO
215             ENDDO
216          ENDDO
217       ENDIF
218
219!
220!--    Horizontally averaged profiles of virtual potential temperature,
221!--    total water content, specific humidity and liquid water potential
222!--    temperature
223       IF ( humidity )  THEN
224          !$OMP DO
225          DO  i = nxl, nxr
226             DO  j =  nys, nyn
227                DO  k = nzb_s_inner(j,i), nzt+1
228                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
229                                      vpt(k,j,i) * rmask(j,i,sr)
230                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
231                                      q(k,j,i) * rmask(j,i,sr)
232                ENDDO
233             ENDDO
234          ENDDO
235          IF ( cloud_physics )  THEN
236             !$OMP DO
237             DO  i = nxl, nxr
238                DO  j =  nys, nyn
239                   DO  k = nzb_s_inner(j,i), nzt+1
240                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
241                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
242                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
243                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
244                                                          ) * rmask(j,i,sr)
245                   ENDDO
246                ENDDO
247             ENDDO
248          ENDIF
249       ENDIF
250
251!
252!--    Horizontally averaged profiles of passive scalar
253       IF ( passive_scalar )  THEN
254          !$OMP DO
255          DO  i = nxl, nxr
256             DO  j =  nys, nyn
257                DO  k = nzb_s_inner(j,i), nzt+1
258                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
259                ENDDO
260             ENDDO
261          ENDDO
262       ENDIF
263       !$OMP END PARALLEL
264!
265!--    Summation of thread sums
266       IF ( threads_per_task > 1 )  THEN
267          DO  i = 1, threads_per_task-1
268             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
269             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
270             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
271             IF ( ocean )  THEN
272                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
273             ENDIF
274             IF ( humidity )  THEN
275                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
276                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
277                IF ( cloud_physics )  THEN
278                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
279                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
280                ENDIF
281             ENDIF
282             IF ( passive_scalar )  THEN
283                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
284             ENDIF
285          ENDDO
286       ENDIF
287
288#if defined( __parallel )
289!
290!--    Compute total sum from local sums
291       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
292       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
293                           MPI_SUM, comm2d, ierr )
294       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
295       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
296                           MPI_SUM, comm2d, ierr )
297       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
298       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
299                           MPI_SUM, comm2d, ierr )
300       IF ( ocean )  THEN
301          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
302          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
303                              MPI_REAL, MPI_SUM, comm2d, ierr )
304       ENDIF
305       IF ( humidity ) THEN
306          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
307          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
308                              MPI_REAL, MPI_SUM, comm2d, ierr )
309          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
310          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
311                              MPI_REAL, MPI_SUM, comm2d, ierr )
312          IF ( cloud_physics ) THEN
313             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
314             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
315                                 MPI_REAL, MPI_SUM, comm2d, ierr )
316             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
317             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
318                                 MPI_REAL, MPI_SUM, comm2d, ierr )
319          ENDIF
320       ENDIF
321
322       IF ( passive_scalar )  THEN
323          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
324          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
325                              MPI_REAL, MPI_SUM, comm2d, ierr )
326       ENDIF
327#else
328       sums(:,1) = sums_l(:,1,0)
329       sums(:,2) = sums_l(:,2,0)
330       sums(:,4) = sums_l(:,4,0)
331       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
332       IF ( humidity ) THEN
333          sums(:,44) = sums_l(:,44,0)
334          sums(:,41) = sums_l(:,41,0)
335          IF ( cloud_physics ) THEN
336             sums(:,42) = sums_l(:,42,0)
337             sums(:,43) = sums_l(:,43,0)
338          ENDIF
339       ENDIF
340       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
341#endif
342
343!
344!--    Final values are obtained by division by the total number of grid points
345!--    used for summation. After that store profiles.
346       sums(:,1) = sums(:,1) / ngp_2dh(sr)
347       sums(:,2) = sums(:,2) / ngp_2dh(sr)
348       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
349       hom(:,1,1,sr) = sums(:,1)             ! u
350       hom(:,1,2,sr) = sums(:,2)             ! v
351       hom(:,1,4,sr) = sums(:,4)             ! pt
352
353
354!
355!--    Salinity
356       IF ( ocean )  THEN
357          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
358          hom(:,1,23,sr) = sums(:,23)             ! sa
359       ENDIF
360
361!
362!--    Humidity and cloud parameters
363       IF ( humidity ) THEN
364          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
365          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
366          hom(:,1,44,sr) = sums(:,44)             ! vpt
367          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
368          IF ( cloud_physics ) THEN
369             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
370             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
371             hom(:,1,42,sr) = sums(:,42)             ! qv
372             hom(:,1,43,sr) = sums(:,43)             ! pt
373          ENDIF
374       ENDIF
375
376!
377!--    Passive scalar
378       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /  &
379            ngp_2dh_s_inner(:,sr)                    ! s (q)
380
381!
382!--    Horizontally averaged profiles of the remaining prognostic variables,
383!--    variances, the total and the perturbation energy (single values in last
384!--    column of sums_l) and some diagnostic quantities.
385!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
386!--    ----  speaking the following k-loop would have to be split up and
387!--          rearranged according to the staggered grid.
388!--          However, this implies no error since staggered velocity components
389!--          are zero at the walls and inside buildings.
390       tn = 0
391#if defined( __intel_openmp_bug )
392       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
393       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
394       tn = omp_get_thread_num()
395#else
396       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
397!$     tn = omp_get_thread_num()
398#endif
399       !$OMP DO
400       DO  i = nxl, nxr
401          DO  j =  nys, nyn
402             sums_l_etot = 0.0
403             DO  k = nzb_s_inner(j,i), nzt+1
404!
405!--             Prognostic and diagnostic variables
406                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
407                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
408                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
409                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
410                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
411
412                sums_l(k,33,tn) = sums_l(k,33,tn) + &
413                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
414
415                IF ( humidity )  THEN
416                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
417                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
418                ENDIF
419
420                sums_l_etot  = sums_l_etot + &
421                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
422                                        w(k,j,i)**2 ) * rmask(j,i,sr)
423             ENDDO
424!
425!--          Total and perturbation energy for the total domain (being
426!--          collected in the last column of sums_l). Summation of these
427!--          quantities is seperated from the previous loop in order to
428!--          allow vectorization of that loop.
429             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
430!
431!--          2D-arrays (being collected in the last column of sums_l)
432             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +   &
433                                        us(j,i)   * rmask(j,i,sr)
434             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
435                                        usws(j,i) * rmask(j,i,sr)
436             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
437                                        vsws(j,i) * rmask(j,i,sr)
438             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
439                                        ts(j,i)   * rmask(j,i,sr)
440             IF ( humidity )  THEN
441                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
442                                            qs(j,i)   * rmask(j,i,sr)
443             ENDIF
444          ENDDO
445       ENDDO
446
447!
448!--    Computation of statistics when ws-scheme is not used. Else these
449!--    quantities are evaluated in the advection routines.
450       IF ( .NOT. ws_scheme_mom )  THEN
451          !$OMP DO
452          DO  i = nxl, nxr
453             DO  j =  nys, nyn
454                sums_l_eper = 0.0
455                DO  k = nzb_s_inner(j,i), nzt+1
456                   u2   = u(k,j,i)**2
457                   v2   = v(k,j,i)**2
458                   w2   = w(k,j,i)**2
459                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
460                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
461
462                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
463                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
464                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
465!
466!--   Higher moments
467!--  (Computation of the skewness of w further below)
468                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
469                                                    rmask(j,i,sr)
470!
471!--             Perturbation energy
472
473                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
474                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
475                   sums_l_eper  = sums_l_eper + &
476                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
477
478                ENDDO
479                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
480                     + sums_l_eper
481             ENDDO
482          ENDDO
483       ELSE
484          !$OMP DO
485          DO  i = nxl, nxr
486             DO  j =  nys, nyn
487                DO  k = nzb_s_inner(j,i), nzt + 1
488                   w2   = w(k,j,i)**2
489!
490!--                Higher moments
491!--                (Computation of the skewness of w further below)
492                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
493                                                    rmask(j,i,sr)
494                ENDDO
495             ENDDO
496          ENDDO
497       ENDIF
498
499!
500!--    Horizontally averaged profiles of the vertical fluxes
501
502       !$OMP DO
503       DO  i = nxl, nxr
504          DO  j = nys, nyn
505!
506!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
507!--          oterwise from k=nzb+1)
508!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
509!--          ----  strictly speaking the following k-loop would have to be
510!--                split up according to the staggered grid.
511!--                However, this implies no error since staggered velocity
512!--                components are zero at the walls and inside buildings.
513
514             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
515!
516!--             Momentum flux w"u"
517                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * (                   &
518                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
519                                                           ) * (               &
520                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
521                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
522                                                               ) * rmask(j,i,sr)
523!
524!--             Momentum flux w"v"
525                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * (                   &
526                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
527                                                           ) * (               &
528                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
529                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
530                                                               ) * rmask(j,i,sr)
531!
532!--             Heat flux w"pt"
533                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
534                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
535                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
536                                               * ddzu(k+1) * rmask(j,i,sr)
537
538
539!
540!--             Salinity flux w"sa"
541                IF ( ocean )  THEN
542                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
543                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
544                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
545                                               * ddzu(k+1) * rmask(j,i,sr)
546                ENDIF
547
548!
549!--             Buoyancy flux, water flux (humidity flux) w"q"
550                IF ( humidity ) THEN
551                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
552                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
553                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
554                                               * ddzu(k+1) * rmask(j,i,sr)
555                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
556                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
557                                               * ( q(k+1,j,i) - q(k,j,i) )     &
558                                               * ddzu(k+1) * rmask(j,i,sr)
559                   IF ( cloud_physics ) THEN
560                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
561                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
562                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
563                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
564                                               * ddzu(k+1) * rmask(j,i,sr) 
565                   ENDIF
566                ENDIF
567
568!
569!--             Passive scalar flux
570                IF ( passive_scalar )  THEN
571                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
572                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
573                                               * ( q(k+1,j,i) - q(k,j,i) )     &
574                                               * ddzu(k+1) * rmask(j,i,sr)
575                ENDIF
576
577             ENDDO
578
579!
580!--          Subgridscale fluxes in the Prandtl layer
581             IF ( use_surface_fluxes )  THEN
582                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
583                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
584                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
585                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
586                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
587                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
588                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
589                                    0.0 * rmask(j,i,sr)           ! u"pt"
590                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
591                                    0.0 * rmask(j,i,sr)           ! v"pt"
592                IF ( ocean )  THEN
593                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
594                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
595                ENDIF
596                IF ( humidity )  THEN
597                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
598                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
599                   IF ( cloud_physics )  THEN
600                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
601                                          ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
602                                          shf(j,i) + 0.61 * pt(nzb,j,i) * &
603                                                     qsws(j,i)            &
604                                                              )
605!
606!--                   Formula does not work if ql(nzb) /= 0.0
607                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
608                                          qsws(j,i) * rmask(j,i,sr)
609                   ENDIF
610                ENDIF
611                IF ( passive_scalar )  THEN
612                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
613                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
614                ENDIF
615             ENDIF
616
617!
618!--          Subgridscale fluxes at the top surface
619             IF ( use_top_fluxes )  THEN
620                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
621                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
622                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
623                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
624                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
625                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
626                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
627                                    0.0 * rmask(j,i,sr)           ! u"pt"
628                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
629                                    0.0 * rmask(j,i,sr)           ! v"pt"
630
631                IF ( ocean )  THEN
632                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
633                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
634                ENDIF
635                IF ( humidity )  THEN
636                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
637                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
638                   IF ( cloud_physics )  THEN
639                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
640                                          ( 1.0 + 0.61 * q(nzt,j,i) ) *   &
641                                          tswst(j,i) + 0.61 * pt(nzt,j,i) * &
642                                                     qsws(j,i)            &
643                                                              )
644!
645!--                   Formula does not work if ql(nzb) /= 0.0
646                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
647                                          qswst(j,i) * rmask(j,i,sr)
648                   ENDIF
649                ENDIF
650                IF ( passive_scalar )  THEN
651                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
652                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
653                ENDIF
654             ENDIF
655
656!
657!--          Resolved fluxes (can be computed for all horizontal points)
658!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
659!--          ----  speaking the following k-loop would have to be split up and
660!--                rearranged according to the staggered grid.
661             DO  k = nzb_s_inner(j,i), nzt
662                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
663                              u(k+1,j,i) - hom(k+1,1,1,sr) )
664                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
665                              v(k+1,j,i) - hom(k+1,1,2,sr) )
666                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
667                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
668
669!--             Higher moments
670                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
671                                                    rmask(j,i,sr)
672                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
673                                                    rmask(j,i,sr)
674
675!
676!--             Salinity flux and density (density does not belong to here,
677!--             but so far there is no other suitable place to calculate)
678                IF ( ocean )  THEN
679                   IF( .NOT. ws_scheme_sca )  THEN
680                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
681                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
682                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
683                                                       rmask(j,i,sr)
684                   ENDIF
685                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
686                                                       rmask(j,i,sr)
687                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
688                                                       rmask(j,i,sr)
689                ENDIF
690
691!
692!--             Buoyancy flux, water flux, humidity flux and liquid water
693!--             content
694                IF ( humidity )  THEN
695                   pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
696                                 vpt(k+1,j,i) - hom(k+1,1,44,sr) )
697                   sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
698                                                       rmask(j,i,sr)
699
700                   IF ( cloud_physics  .OR.  cloud_droplets )  THEN
701                      pts = 0.5 *                                            &
702                             ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) &
703                             + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
704                      sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
705                                                          rmask(j,i,sr)
706                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
707                                                          rmask(j,i,sr)
708                   ENDIF
709                ENDIF
710
711!
712!--             Passive scalar flux
713                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca ))  THEN
714                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
715                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
716                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
717                                                       rmask(j,i,sr)
718                ENDIF
719
720!
721!--             Energy flux w*e*
722!--             has to be adjusted
723                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
724                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
725                                             * rmask(j,i,sr)
726             ENDDO
727          ENDDO
728       ENDDO
729!-     for reasons of speed optimization the loop is splitted, to avoid if-else
730!-     statements inside the loop
731!-     Fluxes which have been computed in part directly inside the advection routines
732!-     treated seperatly.
733!-     First treat the momentum fluxes
734       IF ( .NOT. ws_scheme_mom )  THEN
735         !$OMP DO
736         DO  i = nxl, nxr
737            DO  j = nys, nyn
738               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
739                  ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
740                              u(k+1,j,i) - hom(k+1,1,1,sr) )
741                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
742                              v(k+1,j,i) - hom(k+1,1,2,sr) )
743!                             
744!--               Momentum flux w*u*
745                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
746                                                    ( w(k,j,i-1) + w(k,j,i) ) &
747                                                    * ust * rmask(j,i,sr)
748!
749!--               Momentum flux w*v*
750                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
751                                                    ( w(k,j-1,i) + w(k,j,i) ) &
752                                                    * vst * rmask(j,i,sr)
753               ENDDO
754            ENDDO
755         ENDDO
756
757       ENDIF
758       IF ( .NOT. ws_scheme_sca )  THEN
759         !$OMP DO
760         DO  i = nxl, nxr
761            DO  j = nys, nyn
762               DO  k = nzb_diff_s_inner(j,i) - 1, nzt_diff
763!-                vertical heat flux
764                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
765                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
766                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
767                           * w(k,j,i) * rmask(j,i,sr)
768                  IF ( humidity )  THEN
769                     pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
770                                q(k+1,j,i) - hom(k+1,1,41,sr) )
771                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
772                                                      rmask(j,i,sr)
773                  ENDIF
774               ENDDO
775            ENDDO
776         ENDDO
777
778       ENDIF
779
780
781!
782!--    Density at top follows Neumann condition
783       IF ( ocean )  THEN
784          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
785          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
786       ENDIF
787
788!
789!--    Divergence of vertical flux of resolved scale energy and pressure
790!--    fluctuations as well as flux of pressure fluctuation itself (68).
791!--    First calculate the products, then the divergence.
792!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
793       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
794
795          sums_ll = 0.0  ! local array
796
797          !$OMP DO
798          DO  i = nxl, nxr
799             DO  j = nys, nyn
800                DO  k = nzb_s_inner(j,i)+1, nzt
801
802                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
803                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
804                           - 2.0 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
805                           ) )**2                                          &
806                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
807                           - 2.0 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
808                           ) )**2                                          &
809                   + w(k,j,i)**2                                  )
810
811                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
812                                               * ( p(k,j,i) + p(k+1,j,i) )
813
814                ENDDO
815             ENDDO
816          ENDDO
817          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
818          sums_ll(nzt+1,1) = 0.0
819          sums_ll(0,2)     = 0.0
820          sums_ll(nzt+1,2) = 0.0
821
822          DO  k = nzb_s_inner(j,i)+1, nzt
823             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
824             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
825             sums_l(k,68,tn) = sums_ll(k,2)
826          ENDDO
827          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
828          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
829          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
830
831       ENDIF
832
833!
834!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
835       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
836
837          !$OMP DO
838          DO  i = nxl, nxr
839             DO  j = nys, nyn
840                DO  k = nzb_s_inner(j,i)+1, nzt
841
842                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
843                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
844                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
845                                                             ) * ddzw(k)
846
847                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
848                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
849                                                              )
850
851                ENDDO
852             ENDDO
853          ENDDO
854          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
855          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
856
857       ENDIF
858
859!
860!--    Horizontal heat fluxes (subgrid, resolved, total).
861!--    Do it only, if profiles shall be plotted.
862       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
863
864          !$OMP DO
865          DO  i = nxl, nxr
866             DO  j = nys, nyn
867                DO  k = nzb_s_inner(j,i)+1, nzt
868!
869!--                Subgrid horizontal heat fluxes u"pt", v"pt"
870                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
871                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
872                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
873                                                 * ddx * rmask(j,i,sr)
874                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
875                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
876                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
877                                                 * ddy * rmask(j,i,sr)
878!
879!--                Resolved horizontal heat fluxes u*pt*, v*pt*
880                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
881                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
882                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
883                                                 pt(k,j,i)   - hom(k,1,4,sr) )
884                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
885                                 pt(k,j,i)   - hom(k,1,4,sr) )
886                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
887                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
888                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
889                                                 pt(k,j,i)   - hom(k,1,4,sr) )
890                ENDDO
891             ENDDO
892          ENDDO
893!
894!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
895          sums_l(nzb,58,tn) = 0.0
896          sums_l(nzb,59,tn) = 0.0
897          sums_l(nzb,60,tn) = 0.0
898          sums_l(nzb,61,tn) = 0.0
899          sums_l(nzb,62,tn) = 0.0
900          sums_l(nzb,63,tn) = 0.0
901
902       ENDIF
903
904!
905!--    Calculate the user-defined profiles
906       CALL user_statistics( 'profiles', sr, tn )
907       !$OMP END PARALLEL
908
909!
910!--    Summation of thread sums
911       IF ( threads_per_task > 1 )  THEN
912          DO  i = 1, threads_per_task-1
913             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
914             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
915             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
916                                      sums_l(:,45:pr_palm,i)
917             IF ( max_pr_user > 0 )  THEN
918                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
919                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
920                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
921             ENDIF
922          ENDDO
923       ENDIF
924
925#if defined( __parallel )
926
927!
928!--    Compute total sum from local sums
929       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
930       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
931                           MPI_SUM, comm2d, ierr )
932#else
933       sums = sums_l(:,:,0)
934#endif
935
936!
937!--    Final values are obtained by division by the total number of grid points
938!--    used for summation. After that store profiles.
939!--    Profiles:
940       DO  k = nzb, nzt+1
941          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
942          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
943          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
944          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
945          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
946          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
947          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
948          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
949          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
950          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
951          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
952          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
953          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
954          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
955       ENDDO
956
957!--    Upstream-parts
958       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
959!--    u* and so on
960!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
961!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
962!--    above the topography, they are being divided by ngp_2dh(sr)
963       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
964                                    ngp_2dh(sr)
965       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
966                                    ngp_2dh(sr)
967!--    eges, e*
968       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
969                                    ngp_3d(sr)
970!--    Old and new divergence
971       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
972                                    ngp_3d_inner(sr)
973
974!--    User-defined profiles
975       IF ( max_pr_user > 0 )  THEN
976          DO  k = nzb, nzt+1
977             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
978                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
979                                    ngp_2dh_s_inner(k,sr)
980          ENDDO
981       ENDIF
982!
983!--    Collect horizontal average in hom.
984!--    Compute deduced averages (e.g. total heat flux)
985       hom(:,1,3,sr)  = sums(:,3)      ! w
986       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
987       hom(:,1,9,sr)  = sums(:,9)      ! km
988       hom(:,1,10,sr) = sums(:,10)     ! kh
989       hom(:,1,11,sr) = sums(:,11)     ! l
990       hom(:,1,12,sr) = sums(:,12)     ! w"u"
991       hom(:,1,13,sr) = sums(:,13)     ! w*u*
992       hom(:,1,14,sr) = sums(:,14)     ! w"v"
993       hom(:,1,15,sr) = sums(:,15)     ! w*v*
994       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
995       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
996       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
997       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
998       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
999       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1000       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
1001                                       ! profile 24 is initial profile (sa)
1002                                       ! profiles 25-29 left empty for initial
1003                                       ! profiles
1004       hom(:,1,30,sr) = sums(:,30)     ! u*2
1005       hom(:,1,31,sr) = sums(:,31)     ! v*2
1006       hom(:,1,32,sr) = sums(:,32)     ! w*2
1007       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1008       hom(:,1,34,sr) = sums(:,34)     ! e*
1009       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1010       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1011       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1012       hom(:,1,38,sr) = sums(:,38)     ! w*3
1013       hom(:,1,39,sr) = sums(:,38) / ( sums(:,32) + 1E-20 )**1.5    ! Sw
1014       hom(:,1,40,sr) = sums(:,40)     ! p
1015       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
1016       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1017       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1018       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1019       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1020       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1021       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1022       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1023       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1024       hom(:,1,54,sr) = sums(:,54)     ! ql
1025       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1026       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
1027       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
1028       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1029       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1030       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1031       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1032       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1033       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
1034       hom(:,1,64,sr) = sums(:,64)     ! rho
1035       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1036       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1037       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1038       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1039       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
1040       hom(:,1,70,sr) = sums(:,70)     ! q*2
1041       hom(:,1,71,sr) = sums(:,71)     ! prho
1042       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
1043
1044       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
1045                                       ! upstream-parts u_x, u_y, u_z, v_x,
1046                                       ! v_y, usw. (in last but one profile)
1047       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1048                                       ! u*, w'u', w'v', t* (in last profile)
1049
1050       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1051          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1052                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1053       ENDIF
1054
1055!
1056!--    Determine the boundary layer height using two different schemes.
1057!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1058!--    first relative minimum (maximum) of the total heat flux.
1059!--    The corresponding height is assumed as the boundary layer height, if it
1060!--    is less than 1.5 times the height where the heat flux becomes negative
1061!--    (positive) for the first time.
1062       z_i(1) = 0.0
1063       first = .TRUE.
1064
1065       IF ( ocean )  THEN
1066          DO  k = nzt, nzb+1, -1
1067             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1068                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
1069                first = .FALSE.
1070                height = zw(k)
1071             ENDIF
1072             IF ( hom(k,1,18,sr) < 0.0  .AND. &
1073                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1074                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1075                IF ( zw(k) < 1.5 * height )  THEN
1076                   z_i(1) = zw(k)
1077                ELSE
1078                   z_i(1) = height
1079                ENDIF
1080                EXIT
1081             ENDIF
1082          ENDDO
1083       ELSE
1084          DO  k = nzb, nzt-1
1085             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1086                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
1087                first = .FALSE.
1088                height = zw(k)
1089             ENDIF
1090             IF ( hom(k,1,18,sr) < 0.0  .AND. & 
1091                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1092                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1093                IF ( zw(k) < 1.5 * height )  THEN
1094                   z_i(1) = zw(k)
1095                ELSE
1096                   z_i(1) = height
1097                ENDIF
1098                EXIT
1099             ENDIF
1100          ENDDO
1101       ENDIF
1102
1103!
1104!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1105!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1106!--    maximal local temperature gradient: starting from the second (the last
1107!--    but one) vertical gridpoint, the local gradient must be at least
1108!--    0.2K/100m and greater than the next four gradients.
1109!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1110!--             ocean case!
1111       z_i(2) = 0.0
1112       DO  k = nzb+1, nzt+1
1113          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1114       ENDDO
1115       dptdz_threshold = 0.2 / 100.0
1116
1117       IF ( ocean )  THEN
1118          DO  k = nzt+1, nzb+5, -1
1119             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1120                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1121                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1122                z_i(2) = zw(k-1)
1123                EXIT
1124             ENDIF
1125          ENDDO
1126       ELSE
1127          DO  k = nzb+1, nzt-3
1128             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1129                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1130                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1131                z_i(2) = zw(k-1)
1132                EXIT
1133             ENDIF
1134          ENDDO
1135       ENDIF
1136
1137       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1138       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1139
1140!
1141!--    Computation of both the characteristic vertical velocity and
1142!--    the characteristic convective boundary layer temperature.
1143!--    The horizontal average at nzb+1 is input for the average temperature.
1144       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
1145           .AND.  z_i(1) /= 0.0 )  THEN
1146          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
1147                                       hom(nzb,1,18,sr) *      &
1148                                       ABS( z_i(1) ) )**0.333333333
1149!--       so far this only works if Prandtl layer is used
1150          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
1151       ELSE
1152          hom(nzb+8,1,pr_palm,sr)  = 0.0
1153          hom(nzb+11,1,pr_palm,sr) = 0.0
1154       ENDIF
1155
1156!
1157!--    Collect the time series quantities
1158       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1159       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
1160       ts_value(3,sr) = dt_3d
1161       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1162       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
1163       ts_value(6,sr) = u_max
1164       ts_value(7,sr) = v_max
1165       ts_value(8,sr) = w_max
1166       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1167       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1168       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1169       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1170       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
1171       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1172       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1173       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1174       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1175       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
1176       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1177       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
1178       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1179
1180       IF ( ts_value(5,sr) /= 0.0 )  THEN
1181          ts_value(22,sr) = ts_value(4,sr)**2 / &
1182                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1183       ELSE
1184          ts_value(22,sr) = 10000.0
1185       ENDIF
1186
1187       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1188!
1189!--    Calculate additional statistics provided by the user interface
1190       CALL user_statistics( 'time_series', sr, 0 )
1191
1192    ENDDO    ! loop of the subregions
1193
1194!
1195!-- If required, sum up horizontal averages for subsequent time averaging
1196    IF ( do_sum )  THEN
1197       IF ( average_count_pr == 0 )  hom_sum = 0.0
1198       hom_sum = hom_sum + hom(:,1,:,:)
1199       average_count_pr = average_count_pr + 1
1200       do_sum = .FALSE.
1201    ENDIF
1202
1203!
1204!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1205!-- This flag is reset after each time step in time_integration.
1206    flow_statistics_called = .TRUE.
1207
1208    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1209
1210
1211 END SUBROUTINE flow_statistics
1212
1213
1214
Note: See TracBrowser for help on using the repository browser.