source: palm/trunk/SOURCE/turbulence_closure_mod.f90 @ 4433

Last change on this file since 4433 was 4433, checked in by gronemeier, 4 years ago

remove warning for newly implemented RANS mode (turbulence_closure_mod)

  • Property svn:keywords set to Id
File size: 202.8 KB
Line 
1!> @file turbulence_closure_mod.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 2017-2018 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: turbulence_closure_mod.f90 4433 2020-02-28 22:14:43Z gronemeier $
27! remove warning for newly implemented RANS mode
28!
29! 4370 2020-01-10 14:00:44Z raasch
30! vector directives added to force vectorization on Intel19 compiler
31!
32! 4346 2019-12-18 11:55:56Z motisi
33! Introduction of wall_flags_total_0, which currently sets bits based on static
34! topography information used in wall_flags_static_0
35!
36! 4329 2019-12-10 15:46:36Z motisi
37! Renamed wall_flags_0 to wall_flags_static_0
38!
39! 4182 2019-08-22 15:20:23Z scharf
40! Corrected "Former revisions" section
41!
42! 4177 2019-08-20 14:32:34Z gronemeier
43! add comment
44!
45! 4170 2019-08-19 17:12:31Z gronemeier
46! - add performance optimizations according to K. Ketelsen
47!   to diffusion_e and tcm_diffusivities_default
48! - bugfix in calculating l_wall for vertical walls
49! - bugfix in using l_wall in initialization (consider wall_adjustment_factor)
50! - always initialize diss and save the dissipation to that array
51!
52! 4168 2019-08-16 13:50:17Z suehring
53! Replace function get_topography_top_index by topo_top_ind
54!
55! 4110 2019-07-22 17:05:21Z suehring
56! pass integer flag array as well as boundary flags to WS scalar advection
57! routine
58!
59! 4109 2019-07-22 17:00:34Z suehring
60! - Modularize setting of boundary conditions for TKE and dissipation
61! - Neumann boundary condition for TKE at model top is set also in child domain
62! - Revise setting of Neumann boundary conditions at non-cyclic lateral
63!   boundaries
64! - Bugfix, set Neumann boundary condition for TKE at vertical wall instead of
65!   an implicit Dirichlet boundary condition which implied a sink of TKE
66!   at vertical walls
67!
68! 4048 2019-06-21 21:00:21Z knoop
69! write out preprocessor directives; remove tailing whitespaces
70!
71! 3775 2019-03-04 12:40:20Z gronemeier
72! removed unused variables
73!
74! 3724 2019-02-06 16:28:23Z kanani
75! Correct double-used log_point_s units
76!
77! 3719 2019-02-06 13:10:18Z kanani
78! Changed log_point to log_point_s, otherwise this overlaps with
79! 'all progn.equations' cpu measurement.
80!
81! 3684 2019-01-20 20:20:58Z knoop
82! Remove unused variable simulated_time
83!
84! 2696 2017-12-14 17:12:51Z kanani
85! Initial revision
86!
87!
88! Authors:
89! --------
90! @author Tobias Gronemeier
91! @author Hauke Wurps
92!
93! Description:
94! ------------
95!> This module contains the available turbulence closures for PALM.
96!>
97!>
98!> @todo test initialization for all possibilities
99!> @todo add OpenMP directives whereever possible
100!> @todo Check for random disturbances
101!> @note <Enter notes on the module>
102!-----------------------------------------------------------------------------!
103 MODULE turbulence_closure_mod
104
105
106    USE arrays_3d,                                                            &
107        ONLY:  diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3,   &
108               e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m,          &
109               te_m, tend, u, v, vpt, w
110
111    USE basic_constants_and_equations_mod,                                    &
112        ONLY:  g, kappa, lv_d_cp, lv_d_rd, rd_d_rv
113
114    USE control_parameters,                                                   &
115        ONLY:  bc_dirichlet_l,                                                &
116               bc_dirichlet_n,                                                &
117               bc_dirichlet_r,                                                &
118               bc_dirichlet_s,                                                &
119               bc_radiation_l,                                                &
120               bc_radiation_n,                                                &
121               bc_radiation_r,                                                &
122               bc_radiation_s,                                                &
123               child_domain,                                                  &
124               constant_diffusion, dt_3d, e_init, humidity,                   &
125               initializing_actions, intermediate_timestep_count,             &
126               intermediate_timestep_count_max, km_constant,                  &
127               les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, &
128               pt_reference, rans_mode, rans_tke_e, rans_tke_l,               &
129               timestep_scheme, turbulence_closure,                           &
130               turbulent_inflow, use_upstream_for_tke, vpt_reference,         &
131               ws_scheme_sca, current_timestep_number
132
133    USE advec_ws,                                                             &
134        ONLY:  advec_s_ws
135
136    USE advec_s_bc_mod,                                                       &
137        ONLY:  advec_s_bc
138
139    USE advec_s_pw_mod,                                                       &
140        ONLY:  advec_s_pw
141
142    USE advec_s_up_mod,                                                       &
143        ONLY:  advec_s_up
144
145    USE cpulog,                                                               &
146        ONLY:  cpu_log, log_point_s
147
148    USE indices,                                                              &
149        ONLY:  advc_flags_s,                                                  &
150               nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,    &
151               topo_top_ind,                                                  &
152               wall_flags_total_0
153
154    USE kinds
155
156    USE ocean_mod,                                                            &
157        ONLY:  prho_reference
158
159    USE pegrid
160
161    USE plant_canopy_model_mod,                                               &
162        ONLY:  pcm_tendency
163
164    USE statistics,                                                           &
165        ONLY:  hom, hom_sum, statistic_regions
166
167    USE surface_mod,                                                          &
168        ONLY:  bc_h,                                                          &
169               bc_v,                                                          &
170               surf_def_h,                                                    &
171               surf_def_v,                                                    &
172               surf_lsm_h,                                                    &
173               surf_lsm_v,                                                    &
174               surf_usm_h,                                                    &
175               surf_usm_v
176
177    IMPLICIT NONE
178
179
180    REAL(wp) ::  c_0                !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
181    REAL(wp) ::  c_1                !< model constant for RANS mode
182    REAL(wp) ::  c_2                !< model constant for RANS mode
183    REAL(wp) ::  c_3                !< model constant for RANS mode
184    REAL(wp) ::  c_4                !< model constant for RANS mode
185    REAL(wp) ::  l_max              !< maximum length scale for Blackadar mixing length
186    REAL(wp) ::  dsig_e = 1.0_wp    !< factor to calculate Ke from Km (1/sigma_e)
187    REAL(wp) ::  dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss)
188
189    REAL(wp), DIMENSION(0:4) :: rans_const_c = &       !< model constants for RANS mode (namelist param)
190       (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
191
192    REAL(wp), DIMENSION(2) :: rans_const_sigma = &     !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
193       (/ 1.0_wp, 1.30_wp /)
194
195    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black    !< mixing length according to Blackadar
196
197    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_grid     !< geometric mean of grid sizes dx, dy, dz
198
199    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
200
201!
202!-- Public variables
203    PUBLIC c_0, rans_const_c, rans_const_sigma
204
205    SAVE
206
207    PRIVATE
208!
209!-- Public subroutines
210    PUBLIC                                                                     &
211       tcm_boundary_conds,                                                     &
212       tcm_check_parameters,                                                   &
213       tcm_check_data_output,                                                  &
214       tcm_define_netcdf_grid,                                                 &
215       tcm_init_arrays,                                                        &
216       tcm_init,                                                               &
217       tcm_actions,                                                            &
218       tcm_prognostic_equations,                                               &
219       tcm_swap_timelevel,                                                     &
220       tcm_3d_data_averaging,                                                  &
221       tcm_data_output_2d,                                                     &
222       tcm_data_output_3d,                                                     &
223       tcm_diffusivities
224
225!
226!-- PALM interfaces:
227!-- Boundary conditions for subgrid TKE and dissipation
228    INTERFACE tcm_boundary_conds
229       MODULE PROCEDURE tcm_boundary_conds
230    END INTERFACE tcm_boundary_conds
231!
232!-- Input parameter checks to be done in check_parameters
233    INTERFACE tcm_check_parameters
234       MODULE PROCEDURE tcm_check_parameters
235    END INTERFACE tcm_check_parameters
236
237!
238!-- Data output checks for 2D/3D data to be done in check_parameters
239    INTERFACE tcm_check_data_output
240       MODULE PROCEDURE tcm_check_data_output
241    END INTERFACE tcm_check_data_output
242
243!
244!-- Definition of data output quantities
245    INTERFACE tcm_define_netcdf_grid
246       MODULE PROCEDURE tcm_define_netcdf_grid
247    END INTERFACE tcm_define_netcdf_grid
248
249!
250!-- Initialization of arrays
251    INTERFACE tcm_init_arrays
252       MODULE PROCEDURE tcm_init_arrays
253    END INTERFACE tcm_init_arrays
254
255!
256!-- Initialization actions
257    INTERFACE tcm_init
258       MODULE PROCEDURE tcm_init
259    END INTERFACE tcm_init
260
261!
262!-- Location specific actions
263    INTERFACE tcm_actions
264       MODULE PROCEDURE tcm_actions
265       MODULE PROCEDURE tcm_actions_ij
266    END INTERFACE tcm_actions
267
268!
269!-- Prognostic equations for TKE and TKE dissipation rate
270    INTERFACE tcm_prognostic_equations
271       MODULE PROCEDURE tcm_prognostic_equations
272       MODULE PROCEDURE tcm_prognostic_equations_ij
273    END INTERFACE tcm_prognostic_equations
274
275!
276!-- Swapping of time levels (required for prognostic variables)
277    INTERFACE tcm_swap_timelevel
278       MODULE PROCEDURE tcm_swap_timelevel
279    END INTERFACE tcm_swap_timelevel
280
281!
282!-- Averaging of 3D data for output
283    INTERFACE tcm_3d_data_averaging
284       MODULE PROCEDURE tcm_3d_data_averaging
285    END INTERFACE tcm_3d_data_averaging
286
287!
288!-- Data output of 2D quantities
289    INTERFACE tcm_data_output_2d
290       MODULE PROCEDURE tcm_data_output_2d
291    END INTERFACE tcm_data_output_2d
292
293!
294!-- Data output of 3D data
295    INTERFACE tcm_data_output_3d
296       MODULE PROCEDURE tcm_data_output_3d
297    END INTERFACE tcm_data_output_3d
298
299!
300!-- Call tcm_diffusivities_default and tcm_diffusivities_dynamic
301    INTERFACE tcm_diffusivities
302       MODULE PROCEDURE tcm_diffusivities
303    END INTERFACE tcm_diffusivities
304
305
306 CONTAINS
307
308!------------------------------------------------------------------------------!
309! Description:
310! ------------
311!> Check parameters routine for turbulence closure module.
312!------------------------------------------------------------------------------!
313 SUBROUTINE tcm_boundary_conds
314
315    USE pmc_interface,                                                         &
316        ONLY : rans_mode_parent
317
318    IMPLICIT NONE
319
320    INTEGER(iwp) ::  i  !< grid index x direction
321    INTEGER(iwp) ::  j  !< grid index y direction
322    INTEGER(iwp) ::  k  !< grid index z direction
323    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
324    INTEGER(iwp) ::  m  !< running index surface elements
325!
326!-- Boundary conditions for TKE.
327    IF ( .NOT. constant_diffusion )  THEN
328!
329!--    In LES mode, Neumann conditions with de/x_i=0 are assumed at solid walls.
330!--    Note, only TKE is prognostic in this case and dissipation is only
331!--    a diagnostic quantity.
332       IF ( .NOT. rans_mode )  THEN
333!
334!--       Horizontal walls, upward- and downward-facing
335          DO  l = 0, 1
336             !$OMP PARALLEL DO PRIVATE( i, j, k )
337             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
338             !$ACC PRESENT(bc_h, e_p)
339             DO  m = 1, bc_h(l)%ns
340                i = bc_h(l)%i(m)
341                j = bc_h(l)%j(m)
342                k = bc_h(l)%k(m)
343                e_p(k+bc_h(l)%koff,j,i) = e_p(k,j,i)
344             ENDDO
345          ENDDO
346!
347!--       Vertical walls
348          DO  l = 0, 3
349!
350!--          Note concerning missing ACC directive for this loop: Even though
351!--          the data structure bc_v is present, it may not contain any
352!--          allocated arrays in the flat but also in a topography case,
353!--          leading to a runtime error. Therefore, omit ACC directives
354!--          for this loop, in contrast to the bc_h loop.
355             !$OMP PARALLEL DO PRIVATE( i, j, k )
356             DO  m = 1, bc_v(l)%ns
357                i = bc_v(l)%i(m)
358                j = bc_v(l)%j(m)
359                k = bc_v(l)%k(m)
360                e_p(k,j+bc_v(l)%joff,i+bc_v(l)%ioff) = e_p(k,j,i)
361             ENDDO
362          ENDDO
363!
364!--    In RANS mode, wall function is used as boundary condition for TKE
365       ELSE
366!
367!--       Use wall function within constant-flux layer
368!--       Note, grid points listed in bc_h are not included in any calculations in RANS mode and
369!--       are therefore not set here.
370!
371!--       Upward-facing surfaces
372!--       Default surfaces
373          DO  m = 1, surf_def_h(0)%ns
374             i = surf_def_h(0)%i(m)
375             j = surf_def_h(0)%j(m)
376             k = surf_def_h(0)%k(m)
377             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
378          ENDDO
379!
380!--       Natural surfaces
381          DO  m = 1, surf_lsm_h%ns
382             i = surf_lsm_h%i(m)
383             j = surf_lsm_h%j(m)
384             k = surf_lsm_h%k(m)
385             e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2
386          ENDDO
387!
388!--       Urban surfaces
389          DO  m = 1, surf_usm_h%ns
390             i = surf_usm_h%i(m)
391             j = surf_usm_h%j(m)
392             k = surf_usm_h%k(m)
393             e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2
394          ENDDO
395!
396!--       Vertical surfaces
397          DO  l = 0, 3
398!
399!--          Default surfaces
400             DO  m = 1, surf_def_v(l)%ns
401                i = surf_def_v(l)%i(m)
402                j = surf_def_v(l)%j(m)
403                k = surf_def_v(l)%k(m)
404                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
405             ENDDO
406!
407!--          Natural surfaces
408             DO  m = 1, surf_lsm_v(l)%ns
409                i = surf_lsm_v(l)%i(m)
410                j = surf_lsm_v(l)%j(m)
411                k = surf_lsm_v(l)%k(m)
412                e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2
413             ENDDO
414!
415!--          Urban surfaces
416             DO  m = 1, surf_usm_v(l)%ns
417                i = surf_usm_v(l)%i(m)
418                j = surf_usm_v(l)%j(m)
419                k = surf_usm_v(l)%k(m)
420                e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2
421             ENDDO
422          ENDDO
423       ENDIF
424!
425!--    Set Neumann boundary condition for TKE at model top. Do this also
426!--    in case of a nested run.
427       !$ACC KERNELS PRESENT(e_p)
428       e_p(nzt+1,:,:) = e_p(nzt,:,:)
429       !$ACC END KERNELS
430!
431!--    Nesting case: if parent operates in RANS mode and child in LES mode,
432!--    no TKE is transfered. This case, set Neumann conditions at lateral and
433!--    top child boundaries.
434!--    If not ( both either in RANS or in LES mode ), TKE boundary condition
435!--    is treated in the nesting.
436       If ( child_domain )  THEN
437          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
438
439             e_p(nzt+1,:,:) = e_p(nzt,:,:)
440             IF ( bc_dirichlet_l )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
441             IF ( bc_dirichlet_r )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
442             IF ( bc_dirichlet_s )  e_p(:,nys-1,:) = e_p(:,nys,:)
443             IF ( bc_dirichlet_n )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
444
445          ENDIF
446       ENDIF
447!
448!--    At in- and outflow boundaries also set Neumann boundary conditions
449!--    for the SGS-TKE. An exception is made for the child domain if
450!--    both parent and child operate in RANS mode. This case no
451!--    lateral Neumann boundary conditions will be set but Dirichlet
452!--    conditions will be set in the nesting.
453       IF ( .NOT. child_domain  .AND.  .NOT. rans_mode_parent  .AND.           &
454            .NOT. rans_mode )  THEN
455          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
456             e_p(:,nys-1,:) = e_p(:,nys,:)
457             IF ( rans_tke_e )  diss_p(:,nys-1,:) = diss_p(:,nys,:)
458          ENDIF
459          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
460             e_p(:,nyn+1,:) = e_p(:,nyn,:)
461             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)
462          ENDIF
463          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
464             e_p(:,:,nxl-1) = e_p(:,:,nxl)
465             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)
466          ENDIF
467          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
468             e_p(:,:,nxr+1) = e_p(:,:,nxr)
469             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)
470          ENDIF
471       ENDIF
472    ENDIF
473
474!
475!-- Boundary conditions for TKE dissipation rate in RANS mode.
476    IF ( rans_tke_e )  THEN
477!
478!--    Use wall function within constant-flux layer
479!--    Upward-facing surfaces
480!--    Default surfaces
481       DO  m = 1, surf_def_h(0)%ns
482          i = surf_def_h(0)%i(m)
483          j = surf_def_h(0)%j(m)
484          k = surf_def_h(0)%k(m)
485          diss_p(k,j,i) = surf_def_h(0)%us(m)**3          &
486                        / ( kappa * surf_def_h(0)%z_mo(m) )
487       ENDDO
488!
489!--    Natural surfaces
490       DO  m = 1, surf_lsm_h%ns
491          i = surf_lsm_h%i(m)
492          j = surf_lsm_h%j(m)
493          k = surf_lsm_h%k(m)
494          diss_p(k,j,i) = surf_lsm_h%us(m)**3          &
495                        / ( kappa * surf_lsm_h%z_mo(m) )
496       ENDDO
497!
498!--    Urban surfaces
499       DO  m = 1, surf_usm_h%ns
500          i = surf_usm_h%i(m)
501          j = surf_usm_h%j(m)
502          k = surf_usm_h%k(m)
503          diss_p(k,j,i) = surf_usm_h%us(m)**3          &
504                        / ( kappa * surf_usm_h%z_mo(m) )
505       ENDDO
506!
507!--    Vertical surfaces
508       DO  l = 0, 3
509!
510!--       Default surfaces
511          DO  m = 1, surf_def_v(l)%ns
512             i = surf_def_v(l)%i(m)
513             j = surf_def_v(l)%j(m)
514             k = surf_def_v(l)%k(m)
515             diss_p(k,j,i) = surf_def_v(l)%us(m)**3          &
516                           / ( kappa * surf_def_v(l)%z_mo(m) )
517          ENDDO
518!
519!--       Natural surfaces
520          DO  m = 1, surf_lsm_v(l)%ns
521             i = surf_lsm_v(l)%i(m)
522             j = surf_lsm_v(l)%j(m)
523             k = surf_lsm_v(l)%k(m)
524             diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3          &
525                           / ( kappa * surf_lsm_v(l)%z_mo(m) )
526          ENDDO
527!
528!--       Urban surfaces
529          DO  m = 1, surf_usm_v(l)%ns
530             i = surf_usm_v(l)%i(m)
531             j = surf_usm_v(l)%j(m)
532             k = surf_usm_v(l)%k(m)
533             diss_p(k,j,i) = surf_usm_v(l)%us(m)**3          &
534                           / ( kappa * surf_usm_v(l)%z_mo(m) )
535          ENDDO
536       ENDDO
537!
538!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
539!--    minimum value
540       DO  i = nxl, nxr
541          DO  j = nys, nyn
542             DO  k = nzb, nzt+1
543                diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i),          &
544                                          2.0_wp * diss(k,j,i) ), &
545                                     0.1_wp * diss(k,j,i),        &
546                                     0.0001_wp )
547             ENDDO
548          ENDDO
549       ENDDO
550
551       diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
552
553    ENDIF
554
555 END SUBROUTINE tcm_boundary_conds
556
557!------------------------------------------------------------------------------!
558! Description:
559! ------------
560!> Check parameters routine for turbulence closure module.
561!------------------------------------------------------------------------------!
562 SUBROUTINE tcm_check_parameters
563
564    USE control_parameters,                                                    &
565        ONLY:  message_string, turbulent_inflow, turbulent_outflow
566
567    IMPLICIT NONE
568
569!
570!-- Define which turbulence closure is going to be used
571    SELECT CASE ( TRIM( turbulence_closure ) )
572
573       CASE ( 'dynamic' )
574          les_dynamic = .TRUE.
575
576       CASE ( 'Moeng_Wyngaard' )
577          les_mw = .TRUE.
578
579       CASE ( 'TKE-l' )
580          rans_tke_l = .TRUE.
581          rans_mode = .TRUE.
582
583       CASE ( 'TKE-e' )
584          rans_tke_e = .TRUE.
585          rans_mode = .TRUE.
586
587       CASE DEFAULT
588          message_string = 'Unknown turbulence closure: ' //                &
589                           TRIM( turbulence_closure )
590          CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
591
592    END SELECT
593!
594!-- Set variables for RANS mode or LES mode
595    IF ( rans_mode )  THEN
596!
597!--    Assign values to constants for RANS mode
598       dsig_e    = 1.0_wp / rans_const_sigma(1)
599       dsig_diss = 1.0_wp / rans_const_sigma(2)
600
601       c_0 = rans_const_c(0)
602       c_1 = rans_const_c(1)
603       c_2 = rans_const_c(2)
604       c_3 = rans_const_c(3)   !> @todo clarify how to switch between different models
605       c_4 = rans_const_c(4)
606
607       IF ( turbulent_inflow .OR. turbulent_outflow )  THEN
608          message_string = 'turbulent inflow/outflow is not yet '//            &
609                           'implemented for RANS mode'
610          CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 )
611       ENDIF
612
613    ELSE
614!
615!--    LES mode
616       c_0 = 0.1_wp    !according to Lilly (1967) and Deardorff (1980)
617
618       dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead
619                       !of K_e which is used in RANS mode
620
621    ENDIF
622
623 END SUBROUTINE tcm_check_parameters
624
625!------------------------------------------------------------------------------!
626! Description:
627! ------------
628!> Check data output.
629!------------------------------------------------------------------------------!
630 SUBROUTINE tcm_check_data_output( var, unit )
631
632    IMPLICIT NONE
633
634    CHARACTER (LEN=*) ::  unit     !< unit of output variable
635    CHARACTER (LEN=*) ::  var      !< name of output variable
636
637
638    SELECT CASE ( TRIM( var ) )
639
640       CASE ( 'diss' )
641          unit = 'm2/s3'
642
643       CASE ( 'kh', 'km' )
644          unit = 'm2/s'
645
646       CASE DEFAULT
647          unit = 'illegal'
648
649    END SELECT
650
651 END SUBROUTINE tcm_check_data_output
652
653
654!------------------------------------------------------------------------------!
655! Description:
656! ------------
657!> Define appropriate grid for netcdf variables.
658!> It is called out from subroutine netcdf.
659!------------------------------------------------------------------------------!
660 SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
661
662    IMPLICIT NONE
663
664    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
665    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
666    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
667    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< name of output variable
668
669    LOGICAL, INTENT(OUT) ::  found   !< flag if output variable is found
670
671    found  = .TRUE.
672
673!
674!-- Check for the grid
675    SELECT CASE ( TRIM( var ) )
676
677       CASE ( 'diss', 'diss_xy', 'diss_xz', 'diss_yz' )
678          grid_x = 'x'
679          grid_y = 'y'
680          grid_z = 'zu'
681
682       CASE ( 'kh', 'kh_xy', 'kh_xz', 'kh_yz' )
683          grid_x = 'x'
684          grid_y = 'y'
685          grid_z = 'zu'
686
687       CASE ( 'km', 'km_xy', 'km_xz', 'km_yz' )
688          grid_x = 'x'
689          grid_y = 'y'
690          grid_z = 'zu'
691
692       CASE DEFAULT
693          found  = .FALSE.
694          grid_x = 'none'
695          grid_y = 'none'
696          grid_z = 'none'
697
698    END SELECT
699
700 END SUBROUTINE tcm_define_netcdf_grid
701
702
703!------------------------------------------------------------------------------!
704! Description:
705! ------------
706!> Average 3D data.
707!------------------------------------------------------------------------------!
708 SUBROUTINE tcm_3d_data_averaging( mode, variable )
709
710
711    USE averaging,                                                             &
712        ONLY:  diss_av, kh_av, km_av
713
714    USE control_parameters,                                                    &
715        ONLY:  average_count_3d
716
717    IMPLICIT NONE
718
719    CHARACTER (LEN=*) ::  mode       !< flag defining mode 'allocate', 'sum' or 'average'
720    CHARACTER (LEN=*) ::  variable   !< name of variable
721
722    INTEGER(iwp) ::  i   !< loop index
723    INTEGER(iwp) ::  j   !< loop index
724    INTEGER(iwp) ::  k   !< loop index
725
726    IF ( mode == 'allocate' )  THEN
727
728       SELECT CASE ( TRIM( variable ) )
729
730          CASE ( 'diss' )
731             IF ( .NOT. ALLOCATED( diss_av ) )  THEN
732                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
733             ENDIF
734             diss_av = 0.0_wp
735
736          CASE ( 'kh' )
737             IF ( .NOT. ALLOCATED( kh_av ) )  THEN
738                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
739             ENDIF
740             kh_av = 0.0_wp
741
742          CASE ( 'km' )
743             IF ( .NOT. ALLOCATED( km_av ) )  THEN
744                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
745             ENDIF
746             km_av = 0.0_wp
747
748          CASE DEFAULT
749             CONTINUE
750
751       END SELECT
752
753    ELSEIF ( mode == 'sum' )  THEN
754
755       SELECT CASE ( TRIM( variable ) )
756
757          CASE ( 'diss' )
758             IF ( ALLOCATED( diss_av ) ) THEN
759                DO  i = nxlg, nxrg
760                   DO  j = nysg, nyng
761                      DO  k = nzb, nzt+1
762                         diss_av(k,j,i) = diss_av(k,j,i) + diss(k,j,i)
763                      ENDDO
764                   ENDDO
765                ENDDO
766             ENDIF
767
768          CASE ( 'kh' )
769             IF ( ALLOCATED( kh_av ) ) THEN
770                DO  i = nxlg, nxrg
771                   DO  j = nysg, nyng
772                      DO  k = nzb, nzt+1
773                         kh_av(k,j,i) = kh_av(k,j,i) + kh(k,j,i)
774                      ENDDO
775                   ENDDO
776                ENDDO
777             ENDIF
778
779          CASE ( 'km' )
780             IF ( ALLOCATED( km_av ) ) THEN
781                DO  i = nxlg, nxrg
782                   DO  j = nysg, nyng
783                      DO  k = nzb, nzt+1
784                         km_av(k,j,i) = km_av(k,j,i) + km(k,j,i)
785                      ENDDO
786                   ENDDO
787                ENDDO
788             ENDIF
789
790          CASE DEFAULT
791             CONTINUE
792
793       END SELECT
794
795    ELSEIF ( mode == 'average' )  THEN
796
797       SELECT CASE ( TRIM( variable ) )
798
799          CASE ( 'diss' )
800             IF ( ALLOCATED( diss_av ) ) THEN
801                DO  i = nxlg, nxrg
802                   DO  j = nysg, nyng
803                      DO  k = nzb, nzt+1
804                         diss_av(k,j,i) = diss_av(k,j,i)                       &
805                                        / REAL( average_count_3d, KIND=wp )
806                      ENDDO
807                   ENDDO
808                ENDDO
809             ENDIF
810
811          CASE ( 'kh' )
812             IF ( ALLOCATED( kh_av ) ) THEN
813                DO  i = nxlg, nxrg
814                   DO  j = nysg, nyng
815                      DO  k = nzb, nzt+1
816                         kh_av(k,j,i) = kh_av(k,j,i)                           &
817                                        / REAL( average_count_3d, KIND=wp )
818                      ENDDO
819                   ENDDO
820                ENDDO
821             ENDIF
822
823          CASE ( 'km' )
824             IF ( ALLOCATED( km_av ) ) THEN
825                DO  i = nxlg, nxrg
826                   DO  j = nysg, nyng
827                      DO  k = nzb, nzt+1
828                         km_av(k,j,i) = km_av(k,j,i)                           &
829                                        / REAL( average_count_3d, KIND=wp )
830                      ENDDO
831                   ENDDO
832                ENDDO
833             ENDIF
834
835       END SELECT
836
837    ENDIF
838
839 END SUBROUTINE tcm_3d_data_averaging
840
841
842!------------------------------------------------------------------------------!
843! Description:
844! ------------
845!> Define 2D output variables.
846!------------------------------------------------------------------------------!
847 SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
848                                nzb_do, nzt_do )
849
850    USE averaging,                                                             &
851        ONLY:  diss_av, kh_av, km_av
852
853    IMPLICIT NONE
854
855    CHARACTER (LEN=*) ::  grid       !< name of vertical grid
856    CHARACTER (LEN=*) ::  mode       !< either 'xy', 'xz' or 'yz'
857    CHARACTER (LEN=*) ::  variable   !< name of variable
858
859    INTEGER(iwp) ::  av        !< flag for (non-)average output
860    INTEGER(iwp) ::  flag_nr   !< number of masking flag
861    INTEGER(iwp) ::  i         !< loop index
862    INTEGER(iwp) ::  j         !< loop index
863    INTEGER(iwp) ::  k         !< loop index
864    INTEGER(iwp) ::  nzb_do    !< vertical output index (bottom)
865    INTEGER(iwp) ::  nzt_do    !< vertical output index (top)
866
867    LOGICAL ::  found     !< flag if output variable is found
868    LOGICAL ::  resorted  !< flag if output is already resorted
869
870    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
871
872    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< local
873       !< array to which output data is resorted to
874
875    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
876
877    found = .TRUE.
878    resorted = .FALSE.
879!
880!-- Set masking flag for topography for not resorted arrays
881    flag_nr = 0
882
883    SELECT CASE ( TRIM( variable ) )
884
885       CASE ( 'diss_xy', 'diss_xz', 'diss_yz' )
886          IF ( av == 0 )  THEN
887             to_be_resorted => diss
888          ELSE
889             IF ( .NOT. ALLOCATED( diss_av ) ) THEN
890                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
891                diss_av = REAL( fill_value, KIND = wp )
892             ENDIF
893             to_be_resorted => diss_av
894          ENDIF
895          IF ( mode == 'xy' ) grid = 'zu'
896
897       CASE ( 'kh_xy', 'kh_xz', 'kh_yz' )
898          IF ( av == 0 )  THEN
899             to_be_resorted => kh
900          ELSE
901             IF ( .NOT. ALLOCATED( kh_av ) ) THEN
902                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
903                kh_av = REAL( fill_value, KIND = wp )
904             ENDIF
905             to_be_resorted => kh_av
906          ENDIF
907          IF ( mode == 'xy' ) grid = 'zu'
908
909       CASE ( 'km_xy', 'km_xz', 'km_yz' )
910          IF ( av == 0 )  THEN
911             to_be_resorted => km
912          ELSE
913             IF ( .NOT. ALLOCATED( km_av ) ) THEN
914                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
915                km_av = REAL( fill_value, KIND = wp )
916             ENDIF
917             to_be_resorted => km_av
918          ENDIF
919          IF ( mode == 'xy' ) grid = 'zu'
920
921       CASE DEFAULT
922          found = .FALSE.
923          grid  = 'none'
924
925    END SELECT
926
927    IF ( found .AND. .NOT. resorted )  THEN
928       DO  i = nxl, nxr
929          DO  j = nys, nyn
930             DO  k = nzb_do, nzt_do
931                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
932                                  REAL( fill_value, KIND = wp ),               &
933                                  BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
934             ENDDO
935          ENDDO
936       ENDDO
937    ENDIF
938
939 END SUBROUTINE tcm_data_output_2d
940
941
942!------------------------------------------------------------------------------!
943! Description:
944! ------------
945!> Define 3D output variables.
946!------------------------------------------------------------------------------!
947 SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
948
949
950    USE averaging,                                                             &
951        ONLY:  diss_av, kh_av, km_av
952
953    IMPLICIT NONE
954
955    CHARACTER (LEN=*) ::  variable   !< name of variable
956
957    INTEGER(iwp) ::  av        !< flag for (non-)average output
958    INTEGER(iwp) ::  flag_nr   !< number of masking flag
959    INTEGER(iwp) ::  i         !< loop index
960    INTEGER(iwp) ::  j         !< loop index
961    INTEGER(iwp) ::  k         !< loop index
962    INTEGER(iwp) ::  nzb_do    !< lower limit of the data output (usually 0)
963    INTEGER(iwp) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
964
965    LOGICAL ::  found     !< flag if output variable is found
966    LOGICAL ::  resorted  !< flag if output is already resorted
967
968    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
969
970    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
971       !< array to which output data is resorted to
972
973    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
974
975    found = .TRUE.
976    resorted = .FALSE.
977!
978!-- Set masking flag for topography for not resorted arrays
979    flag_nr = 0
980
981    SELECT CASE ( TRIM( variable ) )
982
983       CASE ( 'diss' )
984          IF ( av == 0 )  THEN
985             to_be_resorted => diss
986          ELSE
987             IF ( .NOT. ALLOCATED( diss_av ) ) THEN
988                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
989                diss_av = REAL( fill_value, KIND = wp )
990             ENDIF
991             to_be_resorted => diss_av
992          ENDIF
993
994       CASE ( 'kh' )
995          IF ( av == 0 )  THEN
996             to_be_resorted => kh
997          ELSE
998             IF ( .NOT. ALLOCATED( kh_av ) ) THEN
999                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1000                kh_av = REAL( fill_value, KIND = wp )
1001             ENDIF
1002             to_be_resorted => kh_av
1003          ENDIF
1004
1005       CASE ( 'km' )
1006          IF ( av == 0 )  THEN
1007             to_be_resorted => km
1008          ELSE
1009             IF ( .NOT. ALLOCATED( km_av ) ) THEN
1010                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1011                km_av = REAL( fill_value, KIND = wp )
1012             ENDIF
1013             to_be_resorted => km_av
1014          ENDIF
1015
1016       CASE DEFAULT
1017          found = .FALSE.
1018
1019    END SELECT
1020
1021
1022    IF ( found .AND. .NOT. resorted )  THEN
1023       DO  i = nxl, nxr
1024          DO  j = nys, nyn
1025             DO  k = nzb_do, nzt_do
1026                local_pf(i,j,k) = MERGE(                                 &
1027                                   to_be_resorted(k,j,i),                &
1028                                   REAL( fill_value, KIND = wp ),        &
1029                                   BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
1030             ENDDO
1031          ENDDO
1032       ENDDO
1033       resorted = .TRUE.
1034    ENDIF
1035
1036 END SUBROUTINE tcm_data_output_3d
1037
1038
1039!------------------------------------------------------------------------------!
1040! Description:
1041! ------------
1042!> Allocate arrays and assign pointers.
1043!------------------------------------------------------------------------------!
1044 SUBROUTINE tcm_init_arrays
1045
1046    USE bulk_cloud_model_mod,                                                  &
1047        ONLY:  collision_turbulence
1048
1049    USE pmc_interface,                                                         &
1050        ONLY:  nested_run
1051
1052    IMPLICIT NONE
1053
1054    ALLOCATE( kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1055    ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1056
1057    ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1058    ALLOCATE( e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1059    ALLOCATE( e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1060
1061!
1062!-- Allocate arrays required for dissipation.
1063!-- Please note, if it is a nested run, arrays need to be allocated even if
1064!-- they do not necessarily need to be transferred, which is attributed to
1065!-- the design of the model coupler which allocates memory for each variable.
1066    ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1067
1068    IF ( rans_tke_e  .OR.  nested_run )  THEN
1069       ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1070       ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1071    ENDIF
1072
1073!
1074!-- Initial assignment of pointers
1075    e  => e_1;   e_p  => e_2;   te_m  => e_3
1076
1077    diss => diss_1
1078    IF ( rans_tke_e  .OR.  nested_run )  THEN
1079       diss_p => diss_2; tdiss_m => diss_3
1080    ENDIF
1081
1082 END SUBROUTINE tcm_init_arrays
1083
1084
1085!------------------------------------------------------------------------------!
1086! Description:
1087! ------------
1088!> Initialization of turbulence closure module.
1089!------------------------------------------------------------------------------!
1090 SUBROUTINE tcm_init
1091
1092    USE control_parameters,                                                    &
1093        ONLY:  bc_dirichlet_l, complex_terrain, topography
1094
1095    USE model_1d_mod,                                                          &
1096        ONLY:  e1d, kh1d, km1d
1097
1098    IMPLICIT NONE
1099
1100    INTEGER(iwp) :: i            !< loop index
1101    INTEGER(iwp) :: j            !< loop index
1102    INTEGER(iwp) :: k            !< loop index
1103    INTEGER(iwp) :: nz_s_shift   !< lower shift index for scalars
1104    INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow
1105
1106!
1107!-- Initialize mixing length
1108    CALL tcm_init_mixing_length
1109
1110!
1111!-- Actions for initial runs
1112    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1113         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1114
1115       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1116
1117          IF ( .NOT. rans_tke_e ) THEN
1118!
1119!--          Transfer initial profiles to the arrays of the 3D model
1120             DO  i = nxlg, nxrg
1121                DO  j = nysg, nyng
1122                   e(:,j,i)  = e1d
1123                   kh(:,j,i) = kh1d
1124                   km(:,j,i) = km1d
1125                ENDDO
1126             ENDDO
1127
1128             IF ( constant_diffusion )  THEN
1129                e = 0.0_wp
1130             ENDIF
1131
1132             diss = 0.0_wp
1133
1134          ELSE
1135!
1136!--          In case of TKE-e closure in RANS mode, do not use e, diss, and km
1137!--          profiles from 1D model. Instead, initialize with constant profiles
1138             IF ( constant_diffusion )  THEN
1139                km = km_constant
1140                kh = km / prandtl_number
1141                e  = 0.0_wp
1142             ELSEIF ( e_init > 0.0_wp )  THEN
1143                DO  i = nxlg, nxrg
1144                   DO  j = nysg, nyng
1145                      DO  k = nzb+1, nzt
1146                         km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
1147                      ENDDO
1148                   ENDDO
1149                ENDDO
1150                km(nzb,:,:)   = km(nzb+1,:,:)
1151                km(nzt+1,:,:) = km(nzt,:,:)
1152                kh = km / prandtl_number
1153                e  = e_init
1154             ELSE
1155                IF ( .NOT. ocean_mode )  THEN
1156                   kh   = 0.01_wp   ! there must exist an initial diffusion, because
1157                   km   = 0.01_wp   ! otherwise no TKE would be produced by the
1158                                    ! production terms, as long as not yet
1159                                    ! e = (u*/cm)**2 at k=nzb+1
1160                ELSE
1161                   kh   = 0.00001_wp
1162                   km   = 0.00001_wp
1163                ENDIF
1164                e    = 0.0_wp
1165             ENDIF
1166
1167             DO  i = nxlg, nxrg
1168                DO  j = nysg, nyng
1169                   DO  k = nzb+1, nzt
1170                      diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i)
1171                   ENDDO
1172                ENDDO
1173             ENDDO
1174             diss(nzb,:,:) = diss(nzb+1,:,:)
1175             diss(nzt+1,:,:) = diss(nzt,:,:)
1176
1177          ENDIF
1178
1179       ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. &
1180                INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
1181
1182          IF ( constant_diffusion )  THEN
1183             km = km_constant
1184             kh = km / prandtl_number
1185             e  = 0.0_wp
1186          ELSEIF ( e_init > 0.0_wp )  THEN
1187             DO  i = nxlg, nxrg
1188                DO  j = nysg, nyng
1189                   DO  k = nzb+1, nzt
1190                      km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
1191                   ENDDO
1192                ENDDO
1193             ENDDO
1194             km(nzb,:,:)   = km(nzb+1,:,:)
1195             km(nzt+1,:,:) = km(nzt,:,:)
1196             kh = km / prandtl_number
1197             e  = e_init
1198          ELSE
1199             IF ( .NOT. ocean_mode )  THEN
1200                kh   = 0.01_wp   ! there must exist an initial diffusion, because
1201                km   = 0.01_wp   ! otherwise no TKE would be produced by the
1202                                 ! production terms, as long as not yet
1203                                 ! e = (u*/cm)**2 at k=nzb+1
1204             ELSE
1205                kh   = 0.00001_wp
1206                km   = 0.00001_wp
1207             ENDIF
1208             e    = 0.0_wp
1209          ENDIF
1210
1211          IF ( rans_tke_e )  THEN
1212             DO  i = nxlg, nxrg
1213                DO  j = nysg, nyng
1214                   DO  k = nzb+1, nzt
1215                      diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i)
1216                   ENDDO
1217                ENDDO
1218             ENDDO
1219             diss(nzb,:,:) = diss(nzb+1,:,:)
1220             diss(nzt+1,:,:) = diss(nzt,:,:)
1221          ELSE
1222             diss = 0.0_wp
1223          ENDIF
1224
1225       ENDIF
1226!
1227!--    Store initial profiles for output purposes etc.
1228       hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
1229       hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
1230!
1231!--    Initialize old and new time levels.
1232       te_m = 0.0_wp
1233       e_p = e
1234       IF ( rans_tke_e )  THEN
1235          tdiss_m = 0.0_wp
1236          diss_p = diss
1237       ENDIF
1238
1239    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
1240             TRIM( initializing_actions ) == 'cyclic_fill' )                   &
1241    THEN
1242
1243!
1244!--    In case of complex terrain and cyclic fill method as initialization,
1245!--    shift initial data in the vertical direction for each point in the
1246!--    x-y-plane depending on local surface height
1247       IF ( complex_terrain  .AND.                                             &
1248            TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1249          DO  i = nxlg, nxrg
1250             DO  j = nysg, nyng
1251                nz_s_shift = topo_top_ind(j,i,0)
1252
1253                e(nz_s_shift:nzt+1,j,i)  =  e(0:nzt+1-nz_s_shift,j,i)
1254                km(nz_s_shift:nzt+1,j,i) = km(0:nzt+1-nz_s_shift,j,i)
1255                kh(nz_s_shift:nzt+1,j,i) = kh(0:nzt+1-nz_s_shift,j,i)
1256             ENDDO
1257          ENDDO
1258          IF ( rans_tke_e )  THEN
1259             DO  i = nxlg, nxrg
1260                DO  j = nysg, nyng
1261                   nz_s_shift = topo_top_ind(j,i,0)
1262
1263                   diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1-nz_s_shift,j,i)
1264                ENDDO
1265             ENDDO
1266          ELSE
1267             diss = 0.0_wp
1268          ENDIF
1269       ENDIF
1270
1271!
1272!--    Initialization of the turbulence recycling method
1273       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
1274            turbulent_inflow )  THEN
1275          mean_inflow_profiles(:,5) = hom_sum(:,8,0)   ! e
1276!
1277!--       In case of complex terrain, determine vertical displacement at inflow
1278!--       boundary and adjust mean inflow profiles
1279          IF ( complex_terrain )  THEN
1280             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND.  &
1281                  nysg <= 0 .AND. nyng >= 0        )  THEN
1282                nz_s_shift_l = topo_top_ind(0,0,0)
1283             ELSE
1284                nz_s_shift_l = 0
1285             ENDIF
1286#if defined( __parallel )
1287             CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER,      &
1288                                MPI_MAX, comm2d, ierr)
1289#else
1290             nz_s_shift = nz_s_shift_l
1291#endif
1292             mean_inflow_profiles(nz_s_shift:nzt+1,5) =  &
1293                hom_sum(0:nzt+1-nz_s_shift,8,0)  ! e
1294          ENDIF
1295!
1296!--       Use these mean profiles at the inflow (provided that Dirichlet
1297!--       conditions are used)
1298          IF ( bc_dirichlet_l )  THEN
1299             DO  j = nysg, nyng
1300                DO  k = nzb, nzt+1
1301                   e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
1302                ENDDO
1303             ENDDO
1304          ENDIF
1305       ENDIF
1306!
1307!--    Inside buildings set TKE back to zero
1308       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
1309            topography /= 'flat' )  THEN
1310!
1311!--       Inside buildings set TKE back to zero.
1312!--       Other scalars (km, kh,...) are ignored at present,
1313!--       maybe revise later.
1314          DO  i = nxlg, nxrg
1315             DO  j = nysg, nyng
1316                DO  k = nzb, nzt
1317                   e(k,j,i)     = MERGE( e(k,j,i), 0.0_wp,                     &
1318                                         BTEST( wall_flags_total_0(k,j,i), 0 ) )
1319                ENDDO
1320             ENDDO
1321          ENDDO
1322
1323          IF ( rans_tke_e )  THEN
1324             DO  i = nxlg, nxrg
1325                DO  j = nysg, nyng
1326                   DO  k = nzb, nzt
1327                      diss(k,j,i)    = MERGE( diss(k,j,i), 0.0_wp,             &
1328                                         BTEST( wall_flags_total_0(k,j,i), 0 ) )
1329                   ENDDO
1330                ENDDO
1331             ENDDO
1332          ENDIF
1333       ENDIF
1334!
1335!--    Initialize new time levels (only done in order to set boundary values
1336!--    including ghost points)
1337       e_p = e
1338!
1339!--    Allthough tendency arrays are set in prognostic_equations, they have
1340!--    to be predefined here because there they are used (but multiplied with 0)
1341!--    before they are set.
1342       te_m = 0.0_wp
1343
1344       IF ( rans_tke_e )  THEN
1345          diss_p = diss
1346          tdiss_m = 0.0_wp
1347       ENDIF
1348
1349    ENDIF
1350
1351 END SUBROUTINE tcm_init
1352
1353
1354!------------------------------------------------------------------------------!
1355! Description:
1356! ------------
1357!> Pre-computation of grid-dependent and near-wall mixing length.
1358!> @todo consider walls in horizontal direction at a distance further than a
1359!>       single grid point (RANS mode)
1360!------------------------------------------------------------------------------!
1361 SUBROUTINE tcm_init_mixing_length
1362
1363    USE arrays_3d,                                                             &
1364        ONLY:  dzw, ug, vg, zu, zw
1365
1366    USE control_parameters,                                                    &
1367        ONLY:  f, message_string, wall_adjustment, wall_adjustment_factor
1368
1369    USE grid_variables,                                                        &
1370        ONLY:  dx, dy
1371
1372    USE indices,                                                               &
1373        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
1374               nzt, wall_flags_total_0
1375
1376    USE kinds
1377
1378
1379    IMPLICIT NONE
1380
1381    INTEGER(iwp) :: dist_dx     !< found distance devided by dx
1382    INTEGER(iwp) :: i           !< index variable along x
1383    INTEGER(iwp) :: ii          !< index variable along x
1384    INTEGER(iwp) :: j           !< index variable along y
1385    INTEGER(iwp) :: k           !< index variable along z
1386    INTEGER(iwp) :: k_max_topo  !< index of maximum topography height
1387    INTEGER(iwp) :: kk          !< index variable along z
1388    INTEGER(iwp) :: rad_i       !< search radius in grid points along x
1389    INTEGER(iwp) :: rad_i_l     !< possible search radius to the left
1390    INTEGER(iwp) :: rad_i_r     !< possible search radius to the right
1391    INTEGER(iwp) :: rad_j       !< search radius in grid points along y
1392    INTEGER(iwp) :: rad_j_n     !< possible search radius to north
1393    INTEGER(iwp) :: rad_j_s     !< possible search radius to south
1394    INTEGER(iwp) :: rad_k       !< search radius in grid points along z
1395    INTEGER(iwp) :: rad_k_b     !< search radius in grid points along negative z
1396    INTEGER(iwp) :: rad_k_t     !< search radius in grid points along positive z
1397
1398    INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
1399
1400    INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
1401
1402    REAL(wp) :: radius          !< search radius in meter
1403
1404    ALLOCATE( l_grid(1:nzt) )
1405    ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1406!
1407!-- Initialize the mixing length in case of an LES-simulation
1408    IF ( .NOT. rans_mode )  THEN
1409!
1410!--    Compute the grid-dependent mixing length.
1411       DO  k = 1, nzt
1412          l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
1413       ENDDO
1414!
1415!--    Initialize near-wall mixing length l_wall only in the vertical direction
1416!--    for the moment, multiplication with wall_adjustment_factor further below
1417       l_wall(nzb,:,:)   = l_grid(1)
1418       DO  k = nzb+1, nzt
1419          l_wall(k,:,:)  = l_grid(k)
1420       ENDDO
1421       l_wall(nzt+1,:,:) = l_grid(nzt)
1422
1423       IF ( wall_adjustment )  THEN
1424
1425          DO  k = 1, nzt
1426             IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.            &
1427                  l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
1428                WRITE( message_string, * ) 'grid anisotropy exceeds ',             &
1429                                           'threshold given by only local',        &
1430                                           ' &horizontal reduction of near_wall ', &
1431                                           'mixing length l_wall',                 &
1432                                           ' &starting from height level k = ', k, &
1433                                           '.'
1434                CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
1435                EXIT
1436             ENDIF
1437          ENDDO
1438!
1439!--       In case of topography: limit near-wall mixing length l_wall further:
1440!--       Go through all points of the subdomain one by one and look for the closest
1441!--       surface.
1442!--       Is this correct in the ocean case?
1443          DO  i = nxl, nxr
1444             DO  j = nys, nyn
1445                DO  k = nzb+1, nzt
1446!
1447!--                Check if current gridpoint belongs to the atmosphere
1448                   IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
1449!
1450!--                   Check for neighbouring grid-points.
1451!--                   Vertical distance, down
1452                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i), 0 ) )             &
1453                         l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
1454!
1455!--                   Vertical distance, up
1456                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )             &
1457                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), zw(k) - zu(k) )
1458!
1459!--                   y-distance
1460                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 )  .OR.         &
1461                           .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) )             &
1462                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
1463!
1464!--                   x-distance
1465                      IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 )  .OR.         &
1466                           .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) )             &
1467                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
1468!
1469!--                   yz-distance (vertical edges, down)
1470                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 )  .OR.       &
1471                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i), 0 )  )          &
1472                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
1473                                              SQRT( 0.25_wp * dy**2 +                  &
1474                                             ( zu(k) - zw(k-1) )**2 ) )
1475!
1476!--                   yz-distance (vertical edges, up)
1477                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 )  .OR.       &
1478                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i), 0 )  )          &
1479                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
1480                                              SQRT( 0.25_wp * dy**2 +                  &
1481                                             ( zw(k) - zu(k) )**2 ) )
1482!
1483!--                   xz-distance (vertical edges, down)
1484                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 )  .OR.       &
1485                           .NOT. BTEST( wall_flags_total_0(k-1,j,i+1), 0 )  )          &
1486                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
1487                                              SQRT( 0.25_wp * dx**2 +                  &
1488                                             ( zu(k) - zw(k-1) )**2 ) )
1489!
1490!--                   xz-distance (vertical edges, up)
1491                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 )  .OR.       &
1492                           .NOT. BTEST( wall_flags_total_0(k+1,j,i+1), 0 )  )          &
1493                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
1494                                              SQRT( 0.25_wp * dx**2 +                  &
1495                                             ( zw(k) - zu(k) )**2 ) )
1496!
1497!--                   xy-distance (horizontal edges)
1498                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 )  .OR.        &
1499                           .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 )  .OR.        &
1500                           .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 )  .OR.        &
1501                           .NOT. BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) )            &
1502                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
1503                                              SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) )
1504!
1505!--                   xyz distance (vertical and horizontal edges, down)
1506                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 )  .OR.      &
1507                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 )  .OR.      &
1508                           .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 )  .OR.      &
1509                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i+1), 0 ) )          &
1510                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
1511                                              SQRT( 0.25_wp * ( dx**2 + dy**2 )         &
1512                                                    +  ( zu(k) - zw(k-1) )**2  ) )
1513!
1514!--                   xyz distance (vertical and horizontal edges, up)
1515                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 )  .OR.      &
1516                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 )  .OR.      &
1517                           .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 )  .OR.      &
1518                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i+1), 0 ) )          &
1519                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
1520                                              SQRT( 0.25_wp * ( dx**2 + dy**2 )         &
1521                                                    +  ( zw(k) - zu(k) )**2  ) )
1522
1523                   ENDIF
1524!
1525!--                Adjust mixing length by wall-adjustment factor and limit it by l_grid
1526                   l_wall(k,j,i) = MIN( l_wall(k,j,i) * wall_adjustment_factor, l_grid(k) )
1527
1528                ENDDO  !k loop
1529             ENDDO  !j loop
1530          ENDDO  !i loop
1531
1532       ENDIF  !if wall_adjustment
1533
1534    ELSE
1535!
1536!--    Initialize the mixing length in case of a RANS simulation
1537       ALLOCATE( l_black(nzb:nzt+1) )
1538
1539!
1540!--    Calculate mixing length according to Blackadar (1962)
1541       IF ( f /= 0.0_wp )  THEN
1542          l_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /            &
1543                  ABS( f ) + 1.0E-10_wp
1544       ELSE
1545          l_max = 30.0_wp
1546       ENDIF
1547
1548       DO  k = nzb, nzt
1549          l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / l_max )
1550       ENDDO
1551
1552       l_black(nzt+1) = l_black(nzt)
1553
1554!
1555!--    Get height level of highest topography within local subdomain
1556       k_max_topo = 0
1557       DO  i = nxlg, nxrg
1558          DO  j = nysg, nyng
1559             DO  k = nzb+1, nzt-1
1560                IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 0 ) .AND.  &
1561                     k > k_max_topo )  &
1562                   k_max_topo = k
1563             ENDDO
1564          ENDDO
1565       ENDDO
1566
1567       l_wall(nzb,:,:) = l_black(nzb)
1568       l_wall(nzt+1,:,:) = l_black(nzt+1)
1569!
1570!--    Limit mixing length to either nearest wall or Blackadar mixing length.
1571!--    For that, analyze each grid point (i/j/k) ("analysed grid point") and
1572!--    search within its vicinity for the shortest distance to a wall by cal-
1573!--    culating the distance between the analysed grid point and the "viewed
1574!--    grid point" if it contains a wall (belongs to topography).
1575       DO  k = nzb+1, nzt
1576
1577          radius = l_black(k)  ! radius within walls are searched
1578!
1579!--       Set l_wall to its default maximum value (l_back)
1580          l_wall(k,:,:) = radius
1581
1582!
1583!--       Compute search radius as number of grid points in all directions
1584          rad_i = CEILING( radius / dx )
1585          rad_j = CEILING( radius / dy )
1586
1587          DO  kk = 0, nzt-k
1588             rad_k_t = kk
1589!
1590!--          Limit upward search radius to height of maximum topography
1591             IF ( zu(k+kk)-zu(k) >= radius .OR. k+kk >= k_max_topo )  EXIT
1592          ENDDO
1593
1594          DO  kk = 0, k
1595             rad_k_b = kk
1596             IF ( zu(k)-zu(k-kk) >= radius )  EXIT
1597          ENDDO
1598
1599!
1600!--       Get maximum vertical radius; necessary for defining arrays
1601          rad_k = MAX( rad_k_b, rad_k_t )
1602!
1603!--       When analysed grid point lies above maximum topography, set search
1604!--       radius to 0 if the distance between the analysed grid point and max
1605!--       topography height is larger than the maximum search radius
1606          IF ( zu(k-rad_k_b) > zu(k_max_topo) )  rad_k_b = 0
1607!
1608!--       Search within vicinity only if the vertical search radius is >0
1609          IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 )  THEN
1610
1611             !> @note shape of vicinity is larger in z direction
1612             !>   Shape of vicinity is two grid points larger than actual search
1613             !>   radius in vertical direction. The first and last grid point is
1614             !>   always set to 1 to asure correct detection of topography. See
1615             !>   function "shortest_distance" for details.
1616             !>   2018-03-16, gronemeier
1617             ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) )
1618             ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) )
1619
1620             vicinity = 1
1621
1622             DO  i = nxl, nxr
1623                DO  j = nys, nyn
1624!
1625!--                Start search only if (i/j/k) belongs to atmosphere
1626                   IF ( BTEST( wall_flags_total_0(k,j,i), 0 )  )  THEN
1627!
1628!--                   Reset topography within vicinity
1629                      vicinity(-rad_k:rad_k,:,:) = 0
1630!
1631!--                   Copy area surrounding analysed grid point into vicinity.
1632!--                   First, limit size of data copied to vicinity by the domain
1633!--                   border
1634                      !> @note limit copied area to 1 grid point in hor. dir.
1635                      !>   Ignore walls in horizontal direction which are
1636                      !>   further away than a single grid point. This allows to
1637                      !>   only search within local subdomain without the need
1638                      !>   of global topography information.
1639                      !>   The error made by this assumption are acceptable at
1640                      !>   the moment.
1641                      !>   2018-10-01, gronemeier
1642                      rad_i_l = MIN( 1, rad_i, i )
1643                      rad_i_r = MIN( 1, rad_i, nx-i )
1644
1645                      rad_j_s = MIN( 1, rad_j, j )
1646                      rad_j_n = MIN( 1, rad_j, ny-j )
1647
1648                      CALL copy_into_vicinity( k, j, i,           &
1649                                               -rad_k_b, rad_k_t, &
1650                                               -rad_j_s, rad_j_n, &
1651                                               -rad_i_l, rad_i_r  )
1652                      !> @note in case of cyclic boundaries, those parts of the
1653                      !>   topography which lies beyond the domain borders but
1654                      !>   still within the search radius still needs to be
1655                      !>   copied into 'vicinity'. As the effective search
1656                      !>   radius is limited to 1 at the moment, no further
1657                      !>   copying is needed. Old implementation (prior to
1658                      !>   2018-10-01) had this covered but used a global array.
1659                      !>   2018-10-01, gronemeier
1660
1661!
1662!--                   Search for walls only if there is any within vicinity
1663                      IF ( MAXVAL( vicinity(-rad_k:rad_k,:,:) ) /= 0 )  THEN
1664!
1665!--                      Search within first half (positive x)
1666                         dist_dx = rad_i
1667                         DO  ii = 0, dist_dx
1668!
1669!--                         Search along vertical direction only if below
1670!--                         maximum topography
1671                            IF ( rad_k_t > 0 ) THEN
1672!
1673!--                            Search for walls within octant (+++)
1674                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1675                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
1676                                       shortest_distance( vic_yz, .TRUE., ii ) )
1677!
1678!--                            Search for walls within octant (+-+)
1679!--                            Switch order of array so that the analysed grid
1680!--                            point is always located at (0/0) (required by
1681!--                            shortest_distance").
1682                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1683                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
1684                                       shortest_distance( vic_yz, .TRUE., ii ) )
1685
1686                            ENDIF
1687!
1688!--                         Search for walls within octant (+--)
1689                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1690                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
1691                                      shortest_distance( vic_yz, .FALSE., ii ) )
1692!
1693!--                         Search for walls within octant (++-)
1694                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1695                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
1696                                      shortest_distance( vic_yz, .FALSE., ii ) )
1697!
1698!--                         Reduce search along x by already found distance
1699                            dist_dx = CEILING( l_wall(k,j,i) / dx )
1700
1701                         ENDDO
1702!
1703!-                       Search within second half (negative x)
1704                         DO  ii = 0, -dist_dx, -1
1705!
1706!--                         Search along vertical direction only if below
1707!--                         maximum topography
1708                            IF ( rad_k_t > 0 ) THEN
1709!
1710!--                            Search for walls within octant (-++)
1711                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1712                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
1713                                      shortest_distance( vic_yz, .TRUE., -ii ) )
1714!
1715!--                            Search for walls within octant (--+)
1716!--                            Switch order of array so that the analysed grid
1717!--                            point is always located at (0/0) (required by
1718!--                            shortest_distance").
1719                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1720                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
1721                                      shortest_distance( vic_yz, .TRUE., -ii ) )
1722
1723                            ENDIF
1724!
1725!--                         Search for walls within octant (---)
1726                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1727                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
1728                                     shortest_distance( vic_yz, .FALSE., -ii ) )
1729!
1730!--                         Search for walls within octant (-+-)
1731                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1732                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
1733                                     shortest_distance( vic_yz, .FALSE., -ii ) )
1734!
1735!--                         Reduce search along x by already found distance
1736                            dist_dx = CEILING( l_wall(k,j,i) / dx )
1737
1738                         ENDDO
1739
1740                      ENDIF  !Check for any walls within vicinity
1741
1742                   ELSE  !Check if (i,j,k) belongs to atmosphere
1743
1744                      l_wall(k,j,i) = l_black(k)
1745
1746                   ENDIF
1747
1748                ENDDO  !j loop
1749             ENDDO  !i loop
1750
1751             DEALLOCATE( vicinity )
1752             DEALLOCATE( vic_yz )
1753
1754          ENDIF  !check vertical size of vicinity
1755
1756       ENDDO  !k loop
1757
1758       !$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1))
1759
1760    ENDIF  !LES or RANS mode
1761
1762!
1763!-- Set lateral boundary conditions for l_wall
1764    CALL exchange_horiz( l_wall, nbgp )
1765
1766    !$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) &
1767    !$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
1768
1769    CONTAINS
1770!------------------------------------------------------------------------------!
1771! Description:
1772! ------------
1773!> Calculate the shortest distance between position (i/j/k)=(0/0/0) and
1774!> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
1775!> closest to the origin (0/0) of 'array'.
1776!------------------------------------------------------------------------------!
1777    REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i )
1778
1779       IMPLICIT NONE
1780
1781       LOGICAL, INTENT(IN) :: orientation    !< flag if array represents an array oriented upwards (true) or downwards (false)
1782
1783       INTEGER(iwp), INTENT(IN) :: pos_i     !< x position of the yz-plane 'array'
1784
1785       INTEGER(iwp) :: a                     !< loop index
1786       INTEGER(iwp) :: b                     !< loop index
1787       INTEGER(iwp) :: jj                    !< loop index
1788
1789       INTEGER(KIND=1) :: maximum            !< maximum of array along z dimension
1790
1791       INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
1792
1793       INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yz-plane at position pos_i
1794
1795!
1796!--    Get coordinate of first maximum along vertical dimension
1797!--    at each y position of array (similar to function maxloc but more stable).
1798       DO  a = 0, rad_j
1799          loc_k(a) = rad_k+1
1800          maximum = MAXVAL( array(:,a) )
1801          DO  b = 0, rad_k+1
1802             IF ( array(b,a) == maximum )  THEN
1803                loc_k(a) = b
1804                EXIT
1805             ENDIF
1806          ENDDO
1807       ENDDO
1808!
1809!--    Set distance to the default maximum value (=search radius)
1810       shortest_distance = radius
1811!
1812!--    Calculate distance between position (0/0/0) and
1813!--    position (pos_i/jj/loc(jj)) and only save the shortest distance.
1814       IF ( orientation ) THEN  !if array is oriented upwards
1815          DO  jj = 0, rad_j
1816             shortest_distance =                                               &
1817                MIN( shortest_distance,                                        &
1818                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
1819                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
1820                         + MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2             &
1821                         )                                                     &
1822                   )
1823          ENDDO
1824       ELSE  !if array is oriented downwards
1825          !> @note MAX within zw required to circumvent error at domain border
1826          !>   At the domain border, if non-cyclic boundary is present, the
1827          !>   index for zw could be -1, which will be errorneous (zw(-1) does
1828          !>   not exist). The MAX function limits the index to be at least 0.
1829          DO  jj = 0, rad_j
1830             shortest_distance =                                               &
1831                MIN( shortest_distance,                                        &
1832                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
1833                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
1834                         + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2    &
1835                         )                                                     &
1836                   )
1837          ENDDO
1838       ENDIF
1839
1840    END FUNCTION
1841
1842!------------------------------------------------------------------------------!
1843! Description:
1844! ------------
1845!> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point
1846!> (kp,jp,ip) containing the first bit of wall_flags_total_0 into the array
1847!> 'vicinity'. Only copy first bit as this indicates the presence of topography.
1848!------------------------------------------------------------------------------!
1849    SUBROUTINE copy_into_vicinity( kp, jp, ip, kb, kt, js, jn, il, ir )
1850
1851       IMPLICIT NONE
1852
1853       INTEGER(iwp), INTENT(IN) :: il !< left loop boundary
1854       INTEGER(iwp), INTENT(IN) :: ip !< center position in x-direction
1855       INTEGER(iwp), INTENT(IN) :: ir !< right loop boundary
1856       INTEGER(iwp), INTENT(IN) :: jn !< northern loop boundary
1857       INTEGER(iwp), INTENT(IN) :: jp !< center position in y-direction
1858       INTEGER(iwp), INTENT(IN) :: js !< southern loop boundary
1859       INTEGER(iwp), INTENT(IN) :: kb !< bottom loop boundary
1860       INTEGER(iwp), INTENT(IN) :: kp !< center position in z-direction
1861       INTEGER(iwp), INTENT(IN) :: kt !< top loop boundary
1862
1863       INTEGER(iwp) :: i   !< loop index
1864       INTEGER(iwp) :: j   !< loop index
1865       INTEGER(iwp) :: k   !< loop index
1866
1867       DO  i = il, ir
1868          DO  j = js, jn
1869             DO  k = kb, kt
1870                vicinity(k,j,i) = MERGE( 0, 1,               &
1871                       BTEST( wall_flags_total_0(kp+k,jp+j,ip+i), 0 ) )
1872             ENDDO
1873          ENDDO
1874       ENDDO
1875
1876    END SUBROUTINE copy_into_vicinity
1877
1878 END SUBROUTINE tcm_init_mixing_length
1879
1880
1881!------------------------------------------------------------------------------!
1882! Description:
1883! ------------
1884!> Initialize virtual velocities used later in production_e.
1885!------------------------------------------------------------------------------!
1886 SUBROUTINE production_e_init
1887
1888    USE arrays_3d,                                                             &
1889        ONLY:  drho_air_zw, zu
1890
1891    USE control_parameters,                                                    &
1892        ONLY:  constant_flux_layer
1893
1894    USE surface_layer_fluxes_mod,                                              &
1895        ONLY:  phi_m
1896
1897    IMPLICIT NONE
1898
1899    INTEGER(iwp) ::  i      !< grid index x-direction
1900    INTEGER(iwp) ::  j      !< grid index y-direction
1901    INTEGER(iwp) ::  k      !< grid index z-direction
1902    INTEGER(iwp) ::  m      !< running index surface elements
1903
1904    REAL(wp) ::  km_sfc     !< diffusion coefficient, used to compute virtual velocities
1905
1906    IF ( constant_flux_layer )  THEN
1907!
1908!--    Calculate a virtual velocity at the surface in a way that the
1909!--    vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
1910!--    Prandtl law (-w'u'/km). This gradient is used in the TKE shear
1911!--    production term at k=1 (see production_e_ij).
1912!--    The velocity gradient has to be limited in case of too small km
1913!--    (otherwise the timestep may be significantly reduced by large
1914!--    surface winds).
1915!--    not available in case of non-cyclic boundary conditions.
1916!--    Default surfaces, upward-facing
1917       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
1918       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
1919       !$ACC PRESENT(surf_def_h(0), u, v, drho_air_zw, zu)
1920       DO  m = 1, surf_def_h(0)%ns
1921
1922          i = surf_def_h(0)%i(m)
1923          j = surf_def_h(0)%j(m)
1924          k = surf_def_h(0)%k(m)
1925!
1926!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1927!--       and km are not on the same grid. Actually, a further
1928!--       interpolation of km onto the u/v-grid is necessary. However, the
1929!--       effect of this error is negligible.
1930          km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) /       &
1931                   phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
1932
1933          surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *          &
1934                                     drho_air_zw(k-1)               *          &
1935                                     ( zu(k+1) - zu(k-1)    )       /          &
1936                                     ( km_sfc  + 1.0E-20_wp )
1937          surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *          &
1938                                     drho_air_zw(k-1)               *          &
1939                                     ( zu(k+1) - zu(k-1)    )       /          &
1940                                     ( km_sfc  + 1.0E-20_wp )
1941
1942          IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
1943               ABS( u(k+1,j,i) - u(k-1,j,i)           )                        &
1944             )  surf_def_h(0)%u_0(m) = u(k-1,j,i)
1945
1946          IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) )  >                     &
1947               ABS( v(k+1,j,i) - v(k-1,j,i)           )                        &
1948             )  surf_def_h(0)%v_0(m) = v(k-1,j,i)
1949
1950       ENDDO
1951!
1952!--    Default surfaces, downward-facing surfaces
1953       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
1954       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
1955       !$ACC PRESENT(surf_def_h(1), u, v, drho_air_zw, zu, km)
1956       DO  m = 1, surf_def_h(1)%ns
1957
1958          i = surf_def_h(1)%i(m)
1959          j = surf_def_h(1)%j(m)
1960          k = surf_def_h(1)%k(m)
1961!
1962!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1963!--       and km are not on the same grid. Actually, a further
1964!--       interpolation of km onto the u/v-grid is necessary. However, the
1965!--       effect of this error is negligible.
1966          surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *          &
1967                                     drho_air_zw(k-1) *                        &
1968                                     ( zu(k+1)    - zu(k-1)    )  /            &
1969                                     ( km(k,j,i)  + 1.0E-20_wp )
1970          surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *          &
1971                                     drho_air_zw(k-1) *                        &
1972                                     ( zu(k+1)    - zu(k-1)    )  /            &
1973                                     ( km(k,j,i)  + 1.0E-20_wp )
1974
1975          IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >                     &
1976               ABS( u(k+1,j,i)           - u(k-1,j,i) )                        &
1977             )  surf_def_h(1)%u_0(m) = u(k+1,j,i)
1978
1979          IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >                     &
1980               ABS( v(k+1,j,i)           - v(k-1,j,i) )                        &
1981             )  surf_def_h(1)%v_0(m) = v(k+1,j,i)
1982
1983       ENDDO
1984!
1985!--    Natural surfaces, upward-facing
1986       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
1987       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
1988       !$ACC PRESENT(surf_lsm_h, u, v, drho_air_zw, zu)
1989       DO  m = 1, surf_lsm_h%ns
1990
1991          i = surf_lsm_h%i(m)
1992          j = surf_lsm_h%j(m)
1993          k = surf_lsm_h%k(m)
1994!
1995!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1996!--       and km are not on the same grid. Actually, a further
1997!--       interpolation of km onto the u/v-grid is necessary. However, the
1998!--       effect of this error is negligible.
1999          km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) /             &
2000                   phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
2001
2002          surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)    *             &
2003                                        drho_air_zw(k-1)         *             &
2004                                        ( zu(k+1) - zu(k-1)    ) /             &
2005                                        ( km_sfc  + 1.0E-20_wp )
2006          surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)    *             &
2007                                        drho_air_zw(k-1)         *             &
2008                                        ( zu(k+1) - zu(k-1)    ) /             &
2009                                        ( km_sfc  + 1.0E-20_wp )
2010
2011          IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
2012               ABS( u(k+1,j,i) - u(k-1,j,i)   )                                &
2013             )  surf_lsm_h%u_0(m) = u(k-1,j,i)
2014
2015          IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) )  >                        &
2016               ABS( v(k+1,j,i) - v(k-1,j,i)   )                                &
2017             )  surf_lsm_h%v_0(m) = v(k-1,j,i)
2018
2019       ENDDO
2020!
2021!--    Urban surfaces, upward-facing
2022       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
2023       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
2024       !$ACC PRESENT(surf_usm_h, u, v, drho_air_zw, zu)
2025       DO  m = 1, surf_usm_h%ns
2026
2027          i = surf_usm_h%i(m)
2028          j = surf_usm_h%j(m)
2029          k = surf_usm_h%k(m)
2030!
2031!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
2032!--       and km are not on the same grid. Actually, a further
2033!--       interpolation of km onto the u/v-grid is necessary. However, the
2034!--       effect of this error is negligible.
2035          km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) /             &
2036                   phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
2037
2038          surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)    *             &
2039                                        drho_air_zw(k-1)         *             &
2040                                        ( zu(k+1) - zu(k-1)    ) /             &
2041                                        ( km_sfc  + 1.0E-20_wp )
2042          surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)    *             &
2043                                        drho_air_zw(k-1)         *             &
2044                                        ( zu(k+1) - zu(k-1)    ) /             &
2045                                        ( km_sfc  + 1.0E-20_wp )
2046
2047          IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
2048               ABS( u(k+1,j,i) - u(k-1,j,i)   )                                &
2049             )  surf_usm_h%u_0(m) = u(k-1,j,i)
2050
2051          IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) )  >                        &
2052               ABS( v(k+1,j,i) - v(k-1,j,i)   )                                &
2053             )  surf_usm_h%v_0(m) = v(k-1,j,i)
2054
2055       ENDDO
2056
2057    ENDIF
2058
2059 END SUBROUTINE production_e_init
2060
2061
2062!--------------------------------------------------------------------------------------------------!
2063! Description:
2064! ------------
2065!> Execute module-specific actions for all grid points
2066!--------------------------------------------------------------------------------------------------!
2067 SUBROUTINE tcm_actions( location )
2068
2069
2070    CHARACTER (LEN=*) ::  location !<
2071
2072!    INTEGER(iwp) ::  i !<
2073!    INTEGER(iwp) ::  j !<
2074!    INTEGER(iwp) ::  k !<
2075
2076!
2077!-- Here the module-specific actions follow
2078!-- No calls for single grid points are allowed at locations before and
2079!-- after the timestep, since these calls are not within an i,j-loop
2080    SELECT CASE ( location )
2081
2082       CASE ( 'before_timestep' )
2083
2084
2085       CASE ( 'before_prognostic_equations' )
2086
2087          IF ( .NOT. constant_diffusion )  CALL production_e_init
2088
2089
2090       CASE ( 'after_integration' )
2091
2092
2093       CASE ( 'after_timestep' )
2094
2095
2096       CASE ( 'u-tendency' )
2097
2098
2099       CASE ( 'v-tendency' )
2100
2101
2102       CASE ( 'w-tendency' )
2103
2104
2105       CASE ( 'pt-tendency' )
2106
2107
2108       CASE ( 'sa-tendency' )
2109
2110
2111       CASE ( 'e-tendency' )
2112
2113
2114       CASE ( 'q-tendency' )
2115
2116
2117       CASE ( 's-tendency' )
2118
2119
2120       CASE DEFAULT
2121          CONTINUE
2122
2123    END SELECT
2124
2125 END SUBROUTINE tcm_actions
2126
2127
2128!--------------------------------------------------------------------------------------------------!
2129! Description:
2130! ------------
2131!> Execute module-specific actions for grid point i,j
2132!--------------------------------------------------------------------------------------------------!
2133 SUBROUTINE tcm_actions_ij( i, j, location )
2134
2135
2136    CHARACTER (LEN=*) ::  location
2137
2138    INTEGER(iwp) ::  i
2139    INTEGER(iwp) ::  j
2140
2141!
2142!-- Here the module-specific actions follow
2143    SELECT CASE ( location )
2144
2145       CASE ( 'u-tendency' )
2146
2147!--       Next line is to avoid compiler warning about unused variables. Please remove.
2148          IF ( i +  j < 0 )  CONTINUE
2149
2150       CASE ( 'v-tendency' )
2151
2152
2153       CASE ( 'w-tendency' )
2154
2155
2156       CASE ( 'pt-tendency' )
2157
2158
2159       CASE ( 'sa-tendency' )
2160
2161
2162       CASE ( 'e-tendency' )
2163
2164
2165       CASE ( 'q-tendency' )
2166
2167
2168       CASE ( 's-tendency' )
2169
2170
2171       CASE DEFAULT
2172          CONTINUE
2173
2174    END SELECT
2175
2176 END SUBROUTINE tcm_actions_ij
2177
2178
2179!------------------------------------------------------------------------------!
2180! Description:
2181! ------------
2182!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
2183!> Vector-optimized version
2184!------------------------------------------------------------------------------!
2185 SUBROUTINE tcm_prognostic_equations
2186
2187    USE control_parameters,                                                    &
2188        ONLY:  scalar_advec, tsc
2189
2190    IMPLICIT NONE
2191
2192    INTEGER(iwp) ::  i       !< loop index
2193    INTEGER(iwp) ::  j       !< loop index
2194    INTEGER(iwp) ::  k       !< loop index
2195
2196    REAL(wp)     ::  sbt     !< wheighting factor for sub-time step
2197
2198!
2199!-- If required, compute prognostic equation for turbulent kinetic
2200!-- energy (TKE)
2201    IF ( .NOT. constant_diffusion )  THEN
2202
2203       CALL cpu_log( log_point_s(67), 'tke-equation', 'start' )
2204
2205       sbt = tsc(2)
2206       IF ( .NOT. use_upstream_for_tke )  THEN
2207          IF ( scalar_advec == 'bc-scheme' )  THEN
2208
2209             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2210!
2211!--             Bott-Chlond scheme always uses Euler time step. Thus:
2212                sbt = 1.0_wp
2213             ENDIF
2214             tend = 0.0_wp
2215             CALL advec_s_bc( e, 'e' )
2216
2217          ENDIF
2218       ENDIF
2219
2220!
2221!--    TKE-tendency terms with no communication
2222       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
2223          IF ( use_upstream_for_tke )  THEN
2224             tend = 0.0_wp
2225             CALL advec_s_up( e )
2226          ELSE
2227             !$ACC KERNELS PRESENT(tend)
2228             tend = 0.0_wp
2229             !$ACC END KERNELS
2230             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2231                IF ( ws_scheme_sca )  THEN
2232                   CALL advec_s_ws( advc_flags_s, e, 'e',                      &
2233                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
2234                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
2235                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
2236                                    bc_dirichlet_s  .OR.  bc_radiation_s )
2237                ELSE
2238                   CALL advec_s_pw( e )
2239                ENDIF
2240             ELSE
2241                CALL advec_s_up( e )
2242             ENDIF
2243          ENDIF
2244       ENDIF
2245
2246       CALL production_e( .FALSE. )
2247
2248       IF ( .NOT. humidity )  THEN
2249          IF ( ocean_mode )  THEN
2250             CALL diffusion_e( prho, prho_reference )
2251          ELSE
2252             CALL diffusion_e( pt, pt_reference )
2253          ENDIF
2254       ELSE
2255          CALL diffusion_e( vpt, pt_reference )
2256       ENDIF
2257
2258!
2259!--    Additional sink term for flows through plant canopies
2260       IF ( plant_canopy )  CALL pcm_tendency( 6 )
2261
2262!       CALL user_actions( 'e-tendency' ) ToDo: find general solution for circular dependency between modules
2263
2264!
2265!--    Prognostic equation for TKE.
2266!--    Eliminate negative TKE values, which can occur due to numerical
2267!--    reasons in the course of the integration. In such cases the old TKE
2268!--    value is reduced by 90%.
2269       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2270       !$ACC PRESENT(e, tend, te_m, wall_flags_total_0) &
2271       !$ACC PRESENT(tsc(3:3)) &
2272       !$ACC PRESENT(e_p)
2273       DO  i = nxl, nxr
2274          DO  j = nys, nyn
2275             !following directive is required to vectorize on Intel19
2276             !DIR$ IVDEP
2277             DO  k = nzb+1, nzt
2278                e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2279                                                 tsc(3) * te_m(k,j,i) )        &
2280                                        )                                      &
2281                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2282                                       BTEST( wall_flags_total_0(k,j,i), 0 )   &
2283                                          )
2284                IF ( e_p(k,j,i) < 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
2285             ENDDO
2286          ENDDO
2287       ENDDO
2288
2289!
2290!--    Calculate tendencies for the next Runge-Kutta step
2291       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2292          IF ( intermediate_timestep_count == 1 )  THEN
2293             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2294             !$ACC PRESENT(tend, te_m)
2295             DO  i = nxl, nxr
2296                DO  j = nys, nyn
2297                   DO  k = nzb+1, nzt
2298                      te_m(k,j,i) = tend(k,j,i)
2299                   ENDDO
2300                ENDDO
2301             ENDDO
2302          ELSEIF ( intermediate_timestep_count < &
2303                   intermediate_timestep_count_max )  THEN
2304             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2305             !$ACC PRESENT(tend, te_m)
2306             DO  i = nxl, nxr
2307                DO  j = nys, nyn
2308                   DO  k = nzb+1, nzt
2309                      te_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2310                                     + 5.3125_wp * te_m(k,j,i)
2311                   ENDDO
2312                ENDDO
2313             ENDDO
2314          ENDIF
2315       ENDIF
2316
2317       CALL cpu_log( log_point_s(67), 'tke-equation', 'stop' )
2318
2319    ENDIF   ! TKE equation
2320
2321!
2322!-- If required, compute prognostic equation for TKE dissipation rate
2323    IF ( rans_tke_e )  THEN
2324
2325       CALL cpu_log( log_point_s(64), 'diss-equation', 'start' )
2326
2327       sbt = tsc(2)
2328       IF ( .NOT. use_upstream_for_tke )  THEN
2329          IF ( scalar_advec == 'bc-scheme' )  THEN
2330
2331             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2332!
2333!--             Bott-Chlond scheme always uses Euler time step. Thus:
2334                sbt = 1.0_wp
2335             ENDIF
2336             tend = 0.0_wp
2337             CALL advec_s_bc( diss, 'diss' )
2338
2339          ENDIF
2340       ENDIF
2341
2342!
2343!--    dissipation-tendency terms with no communication
2344       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
2345          IF ( use_upstream_for_tke )  THEN
2346             tend = 0.0_wp
2347             CALL advec_s_up( diss )
2348          ELSE
2349             tend = 0.0_wp
2350             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2351                IF ( ws_scheme_sca )  THEN
2352                   CALL advec_s_ws( advc_flags_s, diss, 'diss',                &
2353                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
2354                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
2355                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
2356                                    bc_dirichlet_s  .OR.  bc_radiation_s )
2357                ELSE
2358                   CALL advec_s_pw( diss )
2359                ENDIF
2360             ELSE
2361                CALL advec_s_up( diss )
2362             ENDIF
2363          ENDIF
2364       ENDIF
2365!
2366!--    Production of TKE dissipation rate
2367       CALL production_e( .TRUE. )
2368!
2369!--    Diffusion term of TKE dissipation rate
2370       CALL diffusion_diss
2371!
2372!--    Additional sink term for flows through plant canopies
2373!        IF ( plant_canopy )  CALL pcm_tendency( ? )         !> @todo not yet implemented
2374
2375!       CALL user_actions( 'e-tendency' ) ToDo: find general solution for circular dependency between modules
2376
2377!
2378!--    Prognostic equation for TKE dissipation.
2379!--    Eliminate negative dissipation values, which can occur due to numerical
2380!--    reasons in the course of the integration. In such cases the old
2381!--    dissipation value is reduced by 90%.
2382       DO  i = nxl, nxr
2383          DO  j = nys, nyn
2384             DO  k = nzb+1, nzt
2385                diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +  &
2386                                                 tsc(3) * tdiss_m(k,j,i) )     &
2387                                        )                                      &
2388                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2389                                       BTEST( wall_flags_total_0(k,j,i), 0 )   &
2390                                          )
2391                IF ( diss_p(k,j,i) < 0.0_wp )                                  &
2392                   diss_p(k,j,i) = 0.1_wp * diss(k,j,i)
2393             ENDDO
2394          ENDDO
2395       ENDDO
2396
2397!
2398!--    Calculate tendencies for the next Runge-Kutta step
2399       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2400          IF ( intermediate_timestep_count == 1 )  THEN
2401             DO  i = nxl, nxr
2402                DO  j = nys, nyn
2403                   DO  k = nzb+1, nzt
2404                      tdiss_m(k,j,i) = tend(k,j,i)
2405                   ENDDO
2406                ENDDO
2407             ENDDO
2408          ELSEIF ( intermediate_timestep_count < &
2409                   intermediate_timestep_count_max )  THEN
2410             DO  i = nxl, nxr
2411                DO  j = nys, nyn
2412                   DO  k = nzb+1, nzt
2413                      tdiss_m(k,j,i) =   -9.5625_wp * tend(k,j,i)              &
2414                                        + 5.3125_wp * tdiss_m(k,j,i)
2415                   ENDDO
2416                ENDDO
2417             ENDDO
2418          ENDIF
2419       ENDIF
2420
2421       CALL cpu_log( log_point_s(64), 'diss-equation', 'stop' )
2422
2423    ENDIF
2424
2425 END SUBROUTINE tcm_prognostic_equations
2426
2427
2428!------------------------------------------------------------------------------!
2429! Description:
2430! ------------
2431!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
2432!> Cache-optimized version
2433!------------------------------------------------------------------------------!
2434 SUBROUTINE tcm_prognostic_equations_ij( i, j, i_omp, tn )
2435
2436    USE arrays_3d,                                                             &
2437        ONLY:  diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, flux_l_diss,      &
2438               flux_l_e, flux_s_diss, flux_s_e
2439
2440    USE control_parameters,                                                    &
2441        ONLY:  tsc
2442
2443    IMPLICIT NONE
2444
2445    INTEGER(iwp) ::  i       !< loop index x direction
2446    INTEGER(iwp) ::  i_omp   !< first loop index of i-loop in prognostic_equations
2447    INTEGER(iwp) ::  j       !< loop index y direction
2448    INTEGER(iwp) ::  k       !< loop index z direction
2449    INTEGER(iwp) ::  tn      !< task number of openmp task
2450
2451!
2452!-- If required, compute prognostic equation for turbulent kinetic
2453!-- energy (TKE)
2454    IF ( .NOT. constant_diffusion )  THEN
2455
2456!
2457!--    Tendency-terms for TKE
2458       tend(:,j,i) = 0.0_wp
2459       IF ( timestep_scheme(1:5) == 'runge'  &
2460           .AND.  .NOT. use_upstream_for_tke )  THEN
2461           IF ( ws_scheme_sca )  THEN
2462               CALL advec_s_ws( advc_flags_s,                                  &
2463                                i, j, e, 'e', flux_s_e, diss_s_e,              &
2464                                flux_l_e, diss_l_e , i_omp, tn,                &
2465                                bc_dirichlet_l  .OR.  bc_radiation_l,          &
2466                                bc_dirichlet_n  .OR.  bc_radiation_n,          &
2467                                bc_dirichlet_r  .OR.  bc_radiation_r,          &
2468                                bc_dirichlet_s  .OR.  bc_radiation_s )
2469           ELSE
2470               CALL advec_s_pw( i, j, e )
2471           ENDIF
2472       ELSE
2473          CALL advec_s_up( i, j, e )
2474       ENDIF
2475
2476       CALL production_e_ij( i, j, .FALSE. )
2477
2478       IF ( .NOT. humidity )  THEN
2479          IF ( ocean_mode )  THEN
2480             CALL diffusion_e_ij( i, j, prho, prho_reference )
2481          ELSE
2482             CALL diffusion_e_ij( i, j, pt, pt_reference )
2483          ENDIF
2484       ELSE
2485          CALL diffusion_e_ij( i, j, vpt, pt_reference )
2486       ENDIF
2487
2488!
2489!--    Additional sink term for flows through plant canopies
2490       IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 )
2491
2492!       CALL user_actions( i, j, 'e-tendency' ) ToDo: find general solution for circular dependency between modules
2493
2494!
2495!--    Prognostic equation for TKE.
2496!--    Eliminate negative TKE values, which can occur due to numerical
2497!--    reasons in the course of the integration. In such cases the old
2498!--    TKE value is reduced by 90%.
2499       DO  k = nzb+1, nzt
2500          e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +           &
2501                                              tsc(3) * te_m(k,j,i) )           &
2502                                  )                                            &
2503                                 * MERGE( 1.0_wp, 0.0_wp,                      &
2504                                    BTEST( wall_flags_total_0(k,j,i), 0 )      &
2505                                        )
2506          IF ( e_p(k,j,i) <= 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
2507       ENDDO
2508
2509!
2510!--    Calculate tendencies for the next Runge-Kutta step
2511       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2512          IF ( intermediate_timestep_count == 1 )  THEN
2513             DO  k = nzb+1, nzt
2514                te_m(k,j,i) = tend(k,j,i)
2515             ENDDO
2516          ELSEIF ( intermediate_timestep_count < &
2517                   intermediate_timestep_count_max )  THEN
2518             DO  k = nzb+1, nzt
2519                te_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +                     &
2520                                 5.3125_wp * te_m(k,j,i)
2521             ENDDO
2522          ENDIF
2523       ENDIF
2524
2525    ENDIF   ! TKE equation
2526
2527!
2528!-- If required, compute prognostic equation for TKE dissipation rate
2529    IF ( rans_tke_e )  THEN
2530!
2531!--    Tendency-terms for dissipation
2532       tend(:,j,i) = 0.0_wp
2533       IF ( timestep_scheme(1:5) == 'runge'  &
2534           .AND.  .NOT. use_upstream_for_tke )  THEN
2535           IF ( ws_scheme_sca )  THEN
2536               CALL advec_s_ws( advc_flags_s,                                  &
2537                                i, j, diss, 'diss', flux_s_diss, diss_s_diss,  &
2538                                flux_l_diss, diss_l_diss, i_omp, tn,           &
2539                                bc_dirichlet_l  .OR.  bc_radiation_l,          &
2540                                bc_dirichlet_n  .OR.  bc_radiation_n,          &
2541                                bc_dirichlet_r  .OR.  bc_radiation_r,          &
2542                                bc_dirichlet_s  .OR.  bc_radiation_s )
2543           ELSE
2544               CALL advec_s_pw( i, j, diss )
2545           ENDIF
2546       ELSE
2547          CALL advec_s_up( i, j, diss )
2548       ENDIF
2549!
2550!--    Production of TKE dissipation rate
2551       CALL production_e_ij( i, j, .TRUE. )
2552!
2553!--    Diffusion term of TKE dissipation rate
2554       CALL diffusion_diss_ij( i, j )
2555!
2556!--    Additional sink term for flows through plant canopies
2557!        IF ( plant_canopy )  CALL pcm_tendency( i, j, ? )     !> @todo not yet implemented
2558
2559!       CALL user_actions( i, j, 'diss-tendency' ) ToDo: find general solution for circular dependency between modules
2560
2561!
2562!--    Prognostic equation for TKE dissipation
2563!--    Eliminate negative dissipation values, which can occur due to
2564!--    numerical reasons in the course of the integration. In such cases
2565!--    the old dissipation value is reduced by 90%.
2566       DO  k = nzb+1, nzt
2567          diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +     &
2568                                                    tsc(3) * tdiss_m(k,j,i) )  &
2569                                        )                                      &
2570                                        * MERGE( 1.0_wp, 0.0_wp,               &
2571                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
2572                                               )
2573       ENDDO
2574
2575!
2576!--    Calculate tendencies for the next Runge-Kutta step
2577       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2578          IF ( intermediate_timestep_count == 1 )  THEN
2579             DO  k = nzb+1, nzt
2580                tdiss_m(k,j,i) = tend(k,j,i)
2581             ENDDO
2582          ELSEIF ( intermediate_timestep_count < &
2583                   intermediate_timestep_count_max )  THEN
2584             DO  k = nzb+1, nzt
2585                tdiss_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +                  &
2586                                    5.3125_wp * tdiss_m(k,j,i)
2587             ENDDO
2588          ENDIF
2589       ENDIF
2590
2591    ENDIF   ! dissipation equation
2592
2593 END SUBROUTINE tcm_prognostic_equations_ij
2594
2595
2596!------------------------------------------------------------------------------!
2597! Description:
2598! ------------
2599!> Production terms (shear + buoyancy) of the TKE.
2600!> Vector-optimized version
2601!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
2602!>          not considered well!
2603!------------------------------------------------------------------------------!
2604 SUBROUTINE production_e( diss_production )
2605
2606    USE arrays_3d,                                                             &
2607        ONLY:  ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
2608
2609    USE control_parameters,                                                    &
2610        ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
2611               rho_reference, use_single_reference_value, use_surface_fluxes,  &
2612               use_top_fluxes
2613
2614    USE grid_variables,                                                        &
2615        ONLY:  ddx, dx, ddy, dy
2616
2617    USE bulk_cloud_model_mod,                                                  &
2618        ONLY:  bulk_cloud_model
2619
2620    IMPLICIT NONE
2621
2622    LOGICAL :: diss_production
2623
2624    INTEGER(iwp) ::  i       !< running index x-direction
2625    INTEGER(iwp) ::  j       !< running index y-direction
2626    INTEGER(iwp) ::  k       !< running index z-direction
2627    INTEGER(iwp) ::  l       !< running index for different surface type orientation
2628    INTEGER(iwp) ::  m       !< running index surface elements
2629    INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
2630    INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
2631    INTEGER(iwp) ::  flag_nr !< number of masking flag
2632
2633    REAL(wp)     ::  def         !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j
2634    REAL(wp)     ::  flag        !< flag to mask topography
2635    REAL(wp)     ::  k1          !< temporary factor
2636    REAL(wp)     ::  k2          !< temporary factor
2637    REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
2638    REAL(wp)     ::  theta       !< virtual potential temperature
2639    REAL(wp)     ::  temp        !< theta * Exner-function
2640    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
2641    REAL(wp)     ::  usvs        !< momentum flux u"v"
2642    REAL(wp)     ::  vsus        !< momentum flux v"u"
2643    REAL(wp)     ::  wsus        !< momentum flux w"u"
2644    REAL(wp)     ::  wsvs        !< momentum flux w"v"
2645
2646    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx  !< Gradient of u-component in x-direction
2647    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy  !< Gradient of u-component in y-direction
2648    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz  !< Gradient of u-component in z-direction
2649    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx  !< Gradient of v-component in x-direction
2650    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy  !< Gradient of v-component in y-direction
2651    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz  !< Gradient of v-component in z-direction
2652    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx  !< Gradient of w-component in x-direction
2653    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy  !< Gradient of w-component in y-direction
2654    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz  !< Gradient of w-component in z-direction
2655    REAL(wp), DIMENSION(nzb+1:nzt) ::  tmp_flux  !< temporary flux-array in z-direction
2656
2657
2658
2659!
2660!-- Calculate TKE production by shear. Calculate gradients at all grid
2661!-- points first, gradients at surface-bounded grid points will be
2662!-- overwritten further below.
2663    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) &
2664    !$ACC PRIVATE(surf_s, surf_e) &
2665    !$ACC PRIVATE(dudx(:), dudy(:), dudz(:), dvdx(:), dvdy(:), dvdz(:), dwdx(:), dwdy(:), dwdz(:)) &
2666    !$ACC PRESENT(e, u, v, w, diss, dd2zu, ddzw, km, wall_flags_total_0) &
2667    !$ACC PRESENT(tend) &
2668    !$ACC PRESENT(surf_def_h(0:1), surf_def_v(0:3)) &
2669    !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
2670    !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3))
2671    DO  i = nxl, nxr
2672       DO  j = nys, nyn
2673          !$ACC LOOP PRIVATE(k)
2674          DO  k = nzb+1, nzt
2675
2676             dudx(k) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
2677             dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                 &
2678                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
2679             dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                 &
2680                                   u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
2681
2682             dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                 &
2683                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
2684             dvdy(k) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
2685             dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                 &
2686                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
2687
2688             dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                 &
2689                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
2690             dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                 &
2691                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
2692             dwdz(k) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
2693
2694          ENDDO
2695
2696
2697          flag_nr = 29
2698
2699
2700          IF ( constant_flux_layer )  THEN
2701!
2702
2703             flag_nr = 0
2704
2705!--          Position beneath wall
2706!--          (2) - Will allways be executed.
2707!--          'bottom and wall: use u_0,v_0 and wall functions'
2708!
2709!--          Compute gradients at north- and south-facing surfaces.
2710!--          First, for default surfaces, then for urban surfaces.
2711!--          Note, so far no natural vertical surfaces implemented
2712             DO  l = 0, 1
2713                surf_s = surf_def_v(l)%start_index(j,i)
2714                surf_e = surf_def_v(l)%end_index(j,i)
2715                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2716                DO  m = surf_s, surf_e
2717                   k           = surf_def_v(l)%k(m)
2718                   usvs        = surf_def_v(l)%mom_flux_tke(0,m)
2719                   wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
2720
2721                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
2722                                   * 0.5_wp * dy
2723!
2724!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2725                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2726                              BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
2727                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2728                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2729                ENDDO
2730!
2731!--             Natural surfaces
2732                surf_s = surf_lsm_v(l)%start_index(j,i)
2733                surf_e = surf_lsm_v(l)%end_index(j,i)
2734                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2735                DO  m = surf_s, surf_e
2736                   k           = surf_lsm_v(l)%k(m)
2737                   usvs        = surf_lsm_v(l)%mom_flux_tke(0,m)
2738                   wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
2739
2740                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
2741                                   * 0.5_wp * dy
2742!
2743!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2744                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2745                              BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
2746                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2747                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2748                ENDDO
2749!
2750!--             Urban surfaces
2751                surf_s = surf_usm_v(l)%start_index(j,i)
2752                surf_e = surf_usm_v(l)%end_index(j,i)
2753                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2754                DO  m = surf_s, surf_e
2755                   k           = surf_usm_v(l)%k(m)
2756                   usvs        = surf_usm_v(l)%mom_flux_tke(0,m)
2757                   wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
2758
2759                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
2760                                   * 0.5_wp * dy
2761!
2762!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2763                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2764                              BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
2765                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2766                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2767                ENDDO
2768             ENDDO
2769!
2770!--          Compute gradients at east- and west-facing walls
2771             DO  l = 2, 3
2772                surf_s = surf_def_v(l)%start_index(j,i)
2773                surf_e = surf_def_v(l)%end_index(j,i)
2774                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2775                DO  m = surf_s, surf_e
2776                   k     = surf_def_v(l)%k(m)
2777                   vsus  = surf_def_v(l)%mom_flux_tke(0,m)
2778                   wsus  = surf_def_v(l)%mom_flux_tke(1,m)
2779
2780                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
2781                                      * 0.5_wp * dx
2782!
2783!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2784                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2785                              BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
2786                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2787                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2788                ENDDO
2789!
2790!--             Natural surfaces
2791                surf_s = surf_lsm_v(l)%start_index(j,i)
2792                surf_e = surf_lsm_v(l)%end_index(j,i)
2793                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2794                DO  m = surf_s, surf_e
2795                   k     = surf_lsm_v(l)%k(m)
2796                   vsus  = surf_lsm_v(l)%mom_flux_tke(0,m)
2797                   wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
2798
2799                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
2800                                      * 0.5_wp * dx
2801!
2802!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2803                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2804                              BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
2805                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2806                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2807                ENDDO
2808!
2809!--             Urban surfaces
2810                surf_s = surf_usm_v(l)%start_index(j,i)
2811                surf_e = surf_usm_v(l)%end_index(j,i)
2812                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2813                DO  m = surf_s, surf_e
2814                   k     = surf_usm_v(l)%k(m)
2815                   vsus  = surf_usm_v(l)%mom_flux_tke(0,m)
2816                   wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
2817
2818                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
2819                                      * 0.5_wp * dx
2820!
2821!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2822                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2823                              BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
2824                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2825                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2826                ENDDO
2827             ENDDO
2828!
2829!--          Compute gradients at upward-facing surfaces
2830             surf_s = surf_def_h(0)%start_index(j,i)
2831             surf_e = surf_def_h(0)%end_index(j,i)
2832             !$ACC LOOP PRIVATE(m, k)
2833             DO  m = surf_s, surf_e
2834                k = surf_def_h(0)%k(m)
2835!
2836!--             Please note, actually, an interpolation of u_0 and v_0
2837!--             onto the grid center would be required. However, this
2838!--             would require several data transfers between 2D-grid and
2839!--             wall type. The effect of this missing interpolation is
2840!--             negligible. (See also production_e_init).
2841                dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
2842                dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
2843
2844             ENDDO
2845!
2846!--          Natural surfaces
2847             surf_s = surf_lsm_h%start_index(j,i)
2848             surf_e = surf_lsm_h%end_index(j,i)
2849             !$ACC LOOP PRIVATE(m, k)
2850             DO  m = surf_s, surf_e
2851                k = surf_lsm_h%k(m)
2852
2853                dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
2854                dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
2855
2856             ENDDO
2857!
2858!--          Urban surfaces
2859             surf_s = surf_usm_h%start_index(j,i)
2860             surf_e = surf_usm_h%end_index(j,i)
2861             !$ACC LOOP PRIVATE(m, k)
2862             DO  m = surf_s, surf_e
2863                k = surf_usm_h%k(m)
2864
2865                dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
2866                dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
2867
2868             ENDDO
2869!
2870!--          Compute gradients at downward-facing walls, only for
2871!--          non-natural default surfaces
2872             surf_s = surf_def_h(1)%start_index(j,i)
2873             surf_e = surf_def_h(1)%end_index(j,i)
2874             !$ACC LOOP PRIVATE(m, k)
2875             DO  m = surf_s, surf_e
2876                k = surf_def_h(1)%k(m)
2877
2878                dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
2879                dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
2880
2881             ENDDO
2882
2883
2884          ENDIF
2885
2886
2887          !$ACC LOOP PRIVATE(k, def, flag)
2888          DO  k = nzb+1, nzt
2889
2890             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
2891                              dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 +           &
2892                              dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +           &
2893                   2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) +              &
2894                              dwdy(k)*dvdz(k) )
2895
2896             IF ( def < 0.0_wp )  def = 0.0_wp
2897
2898             flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),flag_nr) )
2899
2900             IF ( .NOT. diss_production )  THEN
2901
2902!--             Compute tendency for TKE-production from shear
2903                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
2904
2905             ELSE
2906
2907!--             RANS mode: Compute tendency for dissipation-rate-production from shear
2908                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
2909                              diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
2910
2911             ENDIF
2912
2913          ENDDO
2914
2915
2916       ENDDO
2917    ENDDO
2918
2919!
2920!-- If required, calculate TKE production by buoyancy
2921    IF ( .NOT. neutral )  THEN
2922
2923       IF ( .NOT. humidity )  THEN
2924
2925          IF ( ocean_mode )  THEN
2926!
2927!--          So far in the ocean no special treatment of density flux
2928!--          in the bottom and top surface layer
2929             DO  i = nxl, nxr
2930                DO  j = nys, nyn
2931
2932                   DO  k = nzb+1, nzt
2933                      tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
2934                   ENDDO
2935!
2936!--                Treatment of near-surface grid points, at up- and down-
2937!--                ward facing surfaces
2938                   IF ( use_surface_fluxes )  THEN
2939                      DO  l = 0, 1
2940                         surf_s = surf_def_h(l)%start_index(j,i)
2941                         surf_e = surf_def_h(l)%end_index(j,i)
2942                         DO  m = surf_s, surf_e
2943                            k = surf_def_h(l)%k(m)
2944                            tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
2945                         ENDDO
2946                      ENDDO
2947                   ENDIF
2948
2949                   IF ( use_top_fluxes )  THEN
2950                      surf_s = surf_def_h(2)%start_index(j,i)
2951                      surf_e = surf_def_h(2)%end_index(j,i)
2952                      DO  m = surf_s, surf_e
2953                         k = surf_def_h(2)%k(m)
2954                         tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
2955                      ENDDO
2956                   ENDIF
2957
2958                   IF ( .NOT. diss_production )  THEN
2959
2960!--                   Compute tendency for TKE-production from shear
2961                      DO  k = nzb+1, nzt
2962                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
2963                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2964                                       MERGE( rho_reference, prho(k,j,i),       &
2965                                              use_single_reference_value ) )
2966                      ENDDO
2967
2968                   ELSE
2969
2970!--                   RANS mode: Compute tendency for dissipation-rate-production from shear
2971                      DO  k = nzb+1, nzt
2972                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
2973                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2974                                       MERGE( rho_reference, prho(k,j,i),       &
2975                                              use_single_reference_value ) ) *  &
2976                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
2977                                       c_3
2978                      ENDDO
2979
2980                   ENDIF
2981
2982                ENDDO
2983             ENDDO
2984
2985          ELSE ! or IF ( .NOT. ocean_mode )  THEN
2986
2987             !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
2988             !$ACC PRIVATE(surf_s, surf_e) &
2989             !$ACC PRIVATE(tmp_flux(nzb+1:nzt)) &
2990             !$ACC PRESENT(e, diss, kh, pt, dd2zu, drho_air_zw, wall_flags_total_0) &
2991             !$ACC PRESENT(tend) &
2992             !$ACC PRESENT(surf_def_h(0:2)) &
2993             !$ACC PRESENT(surf_lsm_h) &
2994             !$ACC PRESENT(surf_usm_h)
2995             DO  i = nxl, nxr
2996                DO  j = nys, nyn
2997
2998                   !$ACC LOOP PRIVATE(k)
2999                   DO  k = nzb+1, nzt
3000                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
3001                   ENDDO
3002
3003                   IF ( use_surface_fluxes )  THEN
3004!
3005!--                   Default surfaces, up- and downward-facing
3006                      DO  l = 0, 1
3007                         surf_s = surf_def_h(l)%start_index(j,i)
3008                         surf_e = surf_def_h(l)%end_index(j,i)
3009                         !$ACC LOOP PRIVATE(m, k)
3010                         DO  m = surf_s, surf_e
3011                            k = surf_def_h(l)%k(m)
3012                            tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3013                         ENDDO
3014                      ENDDO
3015!
3016!--                   Natural surfaces
3017                      surf_s = surf_lsm_h%start_index(j,i)
3018                      surf_e = surf_lsm_h%end_index(j,i)
3019                      !$ACC LOOP PRIVATE(m, k)
3020                      DO  m = surf_s, surf_e
3021                         k = surf_lsm_h%k(m)
3022                         tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
3023                      ENDDO
3024!
3025!--                   Urban surfaces
3026                      surf_s = surf_usm_h%start_index(j,i)
3027                      surf_e = surf_usm_h%end_index(j,i)
3028                      !$ACC LOOP PRIVATE(m, k)
3029                      DO  m = surf_s, surf_e
3030                         k = surf_usm_h%k(m)
3031                         tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
3032                      ENDDO
3033                   ENDIF
3034
3035                   IF ( use_top_fluxes )  THEN
3036                      surf_s = surf_def_h(2)%start_index(j,i)
3037                      surf_e = surf_def_h(2)%end_index(j,i)
3038                      !$ACC LOOP PRIVATE(m, k)
3039                      DO  m = surf_s, surf_e
3040                         k = surf_def_h(2)%k(m)
3041                         tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3042                      ENDDO
3043                   ENDIF
3044
3045                   IF ( .NOT. diss_production )  THEN
3046
3047!--                   Compute tendency for TKE-production from shear
3048                     !$ACC LOOP PRIVATE(k, flag)
3049                      DO  k = nzb+1, nzt
3050                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3051                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3052                                       MERGE( pt_reference, pt(k,j,i),          &
3053                                              use_single_reference_value ) )
3054                      ENDDO
3055
3056                   ELSE
3057
3058!--                   RANS mode: Compute tendency for dissipation-rate-production from shear
3059                      DO  k = nzb+1, nzt
3060                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3061                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3062                                       MERGE( pt_reference, pt(k,j,i),          &
3063                                              use_single_reference_value ) ) *  &
3064                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
3065                                       c_3
3066                      ENDDO
3067
3068                   ENDIF
3069
3070                ENDDO
3071             ENDDO
3072
3073          ENDIF ! from IF ( .NOT. ocean_mode )
3074
3075       ELSE ! or IF ( humidity )  THEN
3076
3077          DO  i = nxl, nxr
3078             DO  j = nys, nyn
3079
3080                DO  k = nzb+1, nzt
3081
3082                   IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3083                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3084                      k2 = 0.61_wp * pt(k,j,i)
3085                      tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
3086                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3087                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
3088                                      ) * dd2zu(k)
3089                   ELSE IF ( bulk_cloud_model )  THEN
3090                      IF ( ql(k,j,i) == 0.0_wp )  THEN
3091                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3092                         k2 = 0.61_wp * pt(k,j,i)
3093                      ELSE
3094                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3095                         temp  = theta * exner(k)
3096                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
3097                                       ( q(k,j,i) - ql(k,j,i) ) *              &
3098                              ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
3099                              ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
3100                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3101                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3102                      ENDIF
3103                      tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
3104                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3105                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
3106                                      ) * dd2zu(k)
3107                   ELSE IF ( cloud_droplets )  THEN
3108                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3109                      k2 = 0.61_wp * pt(k,j,i)
3110                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3111                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3112                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
3113                                        pt(k,j,i) * ( ql(k+1,j,i) -            &
3114                                        ql(k-1,j,i) ) ) * dd2zu(k)
3115                   ENDIF
3116
3117                ENDDO
3118
3119                IF ( use_surface_fluxes )  THEN
3120!
3121!--                Treat horizontal default surfaces
3122                   DO  l = 0, 1
3123                      surf_s = surf_def_h(l)%start_index(j,i)
3124                      surf_e = surf_def_h(l)%end_index(j,i)
3125                      DO  m = surf_s, surf_e
3126                         k = surf_def_h(l)%k(m)
3127
3128                         IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3129                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3130                            k2 = 0.61_wp * pt(k,j,i)
3131                         ELSE IF ( bulk_cloud_model )  THEN
3132                            IF ( ql(k,j,i) == 0.0_wp )  THEN
3133                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3134                               k2 = 0.61_wp * pt(k,j,i)
3135                            ELSE
3136                               theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3137                               temp  = theta * exner(k)
3138                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
3139                                          ( q(k,j,i) - ql(k,j,i) ) *           &
3140                                 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3141                                 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3142                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3143                               k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3144                            ENDIF
3145                         ELSE IF ( cloud_droplets )  THEN
3146                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3147                            k2 = 0.61_wp * pt(k,j,i)
3148                         ENDIF
3149
3150                         tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
3151                                         k2 * surf_def_h(l)%qsws(m)            &
3152                                       ) * drho_air_zw(k-1)
3153                      ENDDO
3154                   ENDDO
3155!
3156!--                Treat horizontal natural surfaces
3157                   surf_s = surf_lsm_h%start_index(j,i)
3158                   surf_e = surf_lsm_h%end_index(j,i)
3159                   DO  m = surf_s, surf_e
3160                      k = surf_lsm_h%k(m)
3161
3162                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3163                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3164                         k2 = 0.61_wp * pt(k,j,i)
3165                      ELSE IF ( bulk_cloud_model )  THEN
3166                         IF ( ql(k,j,i) == 0.0_wp )  THEN
3167                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3168                            k2 = 0.61_wp * pt(k,j,i)
3169                         ELSE
3170                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3171                            temp  = theta * exner(k)
3172                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3173                                          ( q(k,j,i) - ql(k,j,i) ) *           &
3174                                 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3175                                 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3176                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3177                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3178                         ENDIF
3179                      ELSE IF ( cloud_droplets )  THEN
3180                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3181                         k2 = 0.61_wp * pt(k,j,i)
3182                      ENDIF
3183
3184                      tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
3185                                      k2 * surf_lsm_h%qsws(m)                  &
3186                                    ) * drho_air_zw(k-1)
3187                   ENDDO
3188!
3189!--                Treat horizontal urban surfaces
3190                   surf_s = surf_usm_h%start_index(j,i)
3191                   surf_e = surf_usm_h%end_index(j,i)
3192                   DO  m = surf_s, surf_e
3193                      k = surf_usm_h%k(m)
3194
3195                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3196                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3197                         k2 = 0.61_wp * pt(k,j,i)
3198                      ELSE IF ( bulk_cloud_model )  THEN
3199                         IF ( ql(k,j,i) == 0.0_wp )  THEN
3200                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3201                            k2 = 0.61_wp * pt(k,j,i)
3202                         ELSE
3203                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3204                            temp  = theta * exner(k)
3205                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3206                                          ( q(k,j,i) - ql(k,j,i) ) *           &
3207                                 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3208                                 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3209                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3210                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3211                         ENDIF
3212                      ELSE IF ( cloud_droplets )  THEN
3213                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3214                         k2 = 0.61_wp * pt(k,j,i)
3215                      ENDIF
3216
3217                      tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
3218                                      k2 * surf_usm_h%qsws(m)                  &
3219                                    ) * drho_air_zw(k-1)
3220                   ENDDO
3221
3222                ENDIF ! from IF ( use_surface_fluxes )  THEN
3223
3224                IF ( use_top_fluxes )  THEN
3225
3226                   surf_s = surf_def_h(2)%start_index(j,i)
3227                   surf_e = surf_def_h(2)%end_index(j,i)
3228                   DO  m = surf_s, surf_e
3229                      k = surf_def_h(2)%k(m)
3230
3231                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3232                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3233                         k2 = 0.61_wp * pt(k,j,i)
3234                      ELSE IF ( bulk_cloud_model )  THEN
3235                         IF ( ql(k,j,i) == 0.0_wp )  THEN
3236                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3237                            k2 = 0.61_wp * pt(k,j,i)
3238                         ELSE
3239                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3240                            temp  = theta * exner(k)
3241                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3242                                       ( q(k,j,i) - ql(k,j,i) ) *              &
3243                              ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
3244                              ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
3245                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3246                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3247                         ENDIF
3248                      ELSE IF ( cloud_droplets )  THEN
3249                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3250                         k2 = 0.61_wp * pt(k,j,i)
3251                      ENDIF
3252
3253                      tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
3254                                      k2 * surf_def_h(2)%qsws(m)               &
3255                                    ) * drho_air_zw(k)
3256
3257                   ENDDO
3258
3259                ENDIF ! from IF ( use_top_fluxes )  THEN
3260
3261                IF ( .NOT. diss_production )  THEN
3262
3263!--                Compute tendency for TKE-production from shear
3264                   DO  k = nzb+1, nzt
3265                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3266                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3267                                    MERGE( vpt_reference, vpt(k,j,i),          &
3268                                           use_single_reference_value ) )
3269                   ENDDO
3270
3271                ELSE
3272
3273!--                RANS mode: Compute tendency for dissipation-rate-production from shear
3274                   DO  k = nzb+1, nzt
3275                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3276                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
3277                                    MERGE( vpt_reference, vpt(k,j,i),          &
3278                                           use_single_reference_value ) ) *    &
3279                                    diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
3280                                    c_3
3281                   ENDDO
3282
3283                ENDIF
3284
3285             ENDDO
3286          ENDDO
3287
3288       ENDIF
3289
3290    ENDIF
3291
3292 END SUBROUTINE production_e
3293
3294
3295!------------------------------------------------------------------------------!
3296! Description:
3297! ------------
3298!> Production terms (shear + buoyancy) of the TKE.
3299!> Cache-optimized version
3300!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
3301!>          not considered well!
3302!------------------------------------------------------------------------------!
3303 SUBROUTINE production_e_ij( i, j, diss_production )
3304
3305    USE arrays_3d,                                                             &
3306        ONLY:  ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
3307
3308    USE control_parameters,                                                    &
3309        ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
3310               rho_reference, use_single_reference_value, use_surface_fluxes,  &
3311               use_top_fluxes
3312
3313    USE grid_variables,                                                        &
3314        ONLY:  ddx, dx, ddy, dy
3315
3316    USE bulk_cloud_model_mod,                                                  &
3317        ONLY:  bulk_cloud_model
3318
3319    IMPLICIT NONE
3320
3321    LOGICAL :: diss_production
3322
3323    INTEGER(iwp) ::  i       !< running index x-direction
3324    INTEGER(iwp) ::  j       !< running index y-direction
3325    INTEGER(iwp) ::  k       !< running index z-direction
3326    INTEGER(iwp) ::  l       !< running index for different surface type orientation
3327    INTEGER(iwp) ::  m       !< running index surface elements
3328    INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
3329    INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
3330    INTEGER(iwp) ::  flag_nr !< number of masking flag
3331
3332    REAL(wp)     ::  def         !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j
3333    REAL(wp)     ::  flag        !< flag to mask topography
3334    REAL(wp)     ::  k1          !< temporary factor
3335    REAL(wp)     ::  k2          !< temporary factor
3336    REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
3337    REAL(wp)     ::  theta       !< virtual potential temperature
3338    REAL(wp)     ::  temp        !< theta * Exner-function
3339    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
3340    REAL(wp)     ::  usvs        !< momentum flux u"v"
3341    REAL(wp)     ::  vsus        !< momentum flux v"u"
3342    REAL(wp)     ::  wsus        !< momentum flux w"u"
3343    REAL(wp)     ::  wsvs        !< momentum flux w"v"
3344
3345    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx  !< Gradient of u-component in x-direction
3346    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy  !< Gradient of u-component in y-direction
3347    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz  !< Gradient of u-component in z-direction
3348    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx  !< Gradient of v-component in x-direction
3349    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy  !< Gradient of v-component in y-direction
3350    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz  !< Gradient of v-component in z-direction
3351    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx  !< Gradient of w-component in x-direction
3352    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy  !< Gradient of w-component in y-direction
3353    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz  !< Gradient of w-component in z-direction
3354    REAL(wp), DIMENSION(nzb+1:nzt) ::  tmp_flux  !< temporary flux-array in z-direction
3355
3356
3357
3358!
3359!-- Calculate TKE production by shear. Calculate gradients at all grid
3360!-- points first, gradients at surface-bounded grid points will be
3361!-- overwritten further below.
3362    DO  k = nzb+1, nzt
3363
3364       dudx(k) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
3365       dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                 &
3366                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
3367       dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                 &
3368                             u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
3369
3370       dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                 &
3371                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
3372       dvdy(k) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
3373       dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                 &
3374                             v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
3375
3376       dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                 &
3377                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
3378       dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                 &
3379                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
3380       dwdz(k) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
3381
3382    ENDDO
3383
3384    flag_nr = 29
3385
3386    IF ( constant_flux_layer )  THEN
3387
3388       flag_nr = 0
3389
3390!--    Position beneath wall
3391!--    (2) - Will allways be executed.
3392!--    'bottom and wall: use u_0,v_0 and wall functions'
3393!
3394!--    Compute gradients at north- and south-facing surfaces.
3395!--    First, for default surfaces, then for urban surfaces.
3396!--    Note, so far no natural vertical surfaces implemented
3397       DO  l = 0, 1
3398          surf_s = surf_def_v(l)%start_index(j,i)
3399          surf_e = surf_def_v(l)%end_index(j,i)
3400          DO  m = surf_s, surf_e
3401             k           = surf_def_v(l)%k(m)
3402             usvs        = surf_def_v(l)%mom_flux_tke(0,m)
3403             wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
3404
3405             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
3406                             * 0.5_wp * dy
3407!
3408!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3409             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3410                        BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
3411             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3412             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3413          ENDDO
3414!
3415!--       Natural surfaces
3416          surf_s = surf_lsm_v(l)%start_index(j,i)
3417          surf_e = surf_lsm_v(l)%end_index(j,i)
3418          DO  m = surf_s, surf_e
3419             k           = surf_lsm_v(l)%k(m)
3420             usvs        = surf_lsm_v(l)%mom_flux_tke(0,m)
3421             wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
3422
3423             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
3424                             * 0.5_wp * dy
3425!
3426!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3427             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3428                        BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
3429             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3430             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3431          ENDDO
3432!
3433!--       Urban surfaces
3434          surf_s = surf_usm_v(l)%start_index(j,i)
3435          surf_e = surf_usm_v(l)%end_index(j,i)
3436          DO  m = surf_s, surf_e
3437             k           = surf_usm_v(l)%k(m)
3438             usvs        = surf_usm_v(l)%mom_flux_tke(0,m)
3439             wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
3440
3441             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
3442                             * 0.5_wp * dy
3443!
3444!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3445             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3446                        BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
3447             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3448             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3449          ENDDO
3450       ENDDO
3451!
3452!--    Compute gradients at east- and west-facing walls
3453       DO  l = 2, 3
3454          surf_s = surf_def_v(l)%start_index(j,i)
3455          surf_e = surf_def_v(l)%end_index(j,i)
3456          DO  m = surf_s, surf_e
3457             k     = surf_def_v(l)%k(m)
3458             vsus  = surf_def_v(l)%mom_flux_tke(0,m)
3459             wsus  = surf_def_v(l)%mom_flux_tke(1,m)
3460
3461             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
3462                                * 0.5_wp * dx
3463!
3464!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3465             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3466                        BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
3467             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3468             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3469          ENDDO
3470!
3471!--       Natural surfaces
3472          surf_s = surf_lsm_v(l)%start_index(j,i)
3473          surf_e = surf_lsm_v(l)%end_index(j,i)
3474          DO  m = surf_s, surf_e
3475             k     = surf_lsm_v(l)%k(m)
3476             vsus  = surf_lsm_v(l)%mom_flux_tke(0,m)
3477             wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
3478
3479             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
3480                                * 0.5_wp * dx
3481!
3482!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3483             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3484                        BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
3485             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3486             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3487          ENDDO
3488!
3489!--       Urban surfaces
3490          surf_s = surf_usm_v(l)%start_index(j,i)
3491          surf_e = surf_usm_v(l)%end_index(j,i)
3492          DO  m = surf_s, surf_e
3493             k     = surf_usm_v(l)%k(m)
3494             vsus  = surf_usm_v(l)%mom_flux_tke(0,m)
3495             wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
3496
3497             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
3498                                * 0.5_wp * dx
3499!
3500!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3501             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3502                        BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
3503             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3504             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3505          ENDDO
3506       ENDDO
3507!
3508!--    Compute gradients at upward-facing surfaces
3509       surf_s = surf_def_h(0)%start_index(j,i)
3510       surf_e = surf_def_h(0)%end_index(j,i)
3511       DO  m = surf_s, surf_e
3512          k = surf_def_h(0)%k(m)
3513!
3514!--       Please note, actually, an interpolation of u_0 and v_0
3515!--       onto the grid center would be required. However, this
3516!--       would require several data transfers between 2D-grid and
3517!--       wall type. The effect of this missing interpolation is
3518!--       negligible. (See also production_e_init).
3519          dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
3520          dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
3521
3522       ENDDO
3523!
3524!--    Natural surfaces
3525       surf_s = surf_lsm_h%start_index(j,i)
3526       surf_e = surf_lsm_h%end_index(j,i)
3527       DO  m = surf_s, surf_e
3528          k = surf_lsm_h%k(m)
3529
3530          dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
3531          dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
3532
3533       ENDDO
3534!
3535!--    Urban surfaces
3536       surf_s = surf_usm_h%start_index(j,i)
3537       surf_e = surf_usm_h%end_index(j,i)
3538       DO  m = surf_s, surf_e
3539          k = surf_usm_h%k(m)
3540
3541          dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
3542          dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
3543
3544       ENDDO
3545!
3546!--    Compute gradients at downward-facing walls, only for
3547!--    non-natural default surfaces
3548       surf_s = surf_def_h(1)%start_index(j,i)
3549       surf_e = surf_def_h(1)%end_index(j,i)
3550       DO  m = surf_s, surf_e
3551          k = surf_def_h(1)%k(m)
3552
3553          dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
3554          dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
3555
3556       ENDDO
3557
3558    ENDIF
3559
3560    DO  k = nzb+1, nzt
3561
3562       def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
3563                        dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 +           &
3564                        dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +           &
3565             2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) +              &
3566                        dwdy(k)*dvdz(k) )
3567
3568       IF ( def < 0.0_wp )  def = 0.0_wp
3569
3570       flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),flag_nr) )
3571
3572       IF ( .NOT. diss_production )  THEN
3573
3574!--       Compute tendency for TKE-production from shear
3575          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
3576
3577       ELSE
3578
3579!--       RANS mode: Compute tendency for dissipation-rate-production from shear
3580          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
3581                        diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
3582
3583       ENDIF
3584
3585    ENDDO
3586
3587!
3588!-- If required, calculate TKE production by buoyancy
3589    IF ( .NOT. neutral )  THEN
3590
3591       IF ( .NOT. humidity )  THEN
3592
3593          IF ( ocean_mode )  THEN
3594!
3595!--          So far in the ocean no special treatment of density flux
3596!--          in the bottom and top surface layer
3597             DO  k = nzb+1, nzt
3598                tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
3599             ENDDO
3600!
3601!--          Treatment of near-surface grid points, at up- and down-
3602!--          ward facing surfaces
3603             IF ( use_surface_fluxes )  THEN
3604                DO  l = 0, 1
3605                   surf_s = surf_def_h(l)%start_index(j,i)
3606                   surf_e = surf_def_h(l)%end_index(j,i)
3607                   DO  m = surf_s, surf_e
3608                      k = surf_def_h(l)%k(m)
3609                      tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3610                   ENDDO
3611                ENDDO
3612             ENDIF
3613
3614             IF ( use_top_fluxes )  THEN
3615                surf_s = surf_def_h(2)%start_index(j,i)
3616                surf_e = surf_def_h(2)%end_index(j,i)
3617                DO  m = surf_s, surf_e
3618                   k = surf_def_h(2)%k(m)
3619                   tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3620                ENDDO
3621             ENDIF
3622
3623             IF ( .NOT. diss_production )  THEN
3624
3625!--             Compute tendency for TKE-production from shear
3626                DO  k = nzb+1, nzt
3627                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3628                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3629                                 MERGE( rho_reference, prho(k,j,i),       &
3630                                        use_single_reference_value ) )
3631                ENDDO
3632
3633             ELSE
3634
3635!--             RANS mode: Compute tendency for dissipation-rate-production from shear
3636                DO  k = nzb+1, nzt
3637                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3638                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3639                                 MERGE( rho_reference, prho(k,j,i),       &
3640                                        use_single_reference_value ) ) *  &
3641                                 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
3642                                 c_3
3643                ENDDO
3644
3645             ENDIF
3646
3647
3648          ELSE ! or IF ( .NOT. ocean_mode )  THEN
3649
3650             DO  k = nzb+1, nzt
3651                tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
3652             ENDDO
3653
3654             IF ( use_surface_fluxes )  THEN
3655!
3656!--             Default surfaces, up- and downward-facing
3657                DO  l = 0, 1
3658                   surf_s = surf_def_h(l)%start_index(j,i)
3659                   surf_e = surf_def_h(l)%end_index(j,i)
3660                   DO  m = surf_s, surf_e
3661                      k = surf_def_h(l)%k(m)
3662                      tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3663                   ENDDO
3664                ENDDO
3665!
3666!--             Natural surfaces
3667                surf_s = surf_lsm_h%start_index(j,i)
3668                surf_e = surf_lsm_h%end_index(j,i)
3669                DO  m = surf_s, surf_e
3670                   k = surf_lsm_h%k(m)
3671                   tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
3672                ENDDO
3673!
3674!--             Urban surfaces
3675                surf_s = surf_usm_h%start_index(j,i)
3676                surf_e = surf_usm_h%end_index(j,i)
3677                DO  m = surf_s, surf_e
3678                   k = surf_usm_h%k(m)
3679                   tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
3680                ENDDO
3681             ENDIF
3682
3683             IF ( use_top_fluxes )  THEN
3684                surf_s = surf_def_h(2)%start_index(j,i)
3685                surf_e = surf_def_h(2)%end_index(j,i)
3686                DO  m = surf_s, surf_e
3687                   k = surf_def_h(2)%k(m)
3688                   tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3689                ENDDO
3690             ENDIF
3691
3692             IF ( .NOT. diss_production )  THEN
3693
3694!--             Compute tendency for TKE-production from shear
3695                DO  k = nzb+1, nzt
3696                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3697                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3698                                 MERGE( pt_reference, pt(k,j,i),          &
3699                                        use_single_reference_value ) )
3700                ENDDO
3701
3702             ELSE
3703
3704!--             RANS mode: Compute tendency for dissipation-rate-production from shear
3705                DO  k = nzb+1, nzt
3706                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3707                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3708                                 MERGE( pt_reference, pt(k,j,i),          &
3709                                        use_single_reference_value ) ) *  &
3710                                 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
3711                                 c_3
3712                ENDDO
3713
3714             ENDIF
3715
3716          ENDIF ! from IF ( .NOT. ocean_mode )
3717
3718       ELSE ! or IF ( humidity )  THEN
3719
3720          DO  k = nzb+1, nzt
3721
3722             IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3723                k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3724                k2 = 0.61_wp * pt(k,j,i)
3725                tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
3726                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3727                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
3728                                ) * dd2zu(k)
3729             ELSE IF ( bulk_cloud_model )  THEN
3730                IF ( ql(k,j,i) == 0.0_wp )  THEN
3731                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3732                   k2 = 0.61_wp * pt(k,j,i)
3733                ELSE
3734                   theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3735                   temp  = theta * exner(k)
3736                   k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
3737                                 ( q(k,j,i) - ql(k,j,i) ) *              &
3738                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
3739                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
3740                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3741                   k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3742                ENDIF
3743                tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
3744                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3745                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
3746                                ) * dd2zu(k)
3747             ELSE IF ( cloud_droplets )  THEN
3748                k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3749                k2 = 0.61_wp * pt(k,j,i)
3750                tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3751                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3752                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
3753                                  pt(k,j,i) * ( ql(k+1,j,i) -            &
3754                                  ql(k-1,j,i) ) ) * dd2zu(k)
3755             ENDIF
3756
3757          ENDDO
3758
3759          IF ( use_surface_fluxes )  THEN
3760!
3761!--          Treat horizontal default surfaces
3762             DO  l = 0, 1
3763                surf_s = surf_def_h(l)%start_index(j,i)
3764                surf_e = surf_def_h(l)%end_index(j,i)
3765                DO  m = surf_s, surf_e
3766                   k = surf_def_h(l)%k(m)
3767
3768                   IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3769                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3770                      k2 = 0.61_wp * pt(k,j,i)
3771                   ELSE IF ( bulk_cloud_model )  THEN
3772                      IF ( ql(k,j,i) == 0.0_wp )  THEN
3773                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3774                         k2 = 0.61_wp * pt(k,j,i)
3775                      ELSE
3776                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3777                         temp  = theta * exner(k)
3778                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
3779                                    ( q(k,j,i) - ql(k,j,i) ) *           &
3780                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3781                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3782                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3783                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3784                      ENDIF
3785                   ELSE IF ( cloud_droplets )  THEN
3786                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3787                      k2 = 0.61_wp * pt(k,j,i)
3788                   ENDIF
3789
3790                   tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
3791                                   k2 * surf_def_h(l)%qsws(m)            &
3792                                 ) * drho_air_zw(k-1)
3793                ENDDO
3794             ENDDO
3795!
3796!--          Treat horizontal natural surfaces
3797             surf_s = surf_lsm_h%start_index(j,i)
3798             surf_e = surf_lsm_h%end_index(j,i)
3799             DO  m = surf_s, surf_e
3800                k = surf_lsm_h%k(m)
3801
3802                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3803                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3804                   k2 = 0.61_wp * pt(k,j,i)
3805                ELSE IF ( bulk_cloud_model )  THEN
3806                   IF ( ql(k,j,i) == 0.0_wp )  THEN
3807                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3808                      k2 = 0.61_wp * pt(k,j,i)
3809                   ELSE
3810                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3811                      temp  = theta * exner(k)
3812                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3813                                    ( q(k,j,i) - ql(k,j,i) ) *           &
3814                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3815                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3816                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3817                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3818                   ENDIF
3819                ELSE IF ( cloud_droplets )  THEN
3820                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3821                   k2 = 0.61_wp * pt(k,j,i)
3822                ENDIF
3823
3824                tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
3825                                k2 * surf_lsm_h%qsws(m)                  &
3826                              ) * drho_air_zw(k-1)
3827             ENDDO
3828!
3829!--          Treat horizontal urban surfaces
3830             surf_s = surf_usm_h%start_index(j,i)
3831             surf_e = surf_usm_h%end_index(j,i)
3832             DO  m = surf_s, surf_e
3833                k = surf_usm_h%k(m)
3834
3835                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3836                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3837                   k2 = 0.61_wp * pt(k,j,i)
3838                ELSE IF ( bulk_cloud_model )  THEN
3839                   IF ( ql(k,j,i) == 0.0_wp )  THEN
3840                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3841                      k2 = 0.61_wp * pt(k,j,i)
3842                   ELSE
3843                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3844                      temp  = theta * exner(k)
3845                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3846                                    ( q(k,j,i) - ql(k,j,i) ) *           &
3847                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3848                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3849                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3850                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3851                   ENDIF
3852                ELSE IF ( cloud_droplets )  THEN
3853                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3854                   k2 = 0.61_wp * pt(k,j,i)
3855                ENDIF
3856
3857                tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
3858                                k2 * surf_usm_h%qsws(m)                  &
3859                              ) * drho_air_zw(k-1)
3860             ENDDO
3861
3862          ENDIF ! from IF ( use_surface_fluxes )  THEN
3863
3864          IF ( use_top_fluxes )  THEN
3865
3866             surf_s = surf_def_h(2)%start_index(j,i)
3867             surf_e = surf_def_h(2)%end_index(j,i)
3868             DO  m = surf_s, surf_e
3869                k = surf_def_h(2)%k(m)
3870
3871                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3872                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3873                   k2 = 0.61_wp * pt(k,j,i)
3874                ELSE IF ( bulk_cloud_model )  THEN
3875                   IF ( ql(k,j,i) == 0.0_wp )  THEN
3876                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3877                      k2 = 0.61_wp * pt(k,j,i)
3878                   ELSE
3879                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3880                      temp  = theta * exner(k)
3881                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3882                                 ( q(k,j,i) - ql(k,j,i) ) *              &
3883                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
3884                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
3885                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3886                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3887                   ENDIF
3888                ELSE IF ( cloud_droplets )  THEN
3889                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3890                   k2 = 0.61_wp * pt(k,j,i)
3891                ENDIF
3892
3893                tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
3894                                k2 * surf_def_h(2)%qsws(m)               &
3895                              ) * drho_air_zw(k)
3896
3897             ENDDO
3898
3899          ENDIF ! from IF ( use_top_fluxes )  THEN
3900
3901          IF ( .NOT. diss_production )  THEN
3902
3903!--          Compute tendency for TKE-production from shear
3904             DO  k = nzb+1, nzt
3905                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3906                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3907                              MERGE( vpt_reference, vpt(k,j,i),          &
3908                                     use_single_reference_value ) )
3909             ENDDO
3910
3911          ELSE
3912
3913!--          RANS mode: Compute tendency for dissipation-rate-production from shear
3914             DO  k = nzb+1, nzt
3915                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
3916                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
3917                              MERGE( vpt_reference, vpt(k,j,i),          &
3918                                     use_single_reference_value ) ) *    &
3919                              diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
3920                              c_3
3921             ENDDO
3922
3923          ENDIF
3924
3925       ENDIF
3926
3927    ENDIF
3928
3929 END SUBROUTINE production_e_ij
3930
3931
3932!------------------------------------------------------------------------------!
3933! Description:
3934! ------------
3935!> Diffusion and dissipation terms for the TKE.
3936!> Vector-optimized version
3937!> @todo Try to avoid the usage of the 3d-array 'diss' where possible (case les
3938!>       and rans_tke_l if not wang_kernel, use_sgs_for_particles, or
3939!>       collision_turbulence).
3940!------------------------------------------------------------------------------!
3941 SUBROUTINE diffusion_e( var, var_reference )
3942
3943    USE arrays_3d,                                                             &
3944        ONLY:  dd2zu, ddzu, ddzw, drho_air, rho_air_zw
3945
3946    USE control_parameters,                                                    &
3947        ONLY:  atmos_ocean_sign, use_single_reference_value
3948
3949    USE grid_variables,                                                        &
3950        ONLY:  ddx2, ddy2
3951
3952    USE bulk_cloud_model_mod,                                                  &
3953        ONLY:  collision_turbulence
3954
3955    USE particle_attributes,                                                   &
3956        ONLY:  use_sgs_for_particles, wang_kernel
3957
3958    IMPLICIT NONE
3959
3960    INTEGER(iwp) ::  i              !< running index x direction
3961    INTEGER(iwp) ::  j              !< running index y direction
3962    INTEGER(iwp) ::  k              !< running index z direction
3963    INTEGER(iwp) ::  m              !< running index surface elements
3964
3965    REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
3966    REAL(wp)     ::  dvar_dz        !< vertical gradient of var
3967    REAL(wp)     ::  l              !< mixing length
3968    REAL(wp)     ::  var_reference  !< reference temperature
3969
3970    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
3971    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
3972
3973    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
3974
3975!
3976!-- Calculate the dissipation
3977    IF ( les_dynamic .OR. les_mw )  THEN
3978
3979       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
3980       !$ACC PRIVATE(l, l_stable, dvar_dz) &
3981       !$ACC PRESENT(diss, e, var, wall_flags_total_0) &
3982       !$ACC PRESENT(dd2zu, l_grid, l_wall)
3983       DO  i = nxl, nxr
3984          DO  j = nys, nyn
3985             !$ACC LOOP PRIVATE(k)
3986             DO  k = nzb+1, nzt
3987
3988                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
3989                IF ( dvar_dz > 0.0_wp ) THEN
3990                   IF ( use_single_reference_value )  THEN
3991                     l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
3992                                 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
3993                   ELSE
3994                     l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
3995                                 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
3996                   ENDIF
3997                ELSE
3998                   l_stable(k) = l_grid(k)
3999                ENDIF
4000
4001             ENDDO
4002
4003             !$ACC LOOP PRIVATE(k)
4004             !DIR$ IVDEP
4005             DO  k = nzb+1, nzt
4006
4007                l  = MIN( l_wall(k,j,i), l_stable(k) )
4008
4009                diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )                &
4010                              * e(k,j,i) * SQRT( e(k,j,i) ) / l                        &
4011                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
4012
4013             ENDDO
4014
4015          ENDDO
4016       ENDDO
4017
4018    ELSEIF ( rans_tke_l )  THEN
4019
4020      !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
4021      !$ACC PRIVATE(l_stable, duv2_dz2, rif, dvar_dz) &
4022      !$ACC PRESENT(diss, e, u, v, var, wall_flags_total_0) &
4023      !$ACC PRESENT(dd2zu, l_black, l_wall)
4024       DO  i = nxl, nxr
4025          DO  j = nys, nyn
4026!
4027!--          Calculate Richardson-flux number
4028             IF ( use_single_reference_value )  THEN
4029                !$ACC LOOP PRIVATE(k)
4030                DO  k = nzb+1, nzt
4031
4032                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4033
4034                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i)