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

Last change on this file since 1015 was 1008, checked in by franke, 12 years ago

last commit documented

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