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

Last change on this file since 3121 was 3121, checked in by gronemeier, 3 years ago

consider wall function for Km within production_e_init

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