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

Last change on this file since 2918 was 2918, checked in by gronemeier, 3 years ago

moved initialization of mixing length to turbulence_closure_mod; consider walls in calculation of ml in rans-mode

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