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

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

New:
---
ocean version including prognostic equation for salinity and equation of state for seawater. Routine buoyancy can be used with both temperature and density.
+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux

advec_s_bc, average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, netcdf, parin, production_e, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary

New:
eqn_state_seawater, init_ocean

Changed:


inipar-parameter use_pt_reference renamed use_reference

hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile

format adjustments for the ocean version (run_control)

advec_particles, buoyancy, calc_liquid_water_content, check_parameters, diffusion_e, diffusivities, header, init_cloud_physics, modules, production_e, prognostic_equations, run_control

Errors:


Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).

Bugfix: error in boundary condition for TKE removed (advec_s_bc)

advec_s_bc, diffusion_e, prognostic_equations

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