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

Last change on this file since 3182 was 3182, checked in by suehring, 3 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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