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

Last change on this file since 2718 was 2718, checked in by maronga, 4 years ago

deleting of deprecated files; headers updated where needed

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