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

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

Modularization of all bulk cloud physics code components

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