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

Last change on this file since 254 was 254, checked in by heinze, 12 years ago

Output of messages replaced by message handling routine.

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