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

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

bugfix: set standard value for rans_const_sigma

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