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

Last change on this file since 343 was 343, checked in by maronga, 12 years ago

adjustments for lcxt4 and ibmy, allow user 2d xy cross section output at z=nzb+1

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