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

Last change on this file since 2938 was 2938, checked in by suehring, 7 years ago

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

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