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

Last change on this file since 4581 was 4581, checked in by suehring, 4 years ago

mesoscale nesting: omit explicit pressure forcing via geostrophic wind components; chemistry: enable profile output of vertical fluxes; urban-surface: bugfix in initialization in case of cyclic_fill

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