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

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

calculate diss production same in vector and cache optimization; move boundary condition for e and diss to boundary_conds

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