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

Last change on this file since 3359 was 3359, checked in by knoop, 3 years ago

Restructured loops and ifs in production_e to ease vectorization on GPUs

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