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

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

last commit documented

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