Home

Context Navigation

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

Last change on this file since 2385 was 2320, checked in by suehring, 7 years ago
large-scale forcing and nudging modularized
Property svn:keywords set to `Id`
File size: 95.0 KB

Line
1	!> @file flow_statistics.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: flow_statistics.f90 2320 2017-07-21 12:47:43Z raasch $
27	! Modularize large-scale forcing and nudging
28	!
29	! 2296 2017-06-28 07:53:56Z maronga
30	! Enabled output of radiation quantities for radiation_scheme = 'constant'
31	!
32	! 2292 2017-06-20 09:51:42Z schwenkel
33	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
34	! includes two more prognostic equations for cloud drop concentration (nc)
35	! and cloud water content (qc).
36	!
37	! 2270 2017-06-09 12:18:47Z maronga
38	! Revised numbering (removed 2 timeseries)
39	!
40	! 2252 2017-06-07 09:35:37Z knoop
41	! perturbation pressure now depending on flux_output_mode
42	!
43	! 2233 2017-05-30 18:08:54Z suehring
44	!
45	! 2232 2017-05-30 17:47:52Z suehring
46	! Adjustments to new topography and surface concept
47	!
48	! 2118 2017-01-17 16:38:49Z raasch
49	! OpenACC version of subroutine removed
50	!
51	! 2073 2016-11-30 14:34:05Z raasch
52	! openmp bugfix: large scale forcing calculations cannot be executed thread
53	! parallel
54	!
55	! 2037 2016-10-26 11:15:40Z knoop
56	! Anelastic approximation implemented
57	!
58	! 2031 2016-10-21 15:11:58Z knoop
59	! renamed variable rho to rho_ocean
60	!
61	! 2026 2016-10-18 10:27:02Z suehring
62	! Bugfix, enable output of s*2.
63	! Change, calculation of domain-averaged perturbation energy.
64	! Some formatting adjustments.
65	!
66	! 2000 2016-08-20 18:09:15Z knoop
67	! Forced header and separation lines into 80 columns
68	!
69	! 1976 2016-07-27 13:28:04Z maronga
70	! Removed some unneeded __rrtmg preprocessor directives
71	!
72	! 1960 2016-07-12 16:34:24Z suehring
73	! Separate humidity and passive scalar
74	!
75	! 1918 2016-05-27 14:35:57Z raasch
76	! in case of Wicker-Skamarock scheme, calculate disturbance kinetic energy here,
77	! if flow_statistics is called before the first initial time step
78	!
79	! 1853 2016-04-11 09:00:35Z maronga
80	! Adjusted for use with radiation_scheme = constant
81	!
82	! 1849 2016-04-08 11:33:18Z hoffmann
83	! prr moved to arrays_3d
84	!
85	! 1822 2016-04-07 07:49:42Z hoffmann
86	! Output of bulk microphysics simplified.
87	!
88	! 1815 2016-04-06 13:49:59Z raasch
89	! cpp-directives for intel openmp bug removed
90	!
91	! 1783 2016-03-06 18:36:17Z raasch
92	! +module netcdf_interface
93	!
94	! 1747 2016-02-08 12:25:53Z raasch
95	! small bugfixes for accelerator version
96	!
97	! 1738 2015-12-18 13:56:05Z raasch
98	! bugfixes for calculations in statistical regions which do not contain grid
99	! points in the lowest vertical levels, mean surface level height considered
100	! in the calculation of the characteristic vertical velocity,
101	! old upstream parts removed
102	!
103	! 1709 2015-11-04 14:47:01Z maronga
104	! Updated output of Obukhov length
105	!
106	! 1691 2015-10-26 16:17:44Z maronga
107	! Revised calculation of Obukhov length. Added output of radiative heating >
108	! rates for RRTMG.
109	!
110	! 1682 2015-10-07 23:56:08Z knoop
111	! Code annotations made doxygen readable
112	!
113	! 1658 2015-09-18 10:52:53Z raasch
114	! bugfix: temporary reduction variables in the openacc branch are now
115	! initialized to zero
116	!
117	! 1654 2015-09-17 09:20:17Z raasch
118	! FORTRAN bugfix of r1652
119	!
120	! 1652 2015-09-17 08:12:24Z raasch
121	! bugfix in calculation of energy production by turbulent transport of TKE
122	!
123	! 1593 2015-05-16 13:58:02Z raasch
124	! FORTRAN errors removed from openacc branch
125	!
126	! 1585 2015-04-30 07:05:52Z maronga
127	! Added output of timeseries and profiles for RRTMG
128	!
129	! 1571 2015-03-12 16:12:49Z maronga
130	! Bugfix: output of rad_net and rad_sw_in
131	!
132	! 1567 2015-03-10 17:57:55Z suehring
133	! Reverse modifications made for monotonic limiter.
134	!
135	! 1557 2015-03-05 16:43:04Z suehring
136	! Adjustments for monotonic limiter
137	!
138	! 1555 2015-03-04 17:44:27Z maronga
139	! Added output of r_a and r_s.
140	!
141	! 1551 2015-03-03 14:18:16Z maronga
142	! Added suppport for land surface model and radiation model output.
143	!
144	! 1498 2014-12-03 14:09:51Z suehring
145	! Comments added
146	!
147	! 1482 2014-10-18 12:34:45Z raasch
148	! missing ngp_sums_ls added in accelerator version
149	!
150	! 1450 2014-08-21 07:31:51Z heinze
151	! bugfix: calculate fac only for simulated_time >= 0.0
152	!
153	! 1396 2014-05-06 13:37:41Z raasch
154	! bugfix: "copyin" replaced by "update device" in openacc-branch
155	!
156	! 1386 2014-05-05 13:55:30Z boeske
157	! bugfix: simulated time before the last timestep is needed for the correct
158	! calculation of the profiles of large scale forcing tendencies
159	!
160	! 1382 2014-04-30 12:15:41Z boeske
161	! Renamed variables which store large scale forcing tendencies
162	! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
163	! q_lsa -> td_lsa_q, q_subs -> td_sub_q,
164	! added Neumann boundary conditions for profile data output of large scale
165	! advection and subsidence terms at nzt+1
166	!
167	! 1374 2014-04-25 12:55:07Z raasch
168	! bugfix: syntax errors removed from openacc-branch
169	! missing variables added to ONLY-lists
170	!
171	! 1365 2014-04-22 15:03:56Z boeske
172	! Output of large scale advection, large scale subsidence and nudging tendencies
173	! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
174	!
175	! 1353 2014-04-08 15:21:23Z heinze
176	! REAL constants provided with KIND-attribute
177	!
178	! 1322 2014-03-20 16:38:49Z raasch
179	! REAL constants defined as wp-kind
180	!
181	! 1320 2014-03-20 08:40:49Z raasch
182	! ONLY-attribute added to USE-statements,
183	! kind-parameters added to all INTEGER and REAL declaration statements,
184	! kinds are defined in new module kinds,
185	! revision history before 2012 removed,
186	! comment fields (!:) to be used for variable explanations added to
187	! all variable declaration statements
188	!
189	! 1299 2014-03-06 13:15:21Z heinze
190	! Output of large scale vertical velocity w_subs
191	!
192	! 1257 2013-11-08 15:18:40Z raasch
193	! openacc "end parallel" replaced by "end parallel loop"
194	!
195	! 1241 2013-10-30 11:36:58Z heinze
196	! Output of ug and vg
197	!
198	! 1221 2013-09-10 08:59:13Z raasch
199	! ported for openACC in separate #else branch
200	!
201	! 1179 2013-06-14 05:57:58Z raasch
202	! comment for profile 77 added
203	!
204	! 1115 2013-03-26 18:16:16Z hoffmann
205	! ql is calculated by calc_liquid_water_content
206	!
207	! 1111 2013-03-08 23:54:10Z raasch
208	! openACC directive added
209	!
210	! 1053 2012-11-13 17:11:03Z hoffmann
211	! additions for two-moment cloud physics scheme:
212	! +nr, qr, qc, prr
213	!
214	! 1036 2012-10-22 13:43:42Z raasch
215	! code put under GPL (PALM 3.9)
216	!
217	! 1007 2012-09-19 14:30:36Z franke
218	! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
219	! turbulent fluxes of WS-scheme
220	! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
221	! droplets at nzb and nzt added
222	!
223	! 801 2012-01-10 17:30:36Z suehring
224	! Calculation of turbulent fluxes in advec_ws is now thread-safe.
225	!
226	! Revision 1.1 1997/08/11 06:15:17 raasch
227	! Initial revision
228	!
229	!
230	! Description:
231	! ------------
232	!> Compute average profiles and further average flow quantities for the different
233	!> user-defined (sub-)regions. The region indexed 0 is the total model domain.
234	!>
235	!> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
236	!> lower vertical index for k-loops for all variables, although strictly
237	!> speaking the k-loops would have to be split up according to the staggered
238	!> grid. However, this implies no error since staggered velocity components
239	!> are zero at the walls and inside buildings.
240	!------------------------------------------------------------------------------!
241	SUBROUTINE flow_statistics
242
243
244	USE arrays_3d, &
245	ONLY: ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh, &
246	momentumflux_output_conversion, nc, nr, p, prho, prr, pt, q, &
247	qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s, &
248	sa, u, ug, v, vg, vpt, w, w_subs, waterflux_output_conversion, zw
249
250	USE cloud_parameters, &
251	ONLY: l_d_cp, pt_d_t
252
253	USE control_parameters, &
254	ONLY: average_count_pr, cloud_droplets, cloud_physics, do_sum, &
255	dt_3d, g, humidity, kappa, land_surface, large_scale_forcing, &
256	large_scale_subsidence, max_pr_user, message_string, neutral, &
257	microphysics_morrison, microphysics_seifert, ocean, &
258	passive_scalar, simulated_time, use_subsidence_tendencies, &
259	use_surface_fluxes, use_top_fluxes, ws_scheme_mom, &
260	ws_scheme_sca
261
262	USE cpulog, &
263	ONLY: cpu_log, log_point
264
265	USE grid_variables, &
266	ONLY: ddx, ddy
267
268	USE indices, &
269	ONLY: ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, &
270	ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
271
272	USE kinds
273
274	USE land_surface_model_mod, &
275	ONLY: m_soil_h, nzb_soil, nzt_soil, t_soil_h
276
277	USE lsf_nudging_mod, &
278	ONLY: td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert
279
280	USE netcdf_interface, &
281	ONLY: dots_rad, dots_soil
282
283	USE pegrid
284
285	USE radiation_model_mod, &
286	ONLY: radiation, radiation_scheme, rad_net, &
287	rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr, &
288	rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr
289
290	#if defined ( __rrtmg )
291	USE radiation_model_mod, &
292	ONLY: rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
293	#endif
294
295	USE statistics
296
297	USE surface_mod, &
298	ONLY : surf_def_h, surf_lsm_h, surf_usm_h
299
300
301	IMPLICIT NONE
302
303	INTEGER(iwp) :: i !<
304	INTEGER(iwp) :: j !<
305	INTEGER(iwp) :: k !<
306	INTEGER(iwp) :: ki !<
307	INTEGER(iwp) :: k_surface_level !<
308	INTEGER(iwp) :: m !< loop variable over all horizontal wall elements
309	INTEGER(iwp) :: l !< loop variable over surface facing -- up- or downward-facing
310	INTEGER(iwp) :: nt !<
311	INTEGER(iwp) :: omp_get_thread_num !<
312	INTEGER(iwp) :: sr !<
313	INTEGER(iwp) :: tn !<
314
315	LOGICAL :: first !<
316
317	REAL(wp) :: dptdz_threshold !<
318	REAL(wp) :: fac !<
319	REAL(wp) :: flag !<
320	REAL(wp) :: height !<
321	REAL(wp) :: pts !<
322	REAL(wp) :: sums_l_eper !<
323	REAL(wp) :: sums_l_etot !<
324	REAL(wp) :: ust !<
325	REAL(wp) :: ust2 !<
326	REAL(wp) :: u2 !<
327	REAL(wp) :: vst !<
328	REAL(wp) :: vst2 !<
329	REAL(wp) :: v2 !<
330	REAL(wp) :: w2 !<
331	REAL(wp) :: z_i(2) !<
332
333	REAL(wp) :: dptdz(nzb+1:nzt+1) !<
334	REAL(wp) :: sums_ll(nzb:nzt+1,2) !<
335
336	CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
337
338
339	!
340	!-- To be on the safe side, check whether flow_statistics has already been
341	!-- called once after the current time step
342	IF ( flow_statistics_called ) THEN
343
344	message_string = 'flow_statistics is called two times within one ' // &
345	'timestep'
346	CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
347
348	ENDIF
349
350	!
351	!-- Compute statistics for each (sub-)region
352	DO sr = 0, statistic_regions
353
354	!
355	!-- Initialize (local) summation array
356	sums_l = 0.0_wp
357
358	!
359	!-- Store sums that have been computed in other subroutines in summation
360	!-- array
361	sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities
362	!-- WARNING: next line still has to be adjusted for OpenMP
363	sums_l(:,21,0) = sums_wsts_bc_l(:,sr) * &
364	heatflux_output_conversion ! heat flux from advec_s_bc
365	sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres
366	sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres
367
368	!
369	!-- When calcuating horizontally-averaged total (resolved- plus subgrid-
370	!-- scale) vertical fluxes and velocity variances by using commonly-
371	!-- applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
372	!-- in combination with the 5th order advection scheme, pronounced
373	!-- artificial kinks could be observed in the vertical profiles near the
374	!-- surface. Please note: these kinks were not related to the model truth,
375	!-- i.e. these kinks are just related to an evaluation problem.
376	!-- In order avoid these kinks, vertical fluxes and horizontal as well
377	!-- vertical velocity variances are calculated directly within the advection
378	!-- routines, according to the numerical discretization, to evaluate the
379	!-- statistical quantities as they will appear within the prognostic
380	!-- equations.
381	!-- Copy the turbulent quantities, evaluated in the advection routines to
382	!-- the local array sums_l() for further computations.
383	IF ( ws_scheme_mom .AND. sr == 0 ) THEN
384
385	!
386	!-- According to the Neumann bc for the horizontal velocity components,
387	!-- the corresponding fluxes has to satisfiy the same bc.
388	IF ( ocean ) THEN
389	sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
390	sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
391	ENDIF
392
393	DO i = 0, threads_per_task-1
394	!
395	!-- Swap the turbulent quantities evaluated in advec_ws.
396	sums_l(:,13,i) = sums_wsus_ws_l(:,i) &
397	* momentumflux_output_conversion ! wu
398	sums_l(:,15,i) = sums_wsvs_ws_l(:,i) &
399	* momentumflux_output_conversion ! wv
400	sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2
401	sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2
402	sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2
403	sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp * &
404	( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + &
405	sums_ws2_ws_l(:,i) ) ! e*
406	ENDDO
407
408	ENDIF
409
410	IF ( ws_scheme_sca .AND. sr == 0 ) THEN
411
412	DO i = 0, threads_per_task-1
413	sums_l(:,17,i) = sums_wspts_ws_l(:,i) &
414	* heatflux_output_conversion ! wpt
415	IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! wsa
416	IF ( humidity ) sums_l(:,49,i) = sums_wsqs_ws_l(:,i) &
417	* waterflux_output_conversion ! wq
418	IF ( passive_scalar ) sums_l(:,114,i) = sums_wsss_ws_l(:,i) ! ws
419	ENDDO
420
421	ENDIF
422	!
423	!-- Horizontally averaged profiles of horizontal velocities and temperature.
424	!-- They must have been computed before, because they are already required
425	!-- for other horizontal averages.
426	tn = 0
427	!$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
428	!$ tn = omp_get_thread_num()
429	!$OMP DO
430	DO i = nxl, nxr
431	DO j = nys, nyn
432	DO k = nzb, nzt+1
433	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
434	sums_l(k,1,tn) = sums_l(k,1,tn) + u(k,j,i) * rmask(j,i,sr) &
435	* flag
436	sums_l(k,2,tn) = sums_l(k,2,tn) + v(k,j,i) * rmask(j,i,sr) &
437	* flag
438	sums_l(k,4,tn) = sums_l(k,4,tn) + pt(k,j,i) * rmask(j,i,sr) &
439	* flag
440	ENDDO
441	ENDDO
442	ENDDO
443
444	!
445	!-- Horizontally averaged profile of salinity
446	IF ( ocean ) THEN
447	!$OMP DO
448	DO i = nxl, nxr
449	DO j = nys, nyn
450	DO k = nzb, nzt+1
451	sums_l(k,23,tn) = sums_l(k,23,tn) + sa(k,j,i) &
452	* rmask(j,i,sr) &
453	* MERGE( 1.0_wp, 0.0_wp, &
454	BTEST( wall_flags_0(k,j,i), 22 ) )
455	ENDDO
456	ENDDO
457	ENDDO
458	ENDIF
459
460	!
461	!-- Horizontally averaged profiles of virtual potential temperature,
462	!-- total water content, specific humidity and liquid water potential
463	!-- temperature
464	IF ( humidity ) THEN
465	!$OMP DO
466	DO i = nxl, nxr
467	DO j = nys, nyn
468	DO k = nzb, nzt+1
469	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
470	sums_l(k,44,tn) = sums_l(k,44,tn) + &
471	vpt(k,j,i) * rmask(j,i,sr) * flag
472	sums_l(k,41,tn) = sums_l(k,41,tn) + &
473	q(k,j,i) * rmask(j,i,sr) * flag
474	ENDDO
475	ENDDO
476	ENDDO
477	IF ( cloud_physics ) THEN
478	!$OMP DO
479	DO i = nxl, nxr
480	DO j = nys, nyn
481	DO k = nzb, nzt+1
482	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
483	sums_l(k,42,tn) = sums_l(k,42,tn) + &
484	( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) &
485	* flag
486	sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
487	pt(k,j,i) + l_d_cppt_d_t(k) ql(k,j,i) &
488	) * rmask(j,i,sr) &
489	* flag
490	ENDDO
491	ENDDO
492	ENDDO
493	ENDIF
494	ENDIF
495
496	!
497	!-- Horizontally averaged profiles of passive scalar
498	IF ( passive_scalar ) THEN
499	!$OMP DO
500	DO i = nxl, nxr
501	DO j = nys, nyn
502	DO k = nzb, nzt+1
503	sums_l(k,115,tn) = sums_l(k,115,tn) + s(k,j,i) &
504	* rmask(j,i,sr) &
505	* MERGE( 1.0_wp, 0.0_wp, &
506	BTEST( wall_flags_0(k,j,i), 22 ) )
507	ENDDO
508	ENDDO
509	ENDDO
510	ENDIF
511	!$OMP END PARALLEL
512	!
513	!-- Summation of thread sums
514	IF ( threads_per_task > 1 ) THEN
515	DO i = 1, threads_per_task-1
516	sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
517	sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
518	sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
519	IF ( ocean ) THEN
520	sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
521	ENDIF
522	IF ( humidity ) THEN
523	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
524	sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
525	IF ( cloud_physics ) THEN
526	sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
527	sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
528	ENDIF
529	ENDIF
530	IF ( passive_scalar ) THEN
531	sums_l(:,115,0) = sums_l(:,115,0) + sums_l(:,115,i)
532	ENDIF
533	ENDDO
534	ENDIF
535
536	#if defined( __parallel )
537	!
538	!-- Compute total sum from local sums
539	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
540	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
541	MPI_SUM, comm2d, ierr )
542	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
543	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
544	MPI_SUM, comm2d, ierr )
545	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
546	CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
547	MPI_SUM, comm2d, ierr )
548	IF ( ocean ) THEN
549	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
550	CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
551	MPI_REAL, MPI_SUM, comm2d, ierr )
552	ENDIF
553	IF ( humidity ) THEN
554	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
555	CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
556	MPI_REAL, MPI_SUM, comm2d, ierr )
557	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
558	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
559	MPI_REAL, MPI_SUM, comm2d, ierr )
560	IF ( cloud_physics ) THEN
561	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
562	CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
563	MPI_REAL, MPI_SUM, comm2d, ierr )
564	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
565	CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
566	MPI_REAL, MPI_SUM, comm2d, ierr )
567	ENDIF
568	ENDIF
569
570	IF ( passive_scalar ) THEN
571	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
572	CALL MPI_ALLREDUCE( sums_l(nzb,115,0), sums(nzb,115), nzt+2-nzb, &
573	MPI_REAL, MPI_SUM, comm2d, ierr )
574	ENDIF
575	#else
576	sums(:,1) = sums_l(:,1,0)
577	sums(:,2) = sums_l(:,2,0)
578	sums(:,4) = sums_l(:,4,0)
579	IF ( ocean ) sums(:,23) = sums_l(:,23,0)
580	IF ( humidity ) THEN
581	sums(:,44) = sums_l(:,44,0)
582	sums(:,41) = sums_l(:,41,0)
583	IF ( cloud_physics ) THEN
584	sums(:,42) = sums_l(:,42,0)
585	sums(:,43) = sums_l(:,43,0)
586	ENDIF
587	ENDIF
588	IF ( passive_scalar ) sums(:,115) = sums_l(:,115,0)
589	#endif
590
591	!
592	!-- Final values are obtained by division by the total number of grid points
593	!-- used for summation. After that store profiles.
594	sums(:,1) = sums(:,1) / ngp_2dh(sr)
595	sums(:,2) = sums(:,2) / ngp_2dh(sr)
596	sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
597	hom(:,1,1,sr) = sums(:,1) ! u
598	hom(:,1,2,sr) = sums(:,2) ! v
599	hom(:,1,4,sr) = sums(:,4) ! pt
600
601
602	!
603	!-- Salinity
604	IF ( ocean ) THEN
605	sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
606	hom(:,1,23,sr) = sums(:,23) ! sa
607	ENDIF
608
609	!
610	!-- Humidity and cloud parameters
611	IF ( humidity ) THEN
612	sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
613	sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
614	hom(:,1,44,sr) = sums(:,44) ! vpt
615	hom(:,1,41,sr) = sums(:,41) ! qv (q)
616	IF ( cloud_physics ) THEN
617	sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
618	sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
619	hom(:,1,42,sr) = sums(:,42) ! qv
620	hom(:,1,43,sr) = sums(:,43) ! pt
621	ENDIF
622	ENDIF
623
624	!
625	!-- Passive scalar
626	IF ( passive_scalar ) hom(:,1,115,sr) = sums(:,115) / &
627	ngp_2dh_s_inner(:,sr) ! s
628
629	!
630	!-- Horizontally averaged profiles of the remaining prognostic variables,
631	!-- variances, the total and the perturbation energy (single values in last
632	!-- column of sums_l) and some diagnostic quantities.
633	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
634	!-- ---- speaking the following k-loop would have to be split up and
635	!-- rearranged according to the staggered grid.
636	!-- However, this implies no error since staggered velocity components
637	!-- are zero at the walls and inside buildings.
638	tn = 0
639	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, &
640	!$OMP sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, &
641	!$OMP w2, flag, m, ki, l )
642	!$ tn = omp_get_thread_num()
643	!$OMP DO
644	DO i = nxl, nxr
645	DO j = nys, nyn
646	sums_l_etot = 0.0_wp
647	DO k = nzb, nzt+1
648	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
649	!
650	!-- Prognostic and diagnostic variables
651	sums_l(k,3,tn) = sums_l(k,3,tn) + w(k,j,i) * rmask(j,i,sr) &
652	* flag
653	sums_l(k,8,tn) = sums_l(k,8,tn) + e(k,j,i) * rmask(j,i,sr) &
654	* flag
655	sums_l(k,9,tn) = sums_l(k,9,tn) + km(k,j,i) * rmask(j,i,sr) &
656	* flag
657	sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr) &
658	* flag
659	sums_l(k,40,tn) = sums_l(k,40,tn) + ( p(k,j,i) &
660	/ momentumflux_output_conversion(k) ) &
661	* flag
662
663	sums_l(k,33,tn) = sums_l(k,33,tn) + &
664	( pt(k,j,i)-hom(k,1,4,sr) )*2 rmask(j,i,sr)&
665	* flag
666
667	IF ( humidity ) THEN
668	sums_l(k,70,tn) = sums_l(k,70,tn) + &
669	( q(k,j,i)-hom(k,1,41,sr) )*2 rmask(j,i,sr)&
670	* flag
671	ENDIF
672	IF ( passive_scalar ) THEN
673	sums_l(k,116,tn) = sums_l(k,116,tn) + &
674	( s(k,j,i)-hom(k,1,115,sr) )*2 rmask(j,i,sr)&
675	* flag
676	ENDIF
677	!
678	!-- Higher moments
679	!-- (Computation of the skewness of w further below)
680	sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)*3 rmask(j,i,sr) &
681	* flag
682
683	sums_l_etot = sums_l_etot + &
684	0.5_wp * ( u(k,j,i)2 + v(k,j,i)2 + &
685	w(k,j,i)*2 ) rmask(j,i,sr)&
686	* flag
687	ENDDO
688	!
689	!-- Total and perturbation energy for the total domain (being
690	!-- collected in the last column of sums_l). Summation of these
691	!-- quantities is seperated from the previous loop in order to
692	!-- allow vectorization of that loop.
693	sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
694	!
695	!-- 2D-arrays (being collected in the last column of sums_l)
696	IF ( surf_def_h(0)%ns >= 1 ) THEN
697	m = surf_def_h(0)%start_index(j,i)
698	sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + &
699	surf_def_h(0)%us(m) * rmask(j,i,sr)
700	sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
701	surf_def_h(0)%usws(m) * rmask(j,i,sr)
702	sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
703	surf_def_h(0)%vsws(m) * rmask(j,i,sr)
704	sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
705	surf_def_h(0)%ts(m) * rmask(j,i,sr)
706	IF ( humidity ) THEN
707	sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
708	surf_def_h(0)%qs(m) * rmask(j,i,sr)
709	ENDIF
710	IF ( passive_scalar ) THEN
711	sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) + &
712	surf_def_h(0)%ss(m) * rmask(j,i,sr)
713	ENDIF
714	ENDIF
715	IF ( surf_lsm_h%ns >= 1 ) THEN
716	m = surf_lsm_h%start_index(j,i)
717	sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + &
718	surf_lsm_h%us(m) * rmask(j,i,sr)
719	sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
720	surf_lsm_h%usws(m) * rmask(j,i,sr)
721	sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
722	surf_lsm_h%vsws(m) * rmask(j,i,sr)
723	sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
724	surf_lsm_h%ts(m) * rmask(j,i,sr)
725	IF ( humidity ) THEN
726	sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
727	surf_lsm_h%qs(m) * rmask(j,i,sr)
728	ENDIF
729	IF ( passive_scalar ) THEN
730	sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) + &
731	surf_lsm_h%ss(m) * rmask(j,i,sr)
732	ENDIF
733	ENDIF
734	IF ( surf_usm_h%ns >= 1 ) THEN
735	m = surf_lsm_h%start_index(j,i)
736	sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + &
737	surf_usm_h%us(m) * rmask(j,i,sr)
738	sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
739	surf_usm_h%usws(m) * rmask(j,i,sr)
740	sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
741	surf_usm_h%vsws(m) * rmask(j,i,sr)
742	sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
743	surf_usm_h%ts(m) * rmask(j,i,sr)
744	IF ( humidity ) THEN
745	sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
746	surf_usm_h%qs(m) * rmask(j,i,sr)
747	ENDIF
748	IF ( passive_scalar ) THEN
749	sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) + &
750	surf_usm_h%ss(m) * rmask(j,i,sr)
751	ENDIF
752	ENDIF
753	ENDDO
754	ENDDO
755
756	!
757	!-- Computation of statistics when ws-scheme is not used. Else these
758	!-- quantities are evaluated in the advection routines.
759	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp ) &
760	THEN
761	!$OMP DO
762	DO i = nxl, nxr
763	DO j = nys, nyn
764	DO k = nzb, nzt+1
765	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
766
767	u2 = u(k,j,i)**2
768	v2 = v(k,j,i)**2
769	w2 = w(k,j,i)**2
770	ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
771	vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
772
773	sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr) &
774	* flag
775	sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr) &
776	* flag
777	sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr) &
778	* flag
779	!
780	!-- Perturbation energy
781
782	sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp * &
783	( ust2 + vst2 + w2 ) * rmask(j,i,sr) &
784	* flag
785	ENDDO
786	ENDDO
787	ENDDO
788	ENDIF
789	!
790	!-- Computaion of domain-averaged perturbation energy. Please note,
791	!-- to prevent that perturbation energy is larger (even if only slightly)
792	!-- than the total kinetic energy, calculation is based on deviations from
793	!-- the horizontal mean, instead of spatial descretization of the advection
794	!-- term.
795	!$OMP DO
796	DO i = nxl, nxr
797	DO j = nys, nyn
798	DO k = nzb, nzt+1
799	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
800
801	w2 = w(k,j,i)**2
802	ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
803	vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
804	w2 = w(k,j,i)**2
805
806	sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) &
807	+ 0.5_wp * ( ust2 + vst2 + w2 ) &
808	* rmask(j,i,sr) &
809	* flag
810	ENDDO
811	ENDDO
812	ENDDO
813
814	!
815	!-- Horizontally averaged profiles of the vertical fluxes
816
817	!$OMP DO
818	DO i = nxl, nxr
819	DO j = nys, nyn
820	!
821	!-- Subgridscale fluxes (without Prandtl layer from k=nzb,
822	!-- oterwise from k=nzb+1)
823	!-- NOTE: for simplicity, nzb_diff_s_inner is used below, although
824	!-- ---- strictly speaking the following k-loop would have to be
825	!-- split up according to the staggered grid.
826	!-- However, this implies no error since staggered velocity
827	!-- components are zero at the walls and inside buildings.
828	!-- Flag 23 is used to mask surface fluxes as well as model-top fluxes,
829	!-- which are added further below.
830	DO k = nzb, nzt
831	flag = MERGE( 1.0_wp, 0.0_wp, &
832	BTEST( wall_flags_0(k,j,i), 23 ) ) * &
833	MERGE( 1.0_wp, 0.0_wp, &
834	BTEST( wall_flags_0(k,j,i), 9 ) )
835	!
836	!-- Momentum flux w"u"
837	sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * ( &
838	km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
839	) * ( &
840	( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
841	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
842	) * rmask(j,i,sr) &
843	* rho_air_zw(k) &
844	* momentumflux_output_conversion(k) &
845	* flag
846	!
847	!-- Momentum flux w"v"
848	sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * ( &
849	km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
850	) * ( &
851	( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
852	+ ( w(k,j,i) - w(k,j-1,i) ) * ddy &
853	) * rmask(j,i,sr) &
854	* rho_air_zw(k) &
855	* momentumflux_output_conversion(k) &
856	* flag
857	!
858	!-- Heat flux w"pt"
859	sums_l(k,16,tn) = sums_l(k,16,tn) &
860	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
861	* ( pt(k+1,j,i) - pt(k,j,i) ) &
862	* rho_air_zw(k) &
863	* heatflux_output_conversion(k) &
864	* ddzu(k+1) * rmask(j,i,sr) &
865	* flag
866
867
868	!
869	!-- Salinity flux w"sa"
870	IF ( ocean ) THEN
871	sums_l(k,65,tn) = sums_l(k,65,tn) &
872	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
873	* ( sa(k+1,j,i) - sa(k,j,i) ) &
874	* ddzu(k+1) * rmask(j,i,sr) &
875	* flag
876	ENDIF
877
878	!
879	!-- Buoyancy flux, water flux (humidity flux) w"q"
880	IF ( humidity ) THEN
881	sums_l(k,45,tn) = sums_l(k,45,tn) &
882	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
883	* ( vpt(k+1,j,i) - vpt(k,j,i) ) &
884	* rho_air_zw(k) &
885	* heatflux_output_conversion(k) &
886	* ddzu(k+1) * rmask(j,i,sr) * flag
887	sums_l(k,48,tn) = sums_l(k,48,tn) &
888	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
889	* ( q(k+1,j,i) - q(k,j,i) ) &
890	* rho_air_zw(k) &
891	* waterflux_output_conversion(k)&
892	* ddzu(k+1) * rmask(j,i,sr) * flag
893
894	IF ( cloud_physics ) THEN
895	sums_l(k,51,tn) = sums_l(k,51,tn) &
896	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
897	* ( ( q(k+1,j,i) - ql(k+1,j,i) )&
898	- ( q(k,j,i) - ql(k,j,i) ) ) &
899	* rho_air_zw(k) &
900	* waterflux_output_conversion(k)&
901	* ddzu(k+1) * rmask(j,i,sr) * flag
902	ENDIF
903	ENDIF
904
905	!
906	!-- Passive scalar flux
907	IF ( passive_scalar ) THEN
908	sums_l(k,117,tn) = sums_l(k,117,tn) &
909	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
910	* ( s(k+1,j,i) - s(k,j,i) ) &
911	* ddzu(k+1) * rmask(j,i,sr) &
912	* flag
913	ENDIF
914
915	ENDDO
916
917	!
918	!-- Subgridscale fluxes in the Prandtl layer
919	IF ( use_surface_fluxes ) THEN
920	DO l = 0, 1
921	ki = MERGE( -1, 0, l == 0 )
922	IF ( surf_def_h(l)%ns >= 1 ) THEN
923	DO m = surf_def_h(l)%start_index(j,i), surf_def_h(l)%end_index(j,i)
924	k = surf_def_h(l)%k(m)
925
926	sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + &
927	momentumflux_output_conversion(k+ki) * &
928	surf_def_h(l)%usws(m) * rmask(j,i,sr) ! w"u"
929	sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + &
930	momentumflux_output_conversion(k+ki) * &
931	surf_def_h(l)%vsws(m) * rmask(j,i,sr) ! w"v"
932	sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + &
933	heatflux_output_conversion(k+ki) * &
934	surf_def_h(l)%shf(m) * rmask(j,i,sr) ! w"pt"
935	sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + &
936	0.0_wp * rmask(j,i,sr) ! u"pt"
937	sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
938	0.0_wp * rmask(j,i,sr) ! v"pt"
939	IF ( ocean ) THEN
940	sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
941	surf_def_h(l)%sasws(m) * rmask(j,i,sr) ! w"sa"
942	ENDIF
943	IF ( humidity ) THEN
944	sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) + &
945	waterflux_output_conversion(k+ki) * &
946	surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
947	sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + ( &
948	( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) * &
949	surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) * &
950	surf_def_h(l)%qsws(m) ) &
951	* heatflux_output_conversion(k+ki)
952	IF ( cloud_droplets ) THEN
953	sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + ( &
954	( 1.0_wp + 0.61_wp * q(k+ki,j,i) - &
955	ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) + &
956	0.61_wp * pt(k+ki,j,i) * surf_def_h(l)%qsws(m) ) &
957	* heatflux_output_conversion(k+ki)
958	ENDIF
959	IF ( cloud_physics ) THEN
960	!
961	!-- Formula does not work if ql(k+ki) /= 0.0
962	sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) + &
963	waterflux_output_conversion(k+ki) * &
964	surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
965	ENDIF
966	ENDIF
967	IF ( passive_scalar ) THEN
968	sums_l(k+ki,117,tn) = sums_l(k+ki,117,tn) + &
969	surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
970	ENDIF
971
972	ENDDO
973
974	ENDIF
975	ENDDO
976	IF ( surf_lsm_h%ns >= 1 ) THEN
977	m = surf_lsm_h%start_index(j,i)
978	sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
979	momentumflux_output_conversion(nzb) * &
980	surf_lsm_h%usws(m) * rmask(j,i,sr) ! w"u"
981	sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
982	momentumflux_output_conversion(nzb) * &
983	surf_lsm_h%vsws(m) * rmask(j,i,sr) ! w"v"
984	sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
985	heatflux_output_conversion(nzb) * &
986	surf_lsm_h%shf(m) * rmask(j,i,sr) ! w"pt"
987	sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
988	0.0_wp * rmask(j,i,sr) ! u"pt"
989	sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
990	0.0_wp * rmask(j,i,sr) ! v"pt"
991	IF ( ocean ) THEN
992	sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
993	surf_lsm_h%sasws(m) * rmask(j,i,sr) ! w"sa"
994	ENDIF
995	IF ( humidity ) THEN
996	sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
997	waterflux_output_conversion(nzb) * &
998	surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
999	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
1000	( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * &
1001	surf_lsm_h%shf(m) + 0.61_wp * pt(nzb,j,i) * &
1002	surf_lsm_h%qsws(m) ) &
1003	* heatflux_output_conversion(nzb)
1004	IF ( cloud_droplets ) THEN
1005	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
1006	( 1.0_wp + 0.61_wp * q(nzb,j,i) - &
1007	ql(nzb,j,i) ) * surf_lsm_h%shf(m) + &
1008	0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) ) &
1009	* heatflux_output_conversion(nzb)
1010	ENDIF
1011	IF ( cloud_physics ) THEN
1012	!
1013	!-- Formula does not work if ql(nzb) /= 0.0
1014	sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &
1015	waterflux_output_conversion(nzb) * &
1016	surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1017	ENDIF
1018	ENDIF
1019	IF ( passive_scalar ) THEN
1020	sums_l(nzb,117,tn) = sums_l(nzb,117,tn) + &
1021	surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
1022	ENDIF
1023
1024
1025	ENDIF
1026	IF ( surf_usm_h%ns >= 1 ) THEN
1027	m = surf_usm_h%start_index(j,i)
1028	sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
1029	momentumflux_output_conversion(nzb) * &
1030	surf_usm_h%usws(m) * rmask(j,i,sr) ! w"u"
1031	sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
1032	momentumflux_output_conversion(nzb) * &
1033	surf_usm_h%vsws(m) * rmask(j,i,sr) ! w"v"
1034	sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
1035	heatflux_output_conversion(nzb) * &
1036	surf_usm_h%shf(m) * rmask(j,i,sr) ! w"pt"
1037	sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
1038	0.0_wp * rmask(j,i,sr) ! u"pt"
1039	sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
1040	0.0_wp * rmask(j,i,sr) ! v"pt"
1041	IF ( ocean ) THEN
1042	sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
1043	surf_usm_h%sasws(m) * rmask(j,i,sr) ! w"sa"
1044	ENDIF
1045	IF ( humidity ) THEN
1046	sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
1047	waterflux_output_conversion(nzb) * &
1048	surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1049	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
1050	( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * &
1051	surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) * &
1052	surf_usm_h%qsws(m) ) &
1053	* heatflux_output_conversion(nzb)
1054	IF ( cloud_droplets ) THEN
1055	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
1056	( 1.0_wp + 0.61_wp * q(nzb,j,i) - &
1057	ql(nzb,j,i) ) * surf_usm_h%shf(m) + &
1058	0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) ) &
1059	* heatflux_output_conversion(nzb)
1060	ENDIF
1061	IF ( cloud_physics ) THEN
1062	!
1063	!-- Formula does not work if ql(nzb) /= 0.0
1064	sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &
1065	waterflux_output_conversion(nzb) * &
1066	surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1067	ENDIF
1068	ENDIF
1069	IF ( passive_scalar ) THEN
1070	sums_l(nzb,117,tn) = sums_l(nzb,117,tn) + &
1071	surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
1072	ENDIF
1073
1074
1075	ENDIF
1076
1077	ENDIF
1078
1079	IF ( .NOT. neutral ) THEN
1080	IF ( surf_def_h(0)%ns >= 1 ) THEN
1081	m = surf_def_h(0)%start_index(j,i)
1082	sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + &
1083	surf_def_h(0)%ol(m) * rmask(j,i,sr) ! L
1084	ENDIF
1085	IF ( surf_lsm_h%ns >= 1 ) THEN
1086	m = surf_lsm_h%start_index(j,i)
1087	sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + &
1088	surf_lsm_h%ol(m) * rmask(j,i,sr) ! L
1089	ENDIF
1090	IF ( surf_usm_h%ns >= 1 ) THEN
1091	m = surf_usm_h%start_index(j,i)
1092	sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + &
1093	surf_usm_h%ol(m) * rmask(j,i,sr) ! L
1094	ENDIF
1095	ENDIF
1096
1097	IF ( radiation ) THEN
1098	sums_l(nzb,99,tn) = sums_l(nzb,99,tn) + rad_net(j,i)
1099	sums_l(nzb,100,tn) = sums_l(nzb,100,tn) + rad_lw_in(nzb,j,i)
1100	sums_l(nzb,101,tn) = sums_l(nzb,101,tn) + rad_lw_out(nzb,j,i)
1101	sums_l(nzb,102,tn) = sums_l(nzb,102,tn) + rad_sw_in(nzb,j,i)
1102	sums_l(nzb,103,tn) = sums_l(nzb,103,tn) + rad_sw_out(nzb,j,i)
1103
1104	#if defined ( __rrtmg )
1105	IF ( radiation_scheme == 'rrtmg' ) THEN
1106	sums_l(nzb,108,tn) = sums_l(nzb,108,tn) + rrtm_aldif(0,j,i)
1107	sums_l(nzb,109,tn) = sums_l(nzb,109,tn) + rrtm_aldir(0,j,i)
1108	sums_l(nzb,110,tn) = sums_l(nzb,110,tn) + rrtm_asdif(0,j,i)
1109	sums_l(nzb,111,tn) = sums_l(nzb,111,tn) + rrtm_asdir(0,j,i)
1110	ENDIF
1111	#endif
1112	ENDIF
1113	!
1114	!-- Subgridscale fluxes at the top surface
1115	IF ( use_top_fluxes ) THEN
1116	m = surf_def_h(2)%start_index(j,i)
1117	sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
1118	momentumflux_output_conversion(nzt:nzt+1) * &
1119	surf_def_h(2)%usws(m) * rmask(j,i,sr) ! w"u"
1120	sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
1121	momentumflux_output_conversion(nzt:nzt+1) * &
1122	surf_def_h(2)%vsws(m) * rmask(j,i,sr) ! w"v"
1123	sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
1124	heatflux_output_conversion(nzt:nzt+1) * &
1125	surf_def_h(2)%shf(m) * rmask(j,i,sr) ! w"pt"
1126	sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
1127	0.0_wp * rmask(j,i,sr) ! u"pt"
1128	sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
1129	0.0_wp * rmask(j,i,sr) ! v"pt"
1130
1131	IF ( ocean ) THEN
1132	sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
1133	surf_def_h(2)%sasws(m) * rmask(j,i,sr) ! w"sa"
1134	ENDIF
1135	IF ( humidity ) THEN
1136	sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
1137	waterflux_output_conversion(nzt) * &
1138	surf_def_h(2)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1139	sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( &
1140	( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * &
1141	surf_def_h(2)%shf(m) + &
1142	0.61_wp * pt(nzt,j,i) * &
1143	surf_def_h(2)%qsws(m) ) &
1144	* heatflux_output_conversion(nzt)
1145	IF ( cloud_droplets ) THEN
1146	sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( &
1147	( 1.0_wp + 0.61_wp * q(nzt,j,i) - &
1148	ql(nzt,j,i) ) * &
1149	surf_def_h(2)%shf(m) + &
1150	0.61_wp * pt(nzt,j,i) * &
1151	surf_def_h(2)%qsws(m) )&
1152	* heatflux_output_conversion(nzt)
1153	ENDIF
1154	IF ( cloud_physics ) THEN
1155	!
1156	!-- Formula does not work if ql(nzb) /= 0.0
1157	sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + & ! w"q" (w"qv")
1158	waterflux_output_conversion(nzt) * &
1159	surf_def_h(2)%qsws(m) * rmask(j,i,sr)
1160	ENDIF
1161	ENDIF
1162	IF ( passive_scalar ) THEN
1163	sums_l(nzt,117,tn) = sums_l(nzt,117,tn) + &
1164	surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
1165	ENDIF
1166	ENDIF
1167
1168	!
1169	!-- Resolved fluxes (can be computed for all horizontal points)
1170	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
1171	!-- ---- speaking the following k-loop would have to be split up and
1172	!-- rearranged according to the staggered grid.
1173	DO k = nzb, nzt
1174	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
1175	ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + &
1176	u(k+1,j,i) - hom(k+1,1,1,sr) )
1177	vst = 0.5_wp * ( v(k,j,i) - hom(k,1,2,sr) + &
1178	v(k+1,j,i) - hom(k+1,1,2,sr) )
1179	pts = 0.5_wp * ( pt(k,j,i) - hom(k,1,4,sr) + &
1180	pt(k+1,j,i) - hom(k+1,1,4,sr) )
1181
1182	!-- Higher moments
1183	sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)*2 &
1184	rmask(j,i,sr) * flag
1185	sums_l(k,36,tn) = sums_l(k,36,tn) + pts*2 w(k,j,i) * &
1186	rmask(j,i,sr) * flag
1187
1188	!
1189	!-- Salinity flux and density (density does not belong to here,
1190	!-- but so far there is no other suitable place to calculate)
1191	IF ( ocean ) THEN
1192	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
1193	pts = 0.5_wp * ( sa(k,j,i) - hom(k,1,23,sr) + &
1194	sa(k+1,j,i) - hom(k+1,1,23,sr) )
1195	sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
1196	rmask(j,i,sr) * flag
1197	ENDIF
1198	sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) * &
1199	rmask(j,i,sr) * flag
1200	sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
1201	rmask(j,i,sr) * flag
1202	ENDIF
1203
1204	!
1205	!-- Buoyancy flux, water flux, humidity flux, liquid water
1206	!-- content, rain drop concentration and rain water content
1207	IF ( humidity ) THEN
1208	IF ( cloud_physics .OR. cloud_droplets ) THEN
1209	pts = 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + &
1210	vpt(k+1,j,i) - hom(k+1,1,44,sr) )
1211	sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
1212	heatflux_output_conversion(k) * &
1213	rmask(j,i,sr) * flag
1214	sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) &
1215	* flag
1216
1217	IF ( .NOT. cloud_droplets ) THEN
1218	pts = 0.5_wp * &
1219	( ( q(k,j,i) - ql(k,j,i) ) - &
1220	hom(k,1,42,sr) + &
1221	( q(k+1,j,i) - ql(k+1,j,i) ) - &
1222	hom(k+1,1,42,sr) )
1223	sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
1224	waterflux_output_conversion(k) * &
1225	rmask(j,i,sr) * &
1226	flag
1227	sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) * &
1228	rmask(j,i,sr) * &
1229	flag
1230	sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) * &
1231	rmask(j,i,sr) * &
1232	flag
1233	IF ( microphysics_morrison ) THEN
1234	sums_l(k,123,tn) = sums_l(k,123,tn) + nc(k,j,i) * &
1235	rmask(j,i,sr) *&
1236	flag
1237	ENDIF
1238	IF ( microphysics_seifert ) THEN
1239	sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * &
1240	rmask(j,i,sr) *&
1241	flag
1242	sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * &
1243	rmask(j,i,sr) *&
1244	flag
1245	ENDIF
1246	ENDIF
1247
1248	ELSE
1249	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
1250	pts = 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + &
1251	vpt(k+1,j,i) - hom(k+1,1,44,sr) )
1252	sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
1253	heatflux_output_conversion(k) * &
1254	rmask(j,i,sr) * &
1255	flag
1256	ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN
1257	sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * &
1258	hom(k,1,41,sr) ) * &
1259	sums_l(k,17,tn) + &
1260	0.61_wp * hom(k,1,4,sr) * &
1261	sums_l(k,49,tn) &
1262	) * heatflux_output_conversion(k) * &
1263	flag
1264	END IF
1265	END IF
1266	ENDIF
1267	!
1268	!-- Passive scalar flux
1269	IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca &
1270	.OR. sr /= 0 ) ) THEN
1271	pts = 0.5_wp * ( s(k,j,i) - hom(k,1,115,sr) + &
1272	s(k+1,j,i) - hom(k+1,1,115,sr) )
1273	sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) * &
1274	rmask(j,i,sr) * flag
1275	ENDIF
1276
1277	!
1278	!-- Energy flux we
1279	!-- has to be adjusted
1280	sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp * &
1281	( ust2 + vst2 + w(k,j,i)**2 ) &
1282	* momentumflux_output_conversion(k) &
1283	* rmask(j,i,sr) * flag
1284	ENDDO
1285	ENDDO
1286	ENDDO
1287	!$OMP END PARALLEL
1288	!
1289	!-- Treat land-surface quantities according to new wall model structure.
1290	IF ( land_surface ) THEN
1291	tn = 0
1292	!$OMP PARALLEL PRIVATE( i, j, m, tn )
1293	!$ tn = omp_get_thread_num()
1294	!$OMP DO
1295	DO m = 1, surf_lsm_h%ns
1296	i = surf_lsm_h%i(m)
1297	j = surf_lsm_h%j(m)
1298
1299	IF ( i >= nxl .AND. i <= nxr .AND. &
1300	j >= nys .AND. j <= nyn ) THEN
1301	sums_l(nzb,93,tn) = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m)
1302	sums_l(nzb,94,tn) = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m)
1303	sums_l(nzb,95,tn) = sums_l(nzb,95,tn) + surf_lsm_h%qsws_soil(m)
1304	sums_l(nzb,96,tn) = sums_l(nzb,96,tn) + surf_lsm_h%qsws_veg(m)
1305	sums_l(nzb,97,tn) = sums_l(nzb,97,tn) + surf_lsm_h%r_a(m)
1306	sums_l(nzb,98,tn) = sums_l(nzb,98,tn)+ surf_lsm_h%r_s(m)
1307	ENDIF
1308	ENDDO
1309	!$OMP END PARALLEL
1310
1311	tn = 0
1312	!$OMP PARALLEL PRIVATE( i, j, k, m, tn )
1313	!$ tn = omp_get_thread_num()
1314	!$OMP DO
1315	DO m = 1, surf_lsm_h%ns
1316
1317	i = surf_lsm_h%i(m)
1318	j = surf_lsm_h%j(m)
1319
1320	IF ( i >= nxl .AND. i <= nxr .AND. &
1321	j >= nys .AND. j <= nyn ) THEN
1322
1323	DO k = nzb_soil, nzt_soil
1324	sums_l(k,89,tn) = sums_l(k,89,tn) + t_soil_h%var_2d(k,m) &
1325	* rmask(j,i,sr)
1326	sums_l(k,91,tn) = sums_l(k,91,tn) + m_soil_h%var_2d(k,m) &
1327	* rmask(j,i,sr)
1328	ENDDO
1329	ENDIF
1330	ENDDO
1331	!$OMP END PARALLEL
1332	ENDIF
1333	!
1334	!-- For speed optimization fluxes which have been computed in part directly
1335	!-- inside the WS advection routines are treated seperatly
1336	!-- Momentum fluxes first:
1337
1338	tn = 0
1339	!$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
1340	!$ tn = omp_get_thread_num()
1341	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
1342	!$OMP DO
1343	DO i = nxl, nxr
1344	DO j = nys, nyn
1345	DO k = nzb, nzt
1346	!
1347	!-- Flag 23 is used to mask surface fluxes as well as model-top
1348	!-- fluxes, which are added further below.
1349	flag = MERGE( 1.0_wp, 0.0_wp, &
1350	BTEST( wall_flags_0(k,j,i), 23 ) ) * &
1351	MERGE( 1.0_wp, 0.0_wp, &
1352	BTEST( wall_flags_0(k,j,i), 9 ) )
1353
1354	ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + &
1355	u(k+1,j,i) - hom(k+1,1,1,sr) )
1356	vst = 0.5_wp * ( v(k,j,i) - hom(k,1,2,sr) + &
1357	v(k+1,j,i) - hom(k+1,1,2,sr) )
1358	!
1359	!-- Momentum flux wu
1360	sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp * &
1361	( w(k,j,i-1) + w(k,j,i) ) &
1362	* momentumflux_output_conversion(k) &
1363	* ust * rmask(j,i,sr) &
1364	* flag
1365	!
1366	!-- Momentum flux wv
1367	sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp * &
1368	( w(k,j-1,i) + w(k,j,i) ) &
1369	* momentumflux_output_conversion(k) &
1370	* vst * rmask(j,i,sr) &
1371	* flag
1372	ENDDO
1373	ENDDO
1374	ENDDO
1375
1376	ENDIF
1377	IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
1378	!$OMP DO
1379	DO i = nxl, nxr
1380	DO j = nys, nyn
1381	DO k = nzb, nzt
1382	flag = MERGE( 1.0_wp, 0.0_wp, &
1383	BTEST( wall_flags_0(k,j,i), 23 ) ) * &
1384	MERGE( 1.0_wp, 0.0_wp, &
1385	BTEST( wall_flags_0(k,j,i), 9 ) )
1386	!
1387	!-- Vertical heat flux
1388	sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp * &
1389	( pt(k,j,i) - hom(k,1,4,sr) + &
1390	pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
1391	* heatflux_output_conversion(k) &
1392	* w(k,j,i) * rmask(j,i,sr) * flag
1393	IF ( humidity ) THEN
1394	pts = 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) + &
1395	q(k+1,j,i) - hom(k+1,1,41,sr) )
1396	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
1397	waterflux_output_conversion(k) * &
1398	rmask(j,i,sr) * flag
1399	ENDIF
1400	IF ( passive_scalar ) THEN
1401	pts = 0.5_wp * ( s(k,j,i) - hom(k,1,115,sr) + &
1402	s(k+1,j,i) - hom(k+1,1,115,sr) )
1403	sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) * &
1404	rmask(j,i,sr) * flag
1405	ENDIF
1406	ENDDO
1407	ENDDO
1408	ENDDO
1409
1410	ENDIF
1411
1412	!
1413	!-- Density at top follows Neumann condition
1414	IF ( ocean ) THEN
1415	sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
1416	sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
1417	ENDIF
1418
1419	!
1420	!-- Divergence of vertical flux of resolved scale energy and pressure
1421	!-- fluctuations as well as flux of pressure fluctuation itself (68).
1422	!-- First calculate the products, then the divergence.
1423	!-- Calculation is time consuming. Do it only, if profiles shall be plotted.
1424	IF ( hom(nzb+1,2,55,0) /= 0.0_wp .OR. hom(nzb+1,2,68,0) /= 0.0_wp ) &
1425	THEN
1426	sums_ll = 0.0_wp ! local array
1427
1428	!$OMP DO
1429	DO i = nxl, nxr
1430	DO j = nys, nyn
1431	DO k = nzb+1, nzt
1432	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1433
1434	sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * ( &
1435	( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) ) &
1436	- 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
1437	+ ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) ) &
1438	- 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
1439	+ w(k,j,i)*2 ) flag
1440
1441	sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i) &
1442	* ( ( p(k,j,i) + p(k+1,j,i) ) &
1443	/ momentumflux_output_conversion(k) ) &
1444	* flag
1445
1446	ENDDO
1447	ENDDO
1448	ENDDO
1449	sums_ll(0,1) = 0.0_wp ! because w is zero at the bottom
1450	sums_ll(nzt+1,1) = 0.0_wp
1451	sums_ll(0,2) = 0.0_wp
1452	sums_ll(nzt+1,2) = 0.0_wp
1453
1454	DO k = nzb+1, nzt
1455	sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
1456	sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
1457	sums_l(k,68,tn) = sums_ll(k,2)
1458	ENDDO
1459	sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
1460	sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
1461	sums_l(nzb,68,tn) = 0.0_wp ! because w* = 0 at nzb
1462
1463	ENDIF
1464
1465	!
1466	!-- Divergence of vertical flux of SGS TKE and the flux itself (69)
1467	IF ( hom(nzb+1,2,57,0) /= 0.0_wp .OR. hom(nzb+1,2,69,0) /= 0.0_wp ) &
1468	THEN
1469	!$OMP DO
1470	DO i = nxl, nxr
1471	DO j = nys, nyn
1472	DO k = nzb+1, nzt
1473
1474	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1475
1476	sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * ( &
1477	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1478	- (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) &
1479	) * ddzw(k) &
1480	* flag
1481
1482	sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * ( &
1483	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1484	) * flag
1485
1486	ENDDO
1487	ENDDO
1488	ENDDO
1489	sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
1490	sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
1491
1492	ENDIF
1493
1494	!
1495	!-- Horizontal heat fluxes (subgrid, resolved, total).
1496	!-- Do it only, if profiles shall be plotted.
1497	IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
1498
1499	!$OMP DO
1500	DO i = nxl, nxr
1501	DO j = nys, nyn
1502	DO k = nzb+1, nzt
1503	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1504	!
1505	!-- Subgrid horizontal heat fluxes u"pt", v"pt"
1506	sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp * &
1507	( kh(k,j,i) + kh(k,j,i-1) ) &
1508	* ( pt(k,j,i-1) - pt(k,j,i) ) &
1509	* rho_air_zw(k) &
1510	* heatflux_output_conversion(k) &
1511	* ddx * rmask(j,i,sr) * flag
1512	sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp * &
1513	( kh(k,j,i) + kh(k,j-1,i) ) &
1514	* ( pt(k,j-1,i) - pt(k,j,i) ) &
1515	* rho_air_zw(k) &
1516	* heatflux_output_conversion(k) &
1517	* ddy * rmask(j,i,sr) * flag
1518	!
1519	!-- Resolved horizontal heat fluxes upt, vpt
1520	sums_l(k,59,tn) = sums_l(k,59,tn) + &
1521	( u(k,j,i) - hom(k,1,1,sr) ) &
1522	* 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
1523	pt(k,j,i) - hom(k,1,4,sr) ) &
1524	* heatflux_output_conversion(k) &
1525	* flag
1526	pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1527	pt(k,j,i) - hom(k,1,4,sr) )
1528	sums_l(k,62,tn) = sums_l(k,62,tn) + &
1529	( v(k,j,i) - hom(k,1,2,sr) ) &
1530	* 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1531	pt(k,j,i) - hom(k,1,4,sr) ) &
1532	* heatflux_output_conversion(k) &
1533	* flag
1534	ENDDO
1535	ENDDO
1536	ENDDO
1537	!
1538	!-- Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
1539	sums_l(nzb,58,tn) = 0.0_wp
1540	sums_l(nzb,59,tn) = 0.0_wp
1541	sums_l(nzb,60,tn) = 0.0_wp
1542	sums_l(nzb,61,tn) = 0.0_wp
1543	sums_l(nzb,62,tn) = 0.0_wp
1544	sums_l(nzb,63,tn) = 0.0_wp
1545
1546	ENDIF
1547	!$OMP END PARALLEL
1548
1549	!
1550	!-- Collect current large scale advection and subsidence tendencies for
1551	!-- data output
1552	IF ( large_scale_forcing .AND. ( simulated_time > 0.0_wp ) ) THEN
1553	!
1554	!-- Interpolation in time of LSF_DATA
1555	nt = 1
1556	DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1557	nt = nt + 1
1558	ENDDO
1559	IF ( simulated_time - dt_3d /= time_vert(nt) ) THEN
1560	nt = nt - 1
1561	ENDIF
1562
1563	fac = ( simulated_time - dt_3d - time_vert(nt) ) &
1564	/ ( time_vert(nt+1)-time_vert(nt) )
1565
1566
1567	DO k = nzb, nzt
1568	sums_ls_l(k,0) = td_lsa_lpt(k,nt) &
1569	+ fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1570	sums_ls_l(k,1) = td_lsa_q(k,nt) &
1571	+ fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
1572	ENDDO
1573
1574	sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1575	sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1576
1577	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
1578
1579	DO k = nzb, nzt
1580	sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac * &
1581	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1582	sums_ls_l(k,3) = td_sub_q(k,nt) + fac * &
1583	( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
1584	ENDDO
1585
1586	sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1587	sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1588
1589	ENDIF
1590
1591	ENDIF
1592
1593	tn = 0
1594	!$OMP PARALLEL PRIVATE( i, j, k, tn )
1595	!$ tn = omp_get_thread_num()
1596	IF ( radiation .AND. radiation_scheme == 'rrtmg' ) THEN
1597	!$OMP DO
1598	DO i = nxl, nxr
1599	DO j = nys, nyn
1600	DO k = nzb+1, nzt+1
1601	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1602
1603	sums_l(k,100,tn) = sums_l(k,100,tn) + rad_lw_in(k,j,i) &
1604	* rmask(j,i,sr) * flag
1605	sums_l(k,101,tn) = sums_l(k,101,tn) + rad_lw_out(k,j,i) &
1606	* rmask(j,i,sr) * flag
1607	sums_l(k,102,tn) = sums_l(k,102,tn) + rad_sw_in(k,j,i) &
1608	* rmask(j,i,sr) * flag
1609	sums_l(k,103,tn) = sums_l(k,103,tn) + rad_sw_out(k,j,i) &
1610	* rmask(j,i,sr) * flag
1611	sums_l(k,104,tn) = sums_l(k,104,tn) + rad_lw_cs_hr(k,j,i) &
1612	* rmask(j,i,sr) * flag
1613	sums_l(k,105,tn) = sums_l(k,105,tn) + rad_lw_hr(k,j,i) &
1614	* rmask(j,i,sr) * flag
1615	sums_l(k,106,tn) = sums_l(k,106,tn) + rad_sw_cs_hr(k,j,i) &
1616	* rmask(j,i,sr) * flag
1617	sums_l(k,107,tn) = sums_l(k,107,tn) + rad_sw_hr(k,j,i) &
1618	* rmask(j,i,sr) * flag
1619	ENDDO
1620	ENDDO
1621	ENDDO
1622	ENDIF
1623	!
1624	!-- Calculate the user-defined profiles
1625	CALL user_statistics( 'profiles', sr, tn )
1626	!$OMP END PARALLEL
1627
1628	!
1629	!-- Summation of thread sums
1630	IF ( threads_per_task > 1 ) THEN
1631	DO i = 1, threads_per_task-1
1632	sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i)
1633	sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1634	sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1635	sums_l(:,45:pr_palm,i)
1636	IF ( max_pr_user > 0 ) THEN
1637	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1638	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1639	sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1640	ENDIF
1641	ENDDO
1642	ENDIF
1643
1644	#if defined( __parallel )
1645
1646	!
1647	!-- Compute total sum from local sums
1648	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1649	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
1650	MPI_SUM, comm2d, ierr )
1651	IF ( large_scale_forcing ) THEN
1652	CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls, &
1653	MPI_REAL, MPI_SUM, comm2d, ierr )
1654	ENDIF
1655	#else
1656	sums = sums_l(:,:,0)
1657	IF ( large_scale_forcing ) THEN
1658	sums(:,81:88) = sums_ls_l
1659	ENDIF
1660	#endif
1661
1662	!
1663	!-- Final values are obtained by division by the total number of grid points
1664	!-- used for summation. After that store profiles.
1665	!-- Check, if statistical regions do contain at least one grid point at the
1666	!-- respective k-level, otherwise division by zero will lead to undefined
1667	!-- values, which may cause e.g. problems with NetCDF output
1668	!-- Profiles:
1669	DO k = nzb, nzt+1
1670	sums(k,3) = sums(k,3) / ngp_2dh(sr)
1671	sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr)
1672	sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr)
1673	sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr)
1674	sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr)
1675	sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr)
1676	sums(k,81:88) = sums(k,81:88) / ngp_2dh(sr)
1677	sums(k,89:112) = sums(k,89:112) / ngp_2dh(sr)
1678	sums(k,114) = sums(k,114) / ngp_2dh(sr)
1679	sums(k,117) = sums(k,117) / ngp_2dh(sr)
1680	IF ( ngp_2dh_s_inner(k,sr) /= 0 ) THEN
1681	sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr)
1682	sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr)
1683	sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr)
1684	sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr)
1685	sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr)
1686	sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr)
1687	sums(k,70:80) = sums(k,70:80) / ngp_2dh_s_inner(k,sr)
1688	sums(k,116) = sums(k,116) / ngp_2dh_s_inner(k,sr)
1689	sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
1690	ENDIF
1691	ENDDO
1692
1693	!-- u* and so on
1694	!-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1695	!-- size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1696	!-- above the topography, they are being divided by ngp_2dh(sr)
1697	sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / &
1698	ngp_2dh(sr)
1699	sums(nzb+12,pr_palm) = sums(nzb+12,pr_palm) / & ! qs
1700	ngp_2dh(sr)
1701	sums(nzb+13,pr_palm) = sums(nzb+13,pr_palm) / & ! ss
1702	ngp_2dh(sr)
1703	!-- eges, e*
1704	sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / &
1705	ngp_3d(sr)
1706	!-- Old and new divergence
1707	sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1708	ngp_3d_inner(sr)
1709
1710	!-- User-defined profiles
1711	IF ( max_pr_user > 0 ) THEN
1712	DO k = nzb, nzt+1
1713	sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1714	sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1715	ngp_2dh_s_inner(k,sr)
1716	ENDDO
1717	ENDIF
1718
1719	!
1720	!-- Collect horizontal average in hom.
1721	!-- Compute deduced averages (e.g. total heat flux)
1722	hom(:,1,3,sr) = sums(:,3) ! w
1723	hom(:,1,8,sr) = sums(:,8) ! e profiles 5-7 are initial profiles
1724	hom(:,1,9,sr) = sums(:,9) ! km
1725	hom(:,1,10,sr) = sums(:,10) ! kh
1726	hom(:,1,11,sr) = sums(:,11) ! l
1727	hom(:,1,12,sr) = sums(:,12) ! w"u"
1728	hom(:,1,13,sr) = sums(:,13) ! wu
1729	hom(:,1,14,sr) = sums(:,14) ! w"v"
1730	hom(:,1,15,sr) = sums(:,15) ! wv
1731	hom(:,1,16,sr) = sums(:,16) ! w"pt"
1732	hom(:,1,17,sr) = sums(:,17) ! wpt
1733	hom(:,1,18,sr) = sums(:,16) + sums(:,17) ! wpt
1734	hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu
1735	hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv
1736	hom(:,1,21,sr) = sums(:,21) ! wptBC
1737	hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC
1738	! profile 24 is initial profile (sa)
1739	! profiles 25-29 left empty for initial
1740	! profiles
1741	hom(:,1,30,sr) = sums(:,30) ! u*2
1742	hom(:,1,31,sr) = sums(:,31) ! v*2
1743	hom(:,1,32,sr) = sums(:,32) ! w*2
1744	hom(:,1,33,sr) = sums(:,33) ! pt*2
1745	hom(:,1,34,sr) = sums(:,34) ! e*
1746	hom(:,1,35,sr) = sums(:,35) ! w2pt
1747	hom(:,1,36,sr) = sums(:,36) ! wpt2
1748	hom(:,1,37,sr) = sums(:,37) ! we
1749	hom(:,1,38,sr) = sums(:,38) ! w*3
1750	hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp ! Sw
1751	hom(:,1,40,sr) = sums(:,40) ! p
1752	hom(:,1,45,sr) = sums(:,45) ! w"vpt"
1753	hom(:,1,46,sr) = sums(:,46) ! wvpt
1754	hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt
1755	hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv")
1756	hom(:,1,49,sr) = sums(:,49) ! wq (wqv)
1757	hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv)
1758	hom(:,1,51,sr) = sums(:,51) ! w"qv"
1759	hom(:,1,52,sr) = sums(:,52) ! wqv
1760	hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv)
1761	hom(:,1,54,sr) = sums(:,54) ! ql
1762	hom(:,1,55,sr) = sums(:,55) ! wuu*/dz
1763	hom(:,1,56,sr) = sums(:,56) ! wp/dz
1764	hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho_ocean )/dz
1765	hom(:,1,58,sr) = sums(:,58) ! u"pt"
1766	hom(:,1,59,sr) = sums(:,59) ! upt
1767	hom(:,1,60,sr) = sums(:,58) + sums(:,59) ! upt_t
1768	hom(:,1,61,sr) = sums(:,61) ! v"pt"
1769	hom(:,1,62,sr) = sums(:,62) ! vpt
1770	hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t
1771	hom(:,1,64,sr) = sums(:,64) ! rho_ocean
1772	hom(:,1,65,sr) = sums(:,65) ! w"sa"
1773	hom(:,1,66,sr) = sums(:,66) ! wsa
1774	hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa
1775	hom(:,1,68,sr) = sums(:,68) ! wp
1776	hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho_ocean
1777	hom(:,1,70,sr) = sums(:,70) ! q*2
1778	hom(:,1,71,sr) = sums(:,71) ! prho
1779	hom(:,1,72,sr) = hyp * 1E-2_wp ! hyp in hPa
1780	hom(:,1,123,sr) = sums(:,123) ! nc
1781	hom(:,1,73,sr) = sums(:,73) ! nr
1782	hom(:,1,74,sr) = sums(:,74) ! qr
1783	hom(:,1,75,sr) = sums(:,75) ! qc
1784	hom(:,1,76,sr) = sums(:,76) ! prr (precipitation rate)
1785	! 77 is initial density profile
1786	hom(:,1,78,sr) = ug ! ug
1787	hom(:,1,79,sr) = vg ! vg
1788	hom(:,1,80,sr) = w_subs ! w_subs
1789
1790	IF ( large_scale_forcing ) THEN
1791	hom(:,1,81,sr) = sums_ls_l(:,0) ! td_lsa_lpt
1792	hom(:,1,82,sr) = sums_ls_l(:,1) ! td_lsa_q
1793	IF ( use_subsidence_tendencies ) THEN
1794	hom(:,1,83,sr) = sums_ls_l(:,2) ! td_sub_lpt
1795	hom(:,1,84,sr) = sums_ls_l(:,3) ! td_sub_q
1796	ELSE
1797	hom(:,1,83,sr) = sums(:,83) ! td_sub_lpt
1798	hom(:,1,84,sr) = sums(:,84) ! td_sub_q
1799	ENDIF
1800	hom(:,1,85,sr) = sums(:,85) ! td_nud_lpt
1801	hom(:,1,86,sr) = sums(:,86) ! td_nud_q
1802	hom(:,1,87,sr) = sums(:,87) ! td_nud_u
1803	hom(:,1,88,sr) = sums(:,88) ! td_nud_v
1804	ENDIF
1805
1806	IF ( land_surface ) THEN
1807	hom(:,1,89,sr) = sums(:,89) ! t_soil
1808	! 90 is initial t_soil profile
1809	hom(:,1,91,sr) = sums(:,91) ! m_soil
1810	! 92 is initial m_soil profile
1811	hom(:,1,93,sr) = sums(:,93) ! ghf
1812	hom(:,1,94,sr) = sums(:,94) ! qsws_liq
1813	hom(:,1,95,sr) = sums(:,95) ! qsws_soil
1814	hom(:,1,96,sr) = sums(:,96) ! qsws_veg
1815	hom(:,1,97,sr) = sums(:,97) ! r_a
1816	hom(:,1,98,sr) = sums(:,98) ! r_s
1817
1818	ENDIF
1819
1820	IF ( radiation ) THEN
1821	hom(:,1,99,sr) = sums(:,99) ! rad_net
1822	hom(:,1,100,sr) = sums(:,100) ! rad_lw_in
1823	hom(:,1,101,sr) = sums(:,101) ! rad_lw_out
1824	hom(:,1,102,sr) = sums(:,102) ! rad_sw_in
1825	hom(:,1,103,sr) = sums(:,103) ! rad_sw_out
1826
1827	IF ( radiation_scheme == 'rrtmg' ) THEN
1828	hom(:,1,104,sr) = sums(:,104) ! rad_lw_cs_hr
1829	hom(:,1,105,sr) = sums(:,105) ! rad_lw_hr
1830	hom(:,1,106,sr) = sums(:,106) ! rad_sw_cs_hr
1831	hom(:,1,107,sr) = sums(:,107) ! rad_sw_hr
1832
1833	hom(:,1,108,sr) = sums(:,108) ! rrtm_aldif
1834	hom(:,1,109,sr) = sums(:,109) ! rrtm_aldir
1835	hom(:,1,110,sr) = sums(:,110) ! rrtm_asdif
1836	hom(:,1,111,sr) = sums(:,111) ! rrtm_asdir
1837	ENDIF
1838	ENDIF
1839
1840	hom(:,1,112,sr) = sums(:,112) !: L
1841
1842	IF ( passive_scalar ) THEN
1843	hom(:,1,117,sr) = sums(:,117) ! w"s"
1844	hom(:,1,114,sr) = sums(:,114) ! ws
1845	hom(:,1,118,sr) = sums(:,117) + sums(:,114) ! ws
1846	hom(:,1,116,sr) = sums(:,116) ! s*2
1847	ENDIF
1848
1849	hom(:,1,119,sr) = rho_air ! rho_air in Kg/m^3
1850	hom(:,1,120,sr) = rho_air_zw ! rho_air_zw in Kg/m^3
1851
1852	hom(:,1,pr_palm,sr) = sums(:,pr_palm)
1853	! u, w'u', w'v', t (in last profile)
1854
1855	IF ( max_pr_user > 0 ) THEN ! user-defined profiles
1856	hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1857	sums(:,pr_palm+1:pr_palm+max_pr_user)
1858	ENDIF
1859
1860	!
1861	!-- Determine the boundary layer height using two different schemes.
1862	!-- First scheme: Starting from the Earth's (Ocean's) surface, look for the
1863	!-- first relative minimum (maximum) of the total heat flux.
1864	!-- The corresponding height is assumed as the boundary layer height, if it
1865	!-- is less than 1.5 times the height where the heat flux becomes negative
1866	!-- (positive) for the first time.
1867	z_i(1) = 0.0_wp
1868	first = .TRUE.
1869
1870	IF ( ocean ) THEN
1871	DO k = nzt, nzb+1, -1
1872	IF ( first .AND. hom(k,1,18,sr) < -1.0E-8_wp ) THEN
1873	first = .FALSE.
1874	height = zw(k)
1875	ENDIF
1876	IF ( hom(k,1,18,sr) < -1.0E-8_wp .AND. &
1877	hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN
1878	IF ( zw(k) < 1.5_wp * height ) THEN
1879	z_i(1) = zw(k)
1880	ELSE
1881	z_i(1) = height
1882	ENDIF
1883	EXIT
1884	ENDIF
1885	ENDDO
1886	ELSE
1887	DO k = nzb, nzt-1
1888	IF ( first .AND. hom(k,1,18,sr) < -1.0E-8_wp ) THEN
1889	first = .FALSE.
1890	height = zw(k)
1891	ENDIF
1892	IF ( hom(k,1,18,sr) < -1.0E-8_wp .AND. &
1893	hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN
1894	IF ( zw(k) < 1.5_wp * height ) THEN
1895	z_i(1) = zw(k)
1896	ELSE
1897	z_i(1) = height
1898	ENDIF
1899	EXIT
1900	ENDIF
1901	ENDDO
1902	ENDIF
1903
1904	!
1905	!-- Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1906	!-- by Uhlenbrock(2006). The boundary layer height is the height with the
1907	!-- maximal local temperature gradient: starting from the second (the last
1908	!-- but one) vertical gridpoint, the local gradient must be at least
1909	!-- 0.2K/100m and greater than the next four gradients.
1910	!-- WARNING: The threshold value of 0.2K/100m must be adjusted for the
1911	!-- ocean case!
1912	z_i(2) = 0.0_wp
1913	DO k = nzb+1, nzt+1
1914	dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1915	ENDDO
1916	dptdz_threshold = 0.2_wp / 100.0_wp
1917
1918	IF ( ocean ) THEN
1919	DO k = nzt+1, nzb+5, -1
1920	IF ( dptdz(k) > dptdz_threshold .AND. &
1921	dptdz(k) > dptdz(k-1) .AND. dptdz(k) > dptdz(k-2) .AND. &
1922	dptdz(k) > dptdz(k-3) .AND. dptdz(k) > dptdz(k-4) ) THEN
1923	z_i(2) = zw(k-1)
1924	EXIT
1925	ENDIF
1926	ENDDO
1927	ELSE
1928	DO k = nzb+1, nzt-3
1929	IF ( dptdz(k) > dptdz_threshold .AND. &
1930	dptdz(k) > dptdz(k+1) .AND. dptdz(k) > dptdz(k+2) .AND. &
1931	dptdz(k) > dptdz(k+3) .AND. dptdz(k) > dptdz(k+4) ) THEN
1932	z_i(2) = zw(k-1)
1933	EXIT
1934	ENDIF
1935	ENDDO
1936	ENDIF
1937
1938	hom(nzb+6,1,pr_palm,sr) = z_i(1)
1939	hom(nzb+7,1,pr_palm,sr) = z_i(2)
1940
1941	!
1942	!-- Determine vertical index which is nearest to the mean surface level
1943	!-- height of the respective statistic region
1944	DO k = nzb, nzt
1945	IF ( zw(k) >= mean_surface_level_height(sr) ) THEN
1946	k_surface_level = k
1947	EXIT
1948	ENDIF
1949	ENDDO
1950	!
1951	!-- Computation of both the characteristic vertical velocity and
1952	!-- the characteristic convective boundary layer temperature.
1953	!-- The inversion height entering into the equation is defined with respect
1954	!-- to the mean surface level height of the respective statistic region.
1955	!-- The horizontal average at surface level index + 1 is input for the
1956	!-- average temperature.
1957	IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp .AND. z_i(1) /= 0.0_wp )&
1958	THEN
1959	hom(nzb+8,1,pr_palm,sr) = &
1960	( g / hom(k_surface_level+1,1,4,sr) * &
1961	( hom(k_surface_level,1,18,sr) / &
1962	( heatflux_output_conversion(nzb) * rho_air(nzb) ) ) &
1963	* ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
1964	ELSE
1965	hom(nzb+8,1,pr_palm,sr) = 0.0_wp
1966	ENDIF
1967
1968	!
1969	!-- Collect the time series quantities
1970	ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E
1971	ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*
1972	ts_value(3,sr) = dt_3d
1973	ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u*
1974	ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th*
1975	ts_value(6,sr) = u_max
1976	ts_value(7,sr) = v_max
1977	ts_value(8,sr) = w_max
1978	ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence
1979	ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence
1980	ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1)
1981	ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2)
1982	ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w*
1983	ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0
1984	ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1
1985	ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=1
1986	ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0)
1987	ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp)
1988	ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr) ! u'w' at k=0
1989	ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=0
1990	ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0
1991
1992	IF ( .NOT. neutral ) THEN
1993	ts_value(22,sr) = hom(nzb,1,112,sr) ! L
1994	ELSE
1995	ts_value(22,sr) = 1.0E10_wp
1996	ENDIF
1997
1998	ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q*
1999
2000	IF ( passive_scalar ) THEN
2001	ts_value(24,sr) = hom(nzb+13,1,117,sr) ! w"s" ( to do ! )
2002	ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr) ! s*
2003	ENDIF
2004
2005	!
2006	!-- Collect land surface model timeseries
2007	IF ( land_surface ) THEN
2008	ts_value(dots_soil ,sr) = hom(nzb,1,93,sr) ! ghf
2009	ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr) ! qsws_liq
2010	ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr) ! qsws_soil
2011	ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr) ! qsws_veg
2012	ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr) ! r_a
2013	ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr) ! r_s
2014	ENDIF
2015	!
2016	!-- Collect radiation model timeseries
2017	IF ( radiation ) THEN
2018	ts_value(dots_rad,sr) = hom(nzb,1,99,sr) ! rad_net
2019	ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr) ! rad_lw_in
2020	ts_value(dots_rad+2,sr) = hom(nzb,1,101,sr) ! rad_lw_out
2021	ts_value(dots_rad+3,sr) = hom(nzb,1,102,sr) ! rad_sw_in
2022	ts_value(dots_rad+4,sr) = hom(nzb,1,103,sr) ! rad_sw_out
2023
2024	IF ( radiation_scheme == 'rrtmg' ) THEN
2025	ts_value(dots_rad+5,sr) = hom(nzb,1,108,sr) ! rrtm_aldif
2026	ts_value(dots_rad+6,sr) = hom(nzb,1,109,sr) ! rrtm_aldir
2027	ts_value(dots_rad+7,sr) = hom(nzb,1,110,sr) ! rrtm_asdif
2028	ts_value(dots_rad+8,sr) = hom(nzb,1,111,sr) ! rrtm_asdir
2029	ENDIF
2030
2031	ENDIF
2032
2033	!
2034	!-- Calculate additional statistics provided by the user interface
2035	CALL user_statistics( 'time_series', sr, 0 )
2036
2037	ENDDO ! loop of the subregions
2038
2039	!
2040	!-- If required, sum up horizontal averages for subsequent time averaging.
2041	!-- Do not sum, if flow statistics is called before the first initial time step.
2042	IF ( do_sum .AND. simulated_time /= 0.0_wp ) THEN
2043	IF ( average_count_pr == 0 ) hom_sum = 0.0_wp
2044	hom_sum = hom_sum + hom(:,1,:,:)
2045	average_count_pr = average_count_pr + 1
2046	do_sum = .FALSE.
2047	ENDIF
2048
2049	!
2050	!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
2051	!-- This flag is reset after each time step in time_integration.
2052	flow_statistics_called = .TRUE.
2053
2054	CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
2055
2056
2057	END SUBROUTINE flow_statistics

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |