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

Last change on this file since 1111 was 1111, checked in by raasch, 10 years ago

New:
---

GPU porting of pres, swap_timelevel. Adjustments of openACC directives.
Further porting of poisfft, which now runs completely on GPU without any
host/device data transfer for serial an parallel runs (but parallel runs
require data transfer before and after the MPI transpositions).
GPU-porting of tridiagonal solver:
tridiagonal routines split into extermal subroutines (instead using CONTAINS),
no distinction between parallel/non-parallel in poisfft and tridia any more,
tridia routines moved to end of file because of probable bug in PGI compiler
(otherwise "invalid device function" is indicated during runtime).
(cuda_fft_interfaces, fft_xy, flow_statistics, init_3d_model, palm, poisfft, pres, prognostic_equations, swap_timelevel, time_integration, transpose)
output of accelerator board information. (header)

optimization of tridia routines: constant elements and coefficients of tri are
stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5 to 2,
(init_grid, init_3d_model, modules, palm, poisfft)

poisfft_init is now called internally from poisfft,
(Makefile, Makefile_check, init_pegrid, poisfft, poisfft_hybrid)

CPU-time per grid point and timestep is output to CPU_MEASURES file
(cpu_statistics, modules, time_integration)

Changed:


resorting from/to array work changed, work now has 4 dimensions instead of 1 (transpose)
array diss allocated only if required (init_3d_model)

pressure boundary condition "Neumann+inhomo" removed from the code
(check_parameters, header, poisfft, poisfft_hybrid, pres)

Errors:


bugfix: dependency added for cuda_fft_interfaces (Makefile)
bugfix: CUDA fft plans adjusted for domain decomposition (before they always
used total domain) (fft_xy)

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