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

Last change on this file since 743 was 743, checked in by suehring, 10 years ago

Calculation of turbulent fluxes with WS-scheme only for the whole model domain, not for user-defined subregions

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