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

Last change on this file since 3014 was 3014, checked in by maronga, 3 years ago

series of bugfixes

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