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

Last change on this file since 4170 was 4170, checked in by gronemeier, 2 years ago

changes in turbulence_closure_mod:

  • add performance optimizations according to K. Ketelsen to diffusion_e and tcm_diffusivities_default
  • bugfix in calculating l_wall for vertical walls
  • bugfix in using l_wall in initialization (consider wall_adjustment_factor)
  • always initialize diss and save the dissipation to that array

related changes in time_integration:

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