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

Last change on this file since 3398 was 3398, checked in by knoop, 3 years ago

Refactored production_e and production_e_ij (removed code redundancy)

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