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

Last change on this file since 224 was 220, checked in by raasch, 16 years ago

some small bugfixes in user_module, user_read_restart_data, read_3d_binary, flow_statistics and mrun

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