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

Last change on this file since 668 was 668, checked in by suehring, 11 years ago

last commit documented

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