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

Last change on this file since 3183 was 3183, checked in by suehring, 3 years ago

last commit documented

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