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

Last change on this file since 3300 was 3300, checked in by gronemeier, 3 years ago

bugfix: rename of variable

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