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

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

Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)

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