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

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

write out preprocessor directives; remove tailing whitespaces

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