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

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

Read information from statitic driver for resolved vegetation independently from land- or urban-surface model

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