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

Last change on this file since 2796 was 2764, checked in by gronemeier, 6 years ago

Bugfix: remove duplicate SAVE statements

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