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

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

merge with branch rans

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