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

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

Bugfix: wrong use of surf_lsm_h changed to surf_usm_h

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