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

Last change on this file since 159 was 142, checked in by letzel, 16 years ago

Bugfix in plant_canopy_model: remove IF statement in plant_canopy_model_ij
Bugfix in flow_statistics: divide sums(k,8) (e) and sums(k,34) (e*) by
ngp_2dh_s_inner(k,sr) (like other scalars)

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