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

Last change on this file since 3684 was 3684, checked in by knoop, 2 years ago

Enabled module_interface_actions in time_integration and prognostic_equations

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