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

Last change on this file since 622 was 622, checked in by raasch, 11 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:


Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:


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