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

Last change on this file since 3775 was 3775, checked in by gronemeier, 2 years ago

removed unused variables

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