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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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