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

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

Last commit documented.

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