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

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

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:


Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

  • it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:


Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

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