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

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

renamed tcm_prognostic to tcm_prognostic_equations

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