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

Last change on this file since 1182 was 1182, checked in by raasch, 8 years ago

last commit documented, rc-file for example run updated

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