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

Last change on this file since 696 was 696, checked in by raasch, 10 years ago

adjustments for openmp usage on ibmkisti (mrun, subjob); OpenMP-bugfixes: work_fftx removed from PRIVATE clauses in fftx_tr_xy and tr_yx_fftx (poisfft); Bugfix: Summation of Wicker-Skamarock scheme fluxes and variances for all threads (flow_statistics)

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