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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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