Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 416 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

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