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

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

changes for coupling with independent precursor runs; z_i calculation with Sullivan criterion

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