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

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

calculate km according to MOST within the surface layer

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