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

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

OpenACC port for SPEC

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