Home

Context Navigation

flow_statistics.f90 @ 4177

Last change on this file since 4177 was 4131, checked in by monakurppa, 5 years ago

Several changes in the salsa aerosol module:

Add "salsa_" before each salsa output variable
Add a possibility to output the number (salsa_N_UFP) and mass concentration (salsa_PM0.1) of ultrafine particles, i.e. particles with a diameter smaller than 100 nm
Implement aerosol emission mode "parameterized" which is based on the street type (similar to the chemistry module).
Remove unnecessary nucleation subroutines.
Add the z-dimension for gaseous emissions to correspond the implementation in the chemistry module

Property svn:keywords set to Id

File size: 113.8 KB

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

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

Context Navigation

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

Download in other formats: