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

Last change on this file since 3724 was 3724, checked in by kanani, 2 years ago

Correct double-used log_point_s units (bulk_cloud_model_mod, time_integration, turbulence_closure_mod)

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