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

Last change on this file since 3294 was 3294, checked in by raasch, 3 years ago

modularization of the ocean code

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