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

Last change on this file since 3045 was 3045, checked in by Giersch, 3 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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