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

Last change on this file since 388 was 388, checked in by raasch, 12 years ago

in-situ AND potential density are calculated and used in the ocean version

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