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

Last change on this file since 2701 was 2701, checked in by suehring, 5 years ago

changes from last commit documented

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