Home

Context Navigation

flow_statistics.f90 @ 519

Last change on this file since 519 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

Rev	Line
[1]	1	SUBROUTINE flow_statistics
	2
	3	!------------------------------------------------------------------------------!
[254]	4	! Current revisions:
[1]	5	! -----------------
[392]	6	!
	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: flow_statistics.f90 392 2009-09-24 10:39:14Z raasch $
	11	!
	12	! 388 2009-09-23 09:40:33Z raasch
[388]	13	! Vertical profiles of potential density and hydrostatic pressure are
	14	! calculated.
[343]	15	! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
	16	! end.
[291]	17	! Temperature gradient criterion for estimating the boundary layer height
	18	! replaced by the gradient criterion of Sullivan et al. (1998).
[254]	19	! Output of messages replaced by message handling routine.
[1]	20	!
[198]	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	!
[139]	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	!
[110]	33	! 106 2007-08-16 14:30:26Z raasch
	34	! Prescribed momentum fluxes at the top surface are used,
	35	! profiles for wp and w"e are calculated
	36	!
[98]	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	!
[90]	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	!
[83]	47	! 82 2007-04-16 15:40:52Z raasch
	48	! Cpp-directive lcmuk changed to intel_openmp_bug
	49	!
[77]	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	!
[39]	55	! 19 2007-02-23 04:53:48Z raasch
[77]	56	! fluxes at top modified (tswst, qswst)
[39]	57	!
[3]	58	! RCS Log replace by Id keyword, revision history cleaned up
	59	!
[1]	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	!
[132]	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.
[1]	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
[291]	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)
[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
[254]	105
[274]	106	message_string = 'flow_statistics is called two times within one ' // &
	107	'timestep'
[254]	108	CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
	109
[1]	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
[87]	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
[1]	128
[305]	129	!
[1]	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 )
[82]	135	#if defined( __intel_openmp_bug )
[1]	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
[132]	144	DO k = nzb_s_inner(j,i), nzt+1
[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	!
[96]	153	!-- Horizontally averaged profile of salinity
	154	IF ( ocean ) THEN
	155	!$OMP DO
	156	DO i = nxl, nxr
	157	DO j = nys, nyn
[132]	158	DO k = nzb_s_inner(j,i), nzt+1
[96]	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	!
[1]	167	!-- Horizontally averaged profiles of virtual potential temperature,
	168	!-- total water content, specific humidity and liquid water potential
	169	!-- temperature
[75]	170	IF ( humidity ) THEN
[1]	171	!$OMP DO
	172	DO i = nxl, nxr
	173	DO j = nys, nyn
[132]	174	DO k = nzb_s_inner(j,i), nzt+1
[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
[132]	186	DO k = nzb_s_inner(j,i), nzt+1
[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_cppt_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
[132]	204	DO k = nzb_s_inner(j,i), nzt+1
[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)
[96]	219	IF ( ocean ) THEN
	220	sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
	221	ENDIF
[75]	222	IF ( humidity ) THEN
[1]	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 )
[96]	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
[75]	249	IF ( humidity ) THEN
[1]	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)
[96]	270	IF ( ocean ) sums(:,23) = sums_l(:,23,0)
[75]	271	IF ( humidity ) THEN
[1]	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.
[132]	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)
[1]	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	!
[96]	293	!-- Salinity
	294	IF ( ocean ) THEN
[132]	295	sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
[96]	296	hom(:,1,23,sr) = sums(:,23) ! sa
	297	ENDIF
	298
	299	!
[1]	300	!-- Humidity and cloud parameters
[75]	301	IF ( humidity ) THEN
[132]	302	sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
	303	sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
[1]	304	hom(:,1,44,sr) = sums(:,44) ! vpt
	305	hom(:,1,41,sr) = sums(:,41) ! qv (q)
	306	IF ( cloud_physics ) THEN
[132]	307	sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
	308	sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
[1]	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
[132]	316	IF ( passive_scalar ) hom(:,1,41,sr) = sums(:,41) / &
	317	ngp_2dh_s_inner(:,sr) ! s (q)
[1]	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.
[132]	323	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
[1]	324	!-- ---- speaking the following k-loop would have to be split up and
	325	!-- rearranged according to the staggered grid.
[132]	326	!-- However, this implies no error since staggered velocity components
	327	!-- are zero at the walls and inside buildings.
[1]	328	tn = 0
[82]	329	#if defined( __intel_openmp_bug )
[1]	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
[132]	342	DO k = nzb_s_inner(j,i), nzt+1
[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.
[87]	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
[1]	384	!
	385	!-- 2D-arrays (being collected in the last column of sums_l)
[87]	386	sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + &
[1]	387	us(j,i) * rmask(j,i,sr)
[87]	388	sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
[1]	389	usws(j,i) * rmask(j,i,sr)
[87]	390	sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
[1]	391	vsws(j,i) * rmask(j,i,sr)
[87]	392	sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
[1]	393	ts(j,i) * rmask(j,i,sr)
[197]	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
[1]	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)
[132]	409	!-- NOTE: for simplicity, nzb_diff_s_inner is used below, although
[1]	410	!-- ---- strictly speaking the following k-loop would have to be
	411	!-- split up according to the staggered grid.
[132]	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
[1]	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	!
[96]	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	!
[1]	450	!-- Buoyancy flux, water flux (humidity flux) w"q"
[75]	451	IF ( humidity ) THEN
[1]	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"
[96]	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
[75]	497	IF ( humidity ) THEN
[1]	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	!
[19]	519	!-- Subgridscale fluxes at the top surface
	520	IF ( use_top_fluxes ) THEN
[102]	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"
[19]	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"
[96]	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
[75]	535	IF ( humidity ) THEN
[19]	536	sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
[388]	537	qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
[19]	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) + &
[388]	552	qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
[19]	553	ENDIF
	554	ENDIF
	555
	556	!
[1]	557	!-- Resolved fluxes (can be computed for all horizontal points)
[132]	558	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
[1]	559	!-- ---- speaking the following k-loop would have to be split up and
	560	!-- rearranged according to the staggered grid.
[132]	561	DO k = nzb_s_inner(j,i), nzt
[1]	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 wu
	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 wv
	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 wpt
	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	!
[96]	595	!-- Salinity flux and density (density does not belong to here,
[97]	596	!-- but so far there is no other suitable place to calculate)
[96]	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)
[388]	604	sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
	605	rmask(j,i,sr)
[96]	606	ENDIF
	607
	608	!
[1]	609	!-- Buoyancy flux, water flux, humidity flux and liquid water
	610	!-- content
[75]	611	IF ( humidity ) THEN
[1]	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 we
	642	sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 * &
	643	( ust2 + vst2 + w(k,j,i)**2 )&
	644	* rmask(j,i,sr)
	645
	646	ENDDO
	647	ENDDO
	648	ENDDO
	649
	650	!
[97]	651	!-- Density at top follows Neumann condition
[388]	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
[97]	656
	657	!
[1]	658	!-- Divergence of vertical flux of resolved scale energy and pressure
[106]	659	!-- fluctuations as well as flux of pressure fluctuation itself (68).
	660	!-- First calculate the products, then the divergence.
[1]	661	!-- Calculation is time consuming. Do it only, if profiles shall be plotted.
[106]	662	IF ( hom(nzb+1,2,55,0) /= 0.0 .OR. hom(nzb+1,2,68,0) /= 0.0 ) THEN
[1]	663
	664	sums_ll = 0.0 ! local array
	665
	666	!$OMP DO
	667	DO i = nxl, nxr
	668	DO j = nys, nyn
[132]	669	DO k = nzb_s_inner(j,i)+1, nzt
[1]	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
[132]	691	DO k = nzb_s_inner(j,i)+1, nzt
[1]	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)
[106]	694	sums_l(k,68,tn) = sums_ll(k,2)
[1]	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)
[106]	698	sums_l(nzb,68,tn) = 0.0 ! because w* = 0 at nzb
[1]	699
	700	ENDIF
	701
	702	!
[106]	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
[1]	705
	706	!$OMP DO
	707	DO i = nxl, nxr
	708	DO j = nys, nyn
[132]	709	DO k = nzb_s_inner(j,i)+1, nzt
[1]	710
[106]	711	sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * ( &
[1]	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) &
[106]	714	) * ddzw(k)
[1]	715
[106]	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
[1]	720	ENDDO
	721	ENDDO
	722	ENDDO
	723	sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
[106]	724	sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
[1]	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
[132]	736	DO k = nzb_s_inner(j,i)+1, nzt
[1]	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 upt, vpt
	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)
[97]	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
[1]	770
	771	ENDIF
[87]	772
	773	!
	774	!-- Calculate the user-defined profiles
	775	CALL user_statistics( 'profiles', sr, tn )
[1]	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)
[87]	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
[1]	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
[132]	808	sums(k,3) = sums(k,3) / ngp_2dh(sr)
[142]	809	sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr)
[132]	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)
[142]	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)
[132]	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)
[1]	822	ENDDO
	823	!-- Upstream-parts
[87]	824	sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
[1]	825	!-- u* and so on
[87]	826	!-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
[1]	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)
[87]	829	sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / &
[1]	830	ngp_2dh(sr)
[197]	831	sums(nzb+12,pr_palm) = sums(nzb+12,pr_palm) / & ! qs
	832	ngp_2dh(sr)
[1]	833	!-- eges, e*
[87]	834	sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / &
[132]	835	ngp_3d(sr)
[1]	836	!-- Old and new divergence
[87]	837	sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
[1]	838	ngp_3d_inner(sr)
	839
[87]	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) / &
[132]	845	ngp_2dh_s_inner(k,sr)
[87]	846	ENDDO
	847	ENDIF
	848
[1]	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) ! wu
	859	hom(:,1,14,sr) = sums(:,14) ! w"v"
	860	hom(:,1,15,sr) = sums(:,15) ! wv
	861	hom(:,1,16,sr) = sums(:,16) ! w"pt"
	862	hom(:,1,17,sr) = sums(:,17) ! wpt
	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) ! wptBC
	867	hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC
[96]	868	! profile 24 is initial profile (sa)
	869	! profiles 25-29 left empty for initial
[1]	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) ! w2pt
	877	hom(:,1,36,sr) = sums(:,36) ! wpt2
	878	hom(:,1,37,sr) = sums(:,37) ! we
	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) ! wvpt
	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) ! wq (wqv)
	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) ! wqv
	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) ! wuu*/dz
	893	hom(:,1,56,sr) = sums(:,56) ! wp/dz
[106]	894	hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho )/dz
[1]	895	hom(:,1,58,sr) = sums(:,58) ! u"pt"
	896	hom(:,1,59,sr) = sums(:,59) ! upt
	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) ! vpt
	900	hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t
[96]	901	hom(:,1,64,sr) = sums(:,64) ! rho
	902	hom(:,1,65,sr) = sums(:,65) ! w"sa"
	903	hom(:,1,66,sr) = sums(:,66) ! wsa
	904	hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa
[106]	905	hom(:,1,68,sr) = sums(:,68) ! wp
	906	hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho
[197]	907	hom(:,1,70,sr) = sums(:,70) ! q*2
[388]	908	hom(:,1,71,sr) = sums(:,71) ! prho
	909	hom(:,1,72,sr) = hyp * 1E-4 ! hyp in kPa
[1]	910
[87]	911	hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
[1]	912	! upstream-parts u_x, u_y, u_z, v_x,
	913	! v_y, usw. (in last but one profile)
[87]	914	hom(:,1,pr_palm,sr) = sums(:,pr_palm)
[1]	915	! u, w'u', w'v', t (in last profile)
	916
[87]	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
[1]	922	!
	923	!-- Determine the boundary layer height using two different schemes.
[94]	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.
[1]	929	z_i(1) = 0.0
	930	first = .TRUE.
[97]	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
[94]	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)
[1]	952	ENDIF
[94]	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
[97]	963	ENDIF
[1]	964
	965	!
[291]	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!
[1]	973	z_i(2) = 0.0
[291]	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
[97]	979	IF ( ocean ) THEN
[291]	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)
[97]	985	EXIT
	986	ENDIF
	987	ENDDO
	988	ELSE
[291]	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)
[97]	994	EXIT
	995	ENDIF
	996	ENDDO
	997	ENDIF
[1]	998
[87]	999	hom(nzb+6,1,pr_palm,sr) = z_i(1)
	1000	hom(nzb+7,1,pr_palm,sr) = z_i(2)
[1]	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
[87]	1007	hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * &
[94]	1008	hom(nzb,1,18,sr) * &
	1009	ABS( z_i(1) ) )**0.333333333
[1]	1010	!-- so far this only works if Prandtl layer is used
[87]	1011	hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
[1]	1012	ELSE
[87]	1013	hom(nzb+8,1,pr_palm,sr) = 0.0
	1014	hom(nzb+11,1,pr_palm,sr) = 0.0
[1]	1015	ENDIF
	1016
[48]	1017	!
	1018	!-- Collect the time series quantities
[87]	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*
[48]	1021	ts_value(3,sr) = dt_3d
[87]	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*
[48]	1024	ts_value(6,sr) = u_max
	1025	ts_value(7,sr) = v_max
	1026	ts_value(8,sr) = w_max
[87]	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*
[48]	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)
[197]	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
[343]	1039	ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0
[197]	1040
[48]	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
[1]	1047
[343]	1048	ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q*
[1]	1049	!
[48]	1050	!-- Calculate additional statistics provided by the user interface
[87]	1051	CALL user_statistics( 'time_series', sr, 0 )
[1]	1052
[48]	1053	ENDDO ! loop of the subregions
	1054
[1]	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.

Context Navigation

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

Download in other formats: