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

Last change on this file since 4048 was 4048, checked in by knoop, 2 years ago

Moved turbulence_closure_mod calls into module_interface

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