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

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

removed debug output; removed rans_mode from namelist; altered order of check_parameter-calls of modules

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