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

Last change on this file since 3430 was 3430, checked in by maronga, 3 years ago

added support for buildings with dynamic sgs model

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