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

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