source: palm/tags/release-3.8/SOURCE/flow_statistics.f90 @ 3884

Last change on this file since 3884 was 710, checked in by raasch, 13 years ago

last commit documented

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