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

Last change on this file since 4329 was 4329, checked in by motisi, 18 months ago

Renamed wall_flags_0 to wall_flags_static_0

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