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

Last change on this file since 106 was 106, checked in by raasch, 14 years ago

preliminary update of bugfixes and extensions for non-cyclic BCs

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