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

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

implementation of dynamic sgs model

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