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

Last change on this file since 3004 was 3004, checked in by Giersch, 7 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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