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

Last change on this file since 550 was 550, checked in by maronga, 14 years ago

bugfix: calculation of subgrid-scale flux at nzt for use_top_fluxes=.T.

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