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

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

Bugfix in calculation of vertical velocity skewness.

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