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

Last change on this file since 3646 was 3646, checked in by kanani, 2 years ago

Bugfix: replace simulated_time by time_since_reference_point where required

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