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

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

merge with branch rans: update of rans mode and data output

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