source: palm/trunk/SOURCE/production_e.f90 @ 2604

Last change on this file since 2604 was 2508, checked in by suehring, 7 years ago

Bugfixes in SGS-TKE buoyancy production; revised initialization of vertical-gradient levels in case of ocean runs

  • Property svn:keywords set to Id
File size: 82.0 KB
Line 
1!> @file production_e.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
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 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: production_e.f90 2508 2017-10-02 08:57:09Z schwenkel $
27! - Bugfix in buoyancy production term, wrong base state was set.
28! - Consider use_single_reference_value case if humidity is used.
29! - In case of use_top_fluxes, use correct inverse density at model top
30! - Consider use_surface_fluxes and use_top_fluxes in ocean case 
31!
32! 2478 2017-09-18 13:37:24Z suehring
33! Bugfix, consider case where no constant-flux layer and no surfaces fluxes
34! are used
35!
36! 2329 2017-08-03 14:24:56Z knoop
37! Bugfix: added division by density as kinematic fluxes are needed
38!
39! 2233 2017-05-30 18:08:54Z suehring
40!
41! 2232 2017-05-30 17:47:52Z suehring
42! Adjustments to new surface concept
43!
44! 2126 2017-01-20 15:54:21Z raasch
45! density in ocean case replaced by potential density
46!
47! 2118 2017-01-17 16:38:49Z raasch
48! OpenACC version of subroutine removed
49!
50! 2031 2016-10-21 15:11:58Z knoop
51! renamed variable rho to rho_ocean
52!
53! 2000 2016-08-20 18:09:15Z knoop
54! Forced header and separation lines into 80 columns
55!
56! 1873 2016-04-18 14:50:06Z maronga
57! Module renamed (removed _mod)
58!
59!
60! 1850 2016-04-08 13:29:27Z maronga
61! Module renamed
62!
63!
64! 1691 2015-10-26 16:17:44Z maronga
65! Renamed prandtl_layer to constant_flux_layer.
66!
67! 1682 2015-10-07 23:56:08Z knoop
68! Code annotations made doxygen readable
69!
70! 1374 2014-04-25 12:55:07Z raasch
71! nzb_s_outer removed from acc-present-list
72!
73! 1353 2014-04-08 15:21:23Z heinze
74! REAL constants provided with KIND-attribute
75!
76! 1342 2014-03-26 17:04:47Z kanani
77! REAL constants defined as wp-kind
78!
79! 1320 2014-03-20 08:40:49Z raasch
80! ONLY-attribute added to USE-statements,
81! kind-parameters added to all INTEGER and REAL declaration statements,
82! kinds are defined in new module kinds,
83! old module precision_kind is removed,
84! revision history before 2012 removed,
85! comment fields (!:) to be used for variable explanations added to
86! all variable declaration statements
87!
88! 1257 2013-11-08 15:18:40Z raasch
89! openacc loop and loop vector clauses removed, declare create moved after
90! the FORTRAN declaration statement
91!
92! 1179 2013-06-14 05:57:58Z raasch
93! use_reference renamed use_single_reference_value
94!
95! 1128 2013-04-12 06:19:32Z raasch
96! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
97! j_north
98!
99! 1036 2012-10-22 13:43:42Z raasch
100! code put under GPL (PALM 3.9)
101!
102! 1015 2012-09-27 09:23:24Z raasch
103! accelerator version (*_acc) added
104!
105! 1007 2012-09-19 14:30:36Z franke
106! Bugfix: calculation of buoyancy production has to consider the liquid water
107! mixing ratio in case of cloud droplets
108!
109! 940 2012-07-09 14:31:00Z raasch
110! TKE production by buoyancy can be switched off in case of runs with pure
111! neutral stratification
112!
113! Revision 1.1  1997/09/19 07:45:35  raasch
114! Initial revision
115!
116!
117! Description:
118! ------------
119!> Production terms (shear + buoyancy) of the TKE.
120!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
121!>          not considered well!
122!------------------------------------------------------------------------------!
123 MODULE production_e_mod
124
125    USE kinds
126
127    PRIVATE
128    PUBLIC production_e, production_e_init
129
130    INTERFACE production_e
131       MODULE PROCEDURE production_e
132       MODULE PROCEDURE production_e_ij
133    END INTERFACE production_e
134   
135    INTERFACE production_e_init
136       MODULE PROCEDURE production_e_init
137    END INTERFACE production_e_init
138 
139 CONTAINS
140
141
142!------------------------------------------------------------------------------!
143! Description:
144! ------------
145!> Call for all grid points
146!------------------------------------------------------------------------------!
147    SUBROUTINE production_e
148
149       USE arrays_3d,                                                          &
150           ONLY:  ddzw, dd2zu, drho_air_zw, kh, km, prho, pt, q, ql, tend, u,  &
151                  v, vpt, w
152
153       USE cloud_parameters,                                                   &
154           ONLY:  l_d_cp, l_d_r, pt_d_t, t_d_pt
155
156       USE control_parameters,                                                 &
157           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
158                  humidity, kappa, neutral, ocean, pt_reference,               &
159                  rho_reference, use_single_reference_value,                   &
160                  use_surface_fluxes, use_top_fluxes, vpt_reference
161
162       USE grid_variables,                                                     &
163           ONLY:  ddx, dx, ddy, dy
164
165       USE indices,                                                            &
166           ONLY:  nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
167
168       USE surface_mod,                                                        &
169           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
170                   surf_usm_v
171
172       IMPLICIT NONE
173
174       INTEGER(iwp) ::  i       !< running index x-direction
175       INTEGER(iwp) ::  j       !< running index y-direction
176       INTEGER(iwp) ::  k       !< running index z-direction
177       INTEGER(iwp) ::  l       !< running index for different surface type orientation
178       INTEGER(iwp) ::  m       !< running index surface elements
179       INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
180       INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
181
182       REAL(wp)     ::  def         !<
183       REAL(wp)     ::  flag        !< flag to mask topography
184       REAL(wp)     ::  k1          !<
185       REAL(wp)     ::  k2          !<
186       REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
187       REAL(wp)     ::  theta       !<
188       REAL(wp)     ::  temp        !<
189       REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
190       REAL(wp)     ::  usvs        !< momentum flux u"v"
191       REAL(wp)     ::  vsus        !< momentum flux v"u"
192       REAL(wp)     ::  wsus        !< momentum flux w"u"
193       REAL(wp)     ::  wsvs        !< momentum flux w"v"
194
195       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dudx  !< Gradient of u-component in x-direction
196       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dudy  !< Gradient of u-component in y-direction
197       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dudz  !< Gradient of u-component in z-direction
198       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dvdx  !< Gradient of v-component in x-direction
199       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dvdy  !< Gradient of v-component in y-direction
200       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dvdz  !< Gradient of v-component in z-direction
201       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dwdx  !< Gradient of w-component in x-direction
202       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dwdy  !< Gradient of w-component in y-direction
203       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dwdz  !< Gradient of w-component in z-direction
204
205       DO  i = nxl, nxr
206
207          IF ( constant_flux_layer )  THEN
208
209!
210!--          Calculate TKE production by shear. Calculate gradients at all grid
211!--          points first, gradients at surface-bounded grid points will be
212!--          overwritten further below.
213             DO  j = nys, nyn
214                DO  k = nzb+1, nzt
215
216                   dudx(k,j) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
217                   dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -         &
218                                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
219                   dudz(k,j) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -         &
220                                           u(k-1,j,i) - u(k-1,j,i+1) ) *       &
221                                                            dd2zu(k)
222   
223                   dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -         &
224                                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
225                   dvdy(k,j) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
226                   dvdz(k,j) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -         &
227                                           v(k-1,j,i) - v(k-1,j+1,i) ) *       &
228                                                            dd2zu(k)
229
230                   dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -         &
231                                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
232                   dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -         &
233                                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
234                   dwdz(k,j) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
235
236                ENDDO
237             ENDDO
238
239!
240!--          Position beneath wall
241!--          (2) - Will allways be executed.
242!--          'bottom and wall: use u_0,v_0 and wall functions'
243             DO  j = nys, nyn
244!
245!--             Compute gradients at north- and south-facing surfaces.
246!--             First, for default surfaces, then for urban surfaces.
247!--             Note, so far no natural vertical surfaces implemented
248                DO  l = 0, 1
249                   surf_s = surf_def_v(l)%start_index(j,i)
250                   surf_e = surf_def_v(l)%end_index(j,i)
251                   DO  m = surf_s, surf_e
252                      k           = surf_def_v(l)%k(m)
253                      usvs        = surf_def_v(l)%mom_flux_tke(0,m)
254                      wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
255     
256                      km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp      &
257                                      * 0.5_wp * dy
258!
259!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
260                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
261                                        BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
262                      dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
263                      dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
264                   ENDDO
265!
266!--                Natural surfaces
267                   surf_s = surf_lsm_v(l)%start_index(j,i)
268                   surf_e = surf_lsm_v(l)%end_index(j,i)
269                   DO  m = surf_s, surf_e
270                      k           = surf_lsm_v(l)%k(m)
271                      usvs        = surf_lsm_v(l)%mom_flux_tke(0,m)
272                      wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
273     
274                      km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp      &
275                                      * 0.5_wp * dy
276!
277!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
278                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
279                                        BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
280                      dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
281                      dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
282                   ENDDO 
283!
284!--                Urban surfaces
285                   surf_s = surf_usm_v(l)%start_index(j,i)
286                   surf_e = surf_usm_v(l)%end_index(j,i)
287                   DO  m = surf_s, surf_e
288                      k           = surf_usm_v(l)%k(m)
289                      usvs        = surf_usm_v(l)%mom_flux_tke(0,m)
290                      wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
291     
292                      km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp      &
293                                      * 0.5_wp * dy
294!
295!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
296                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
297                                        BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
298                      dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
299                      dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
300                   ENDDO 
301                ENDDO
302!
303!--             Compute gradients at east- and west-facing walls
304                DO  l = 2, 3
305                   surf_s = surf_def_v(l)%start_index(j,i)
306                   surf_e = surf_def_v(l)%end_index(j,i)
307                   DO  m = surf_s, surf_e
308                      k     = surf_def_v(l)%k(m)
309                      vsus  = surf_def_v(l)%mom_flux_tke(0,m)
310                      wsus  = surf_def_v(l)%mom_flux_tke(1,m)
311
312                      km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp      &
313                                         * 0.5_wp * dx
314!
315!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
316                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
317                                        BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
318                      dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
319                      dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
320                   ENDDO 
321!
322!--                Natural surfaces                   
323                   surf_s = surf_lsm_v(l)%start_index(j,i)
324                   surf_e = surf_lsm_v(l)%end_index(j,i)
325                   DO  m = surf_s, surf_e
326                      k     = surf_lsm_v(l)%k(m)
327                      vsus  = surf_lsm_v(l)%mom_flux_tke(0,m)
328                      wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
329
330                      km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp      &
331                                         * 0.5_wp * dx
332!
333!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
334                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
335                                        BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
336                      dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
337                      dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
338                   ENDDO   
339!
340!--                Urban surfaces                   
341                   surf_s = surf_usm_v(l)%start_index(j,i)
342                   surf_e = surf_usm_v(l)%end_index(j,i)
343                   DO  m = surf_s, surf_e
344                      k     = surf_usm_v(l)%k(m)
345                      vsus  = surf_usm_v(l)%mom_flux_tke(0,m)
346                      wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
347
348                      km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp      &
349                                         * 0.5_wp * dx
350!
351!--                   -1.0 for right-facing wall, 1.0 for left-facing wall
352                      sign_dir = MERGE( 1.0_wp, -1.0_wp,                       &
353                                        BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
354                      dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
355                      dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
356                   ENDDO 
357                ENDDO
358!
359!--             Compute gradients at upward-facing surfaces
360                surf_s = surf_def_h(0)%start_index(j,i)
361                surf_e = surf_def_h(0)%end_index(j,i)
362                DO  m = surf_s, surf_e
363                   k = surf_def_h(0)%k(m)
364!
365!--                Please note, actually, an interpolation of u_0 and v_0
366!--                onto the grid center would be required. However, this
367!--                would require several data transfers between 2D-grid and
368!--                wall type. The effect of this missing interpolation is
369!--                negligible. (See also production_e_init).
370                   dudz(k,j) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)   
371                   dvdz(k,j) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
372           
373                ENDDO
374!
375!--             Natural surfaces
376                surf_s = surf_lsm_h%start_index(j,i)
377                surf_e = surf_lsm_h%end_index(j,i)
378                DO  m = surf_s, surf_e
379                   k = surf_lsm_h%k(m)
380
381                   dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
382                   dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
383           
384                ENDDO
385!
386!--             Urban surfaces
387                surf_s = surf_usm_h%start_index(j,i)
388                surf_e = surf_usm_h%end_index(j,i)
389                DO  m = surf_s, surf_e
390                   k = surf_usm_h%k(m)
391
392                   dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
393                   dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
394           
395                ENDDO
396!
397!--             Compute gradients at downward-facing walls, only for
398!--             non-natural default surfaces
399                surf_s = surf_def_h(1)%start_index(j,i)
400                surf_e = surf_def_h(1)%end_index(j,i)
401                DO  m = surf_s, surf_e
402                   k = surf_def_h(1)%k(m)
403
404                   dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
405                   dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
406
407                ENDDO
408             ENDDO
409
410             DO  j = nys, nyn
411                DO  k = nzb+1, nzt
412
413                   def = 2.0_wp * ( dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 ) + &
414                                    dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 +   &
415                                    dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 +   &
416                         2.0_wp * ( dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j)  +   &
417                                    dwdy(k,j)*dvdz(k,j) )
418
419                   IF ( def < 0.0_wp )  def = 0.0_wp
420
421                   flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
422
423                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
424
425                ENDDO
426             ENDDO
427
428          ELSE
429
430             DO  j = nys, nyn
431!
432!--             Calculate TKE production by shear. Here, no additional
433!--             wall-bounded code is considered.
434!--             Why?
435                DO  k = nzb+1, nzt
436
437                   dudx(k,j)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
438                   dudy(k,j)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -        &
439                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
440                   dudz(k,j)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -        &
441                                            u(k-1,j,i) - u(k-1,j,i+1) ) *      &
442                                                                        dd2zu(k)
443
444                   dvdx(k,j)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -        &
445                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
446                   dvdy(k,j)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
447                   dvdz(k,j)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -        &
448                                            v(k-1,j,i) - v(k-1,j+1,i) ) *      &
449                                                                        dd2zu(k)
450
451                   dwdx(k,j)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -        &
452                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
453                   dwdy(k,j)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -        &
454                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
455                   dwdz(k,j)  =           ( w(k,j,i)   - w(k-1,j,i)   ) *      &
456                                                                        ddzw(k)
457     
458                   def = 2.0_wp * (                                            &
459                                dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2     &
460                                  ) +                                          &
461                                dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 +   &
462                                dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 +   &
463                         2.0_wp * (                                            &
464                                dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j)  +   &
465                                dwdy(k,j)*dvdz(k,j)                            &
466                                  )
467
468                   IF ( def < 0.0_wp )  def = 0.0_wp
469
470                   flag  = MERGE( 1.0_wp, 0.0_wp,                              &
471                                  BTEST( wall_flags_0(k,j,i), 29 ) ) 
472                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
473     
474                ENDDO
475             ENDDO
476
477          ENDIF
478
479!
480!--       If required, calculate TKE production by buoyancy
481          IF ( .NOT. neutral )  THEN
482
483             IF ( .NOT. humidity )  THEN
484
485                IF ( ocean )  THEN
486!
487!--                So far in the ocean no special treatment of density flux
488!--                in the bottom and top surface layer
489                   DO  j = nys, nyn
490                      DO  k = nzb+1, nzt
491                         tend(k,j,i) = tend(k,j,i) +                           &
492                                       kh(k,j,i) * g /                         &
493                              MERGE( rho_reference, prho(k,j,i),               &
494                                     use_single_reference_value ) *            &
495                                       ( prho(k+1,j,i) - prho(k-1,j,i) ) *     &
496                                       dd2zu(k) *                              &
497                                   MERGE( 1.0_wp, 0.0_wp,                      &
498                                          BTEST( wall_flags_0(k,j,i), 30 )     &
499                                        )                            *         &
500                                   MERGE( 1.0_wp, 0.0_wp,                      &
501                                          BTEST( wall_flags_0(k,j,i), 9 )      &
502                                        ) 
503                      ENDDO
504!
505!--                   Treatment of near-surface grid points, at up- and down-
506!--                   ward facing surfaces
507                      IF ( use_surface_fluxes )  THEN
508                         DO  l = 0, 1
509                            surf_s = surf_def_h(l)%start_index(j,i)
510                            surf_e = surf_def_h(l)%end_index(j,i)
511                            DO  m = surf_s, surf_e
512                               k = surf_def_h(l)%k(m)
513                               tend(k,j,i) = tend(k,j,i) + g /                 &
514                                         MERGE( rho_reference, prho(k,j,i),    &
515                                                use_single_reference_value ) * &
516                                         drho_air_zw(k-1) *                    &
517                                         surf_def_h(l)%shf(m)
518                            ENDDO
519                         ENDDO
520
521                      ENDIF
522
523                      IF ( use_top_fluxes )  THEN
524                         surf_s = surf_def_h(2)%start_index(j,i)
525                         surf_e = surf_def_h(2)%end_index(j,i)
526                         DO  m = surf_s, surf_e
527                            k = surf_def_h(2)%k(m)
528                            tend(k,j,i) = tend(k,j,i) + g /                    &
529                                         MERGE( rho_reference, prho(k,j,i),    &
530                                                use_single_reference_value ) * &
531                                         drho_air_zw(k) *                      &
532                                         surf_def_h(2)%shf(m) 
533                         ENDDO
534                      ENDIF
535
536                   ENDDO
537
538                ELSE
539
540                   DO  j = nys, nyn
541                      DO  k = nzb+1, nzt
542!
543!--                      Flag 9 is used to mask top fluxes, flag 30 to mask
544!--                      surface fluxes
545                         tend(k,j,i) = tend(k,j,i) -                           &
546                                       kh(k,j,i) * g /                         &
547                                   MERGE( pt_reference, pt(k,j,i),             &
548                                           use_single_reference_value ) *      &
549                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) *         &
550                                       dd2zu(k)                      *         &
551                                   MERGE( 1.0_wp, 0.0_wp,                      &
552                                          BTEST( wall_flags_0(k,j,i), 30 )     &
553                                        )                            *         &
554                                   MERGE( 1.0_wp, 0.0_wp,                      &
555                                          BTEST( wall_flags_0(k,j,i), 9 )      &
556                                        ) 
557                      ENDDO
558
559                      IF ( use_surface_fluxes )  THEN
560!
561!--                      Default surfaces, up- and downward-facing
562                         DO  l = 0, 1
563                            surf_s = surf_def_h(l)%start_index(j,i)
564                            surf_e = surf_def_h(l)%end_index(j,i)
565                            DO  m = surf_s, surf_e
566                               k = surf_def_h(l)%k(m)
567                               tend(k,j,i) = tend(k,j,i) + g /                 &
568                                    MERGE( pt_reference, pt(k,j,i),            &
569                                           use_single_reference_value )        &
570                                                      * drho_air_zw(k-1)       &
571                                                      * surf_def_h(l)%shf(m)   
572                            ENDDO     
573                         ENDDO
574!
575!--                      Natural surfaces
576                         surf_s = surf_lsm_h%start_index(j,i)
577                         surf_e = surf_lsm_h%end_index(j,i)
578                         DO  m = surf_s, surf_e
579                            k = surf_lsm_h%k(m)
580                            tend(k,j,i) = tend(k,j,i) + g /                    &
581                                    MERGE( pt_reference, pt(k,j,i),            &
582                                           use_single_reference_value )        &
583                                                      * drho_air_zw(k-1)       &
584                                                      * surf_lsm_h%shf(m)   
585                         ENDDO
586!
587!--                      Urban surfaces
588                         surf_s = surf_usm_h%start_index(j,i)
589                         surf_e = surf_usm_h%end_index(j,i)
590                         DO  m = surf_s, surf_e
591                            k = surf_usm_h%k(m)
592                            tend(k,j,i) = tend(k,j,i) + g /                    &
593                                    MERGE( pt_reference, pt(k,j,i),            &
594                                           use_single_reference_value )        &
595                                                      * drho_air_zw(k-1)       &
596                                                      * surf_usm_h%shf(m)   
597                         ENDDO                         
598                      ENDIF
599
600                      IF ( use_top_fluxes )  THEN
601                         surf_s = surf_def_h(2)%start_index(j,i)
602                         surf_e = surf_def_h(2)%end_index(j,i)
603                         DO  m = surf_s, surf_e
604                            k = surf_def_h(2)%k(m)
605                            tend(k,j,i) = tend(k,j,i) + g /                    &
606                                    MERGE( pt_reference, pt(k,j,i),            &
607                                           use_single_reference_value )        &
608                                                      * drho_air_zw(k)         &
609                                                      * surf_def_h(2)%shf(m) 
610                         ENDDO
611                      ENDIF
612                   ENDDO
613
614                ENDIF
615
616             ELSE
617
618                DO  j = nys, nyn
619
620                   DO  k = nzb+1, nzt
621!
622!--                   Flag 9 is used to mask top fluxes, flag 30 to mask
623!--                   surface fluxes
624                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
625                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
626                         k2 = 0.61_wp * pt(k,j,i)
627                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
628                                         g /                                   &
629                                    MERGE( vpt_reference, vpt(k,j,i),          &
630                                           use_single_reference_value ) *      &
631                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
632                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
633                                         ) * dd2zu(k) *                        &
634                                      MERGE( 1.0_wp, 0.0_wp,                   &
635                                             BTEST( wall_flags_0(k,j,i), 30 )  &
636                                           )          *                        &
637                                      MERGE( 1.0_wp, 0.0_wp,                   &
638                                             BTEST( wall_flags_0(k,j,i), 9 )   &
639                                           )
640                      ELSE IF ( cloud_physics )  THEN
641                         IF ( ql(k,j,i) == 0.0_wp )  THEN
642                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
643                            k2 = 0.61_wp * pt(k,j,i)
644                         ELSE
645                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
646                            temp  = theta * t_d_pt(k)
647                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
648                                          ( q(k,j,i) - ql(k,j,i) ) *           &
649                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
650                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
651                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
652                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
653                         ENDIF
654                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
655                                         g /                                   &
656                                    MERGE( vpt_reference, vpt(k,j,i),          &
657                                           use_single_reference_value ) *      &
658                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
659                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
660                                         ) * dd2zu(k) *                        &
661                                      MERGE( 1.0_wp, 0.0_wp,                   &
662                                             BTEST( wall_flags_0(k,j,i), 30 )  &
663                                           )          *                        &
664                                      MERGE( 1.0_wp, 0.0_wp,                   &
665                                             BTEST( wall_flags_0(k,j,i), 9 )   &
666                                           )
667                      ELSE IF ( cloud_droplets )  THEN
668                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
669                         k2 = 0.61_wp * pt(k,j,i)
670                         tend(k,j,i) = tend(k,j,i) -                           & 
671                                       kh(k,j,i) * g /                         &
672                                    MERGE( vpt_reference, vpt(k,j,i),          &
673                                           use_single_reference_value ) *      &
674                                       ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +   &
675                                         k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -   &
676                                         pt(k,j,i) * ( ql(k+1,j,i) -           &
677                                         ql(k-1,j,i) ) ) * dd2zu(k) *          &
678                                      MERGE( 1.0_wp, 0.0_wp,                   &
679                                             BTEST( wall_flags_0(k,j,i), 30 )  &
680                                           )                        *          &
681                                      MERGE( 1.0_wp, 0.0_wp,                   &
682                                             BTEST( wall_flags_0(k,j,i), 9 )   &
683                                           )
684                      ENDIF
685
686                   ENDDO
687
688                ENDDO
689
690                IF ( use_surface_fluxes )  THEN
691
692                   DO  j = nys, nyn
693!
694!--                   Treat horizontal default surfaces
695                      DO  l = 0, 1
696                         surf_s = surf_def_h(l)%start_index(j,i)
697                         surf_e = surf_def_h(l)%end_index(j,i)
698                         DO  m = surf_s, surf_e
699                            k = surf_def_h(l)%k(m)
700
701                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
702                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
703                               k2 = 0.61_wp * pt(k,j,i)
704                            ELSE IF ( cloud_physics )  THEN
705                               IF ( ql(k,j,i) == 0.0_wp )  THEN
706                                  k1 = 1.0_wp + 0.61_wp * q(k,j,i)
707                                  k2 = 0.61_wp * pt(k,j,i)
708                               ELSE
709                                  theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
710                                  temp  = theta * t_d_pt(k)
711                                  k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *         &
712                                             ( q(k,j,i) - ql(k,j,i) ) *        &
713                                    ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /   &
714                                    ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *     &
715                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
716                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
717                               ENDIF
718                            ELSE IF ( cloud_droplets )  THEN
719                               k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
720                               k2 = 0.61_wp * pt(k,j,i)
721                            ENDIF
722
723                            tend(k,j,i) = tend(k,j,i) + g /                    &
724                                    MERGE( vpt_reference, vpt(k,j,i),          &
725                                           use_single_reference_value ) *      &
726                                               ( k1 * surf_def_h(l)%shf(m) +   &
727                                                 k2 * surf_def_h(l)%qsws(m)    &
728                                               ) * drho_air_zw(k-1)
729                         ENDDO
730                      ENDDO
731!
732!--                   Treat horizontal natural surfaces
733                      surf_s = surf_lsm_h%start_index(j,i)
734                      surf_e = surf_lsm_h%end_index(j,i)
735                      DO  m = surf_s, surf_e
736                         k = surf_lsm_h%k(m)
737
738                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
739                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
740                            k2 = 0.61_wp * pt(k,j,i)
741                         ELSE IF ( cloud_physics )  THEN
742                            IF ( ql(k,j,i) == 0.0_wp )  THEN
743                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
744                               k2 = 0.61_wp * pt(k,j,i)
745                            ELSE
746                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
747                               temp  = theta * t_d_pt(k)
748                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
749                                             ( q(k,j,i) - ql(k,j,i) ) *        &
750                                    ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /   &
751                                    ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *     &
752                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
753                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
754                            ENDIF
755                         ELSE IF ( cloud_droplets )  THEN
756                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
757                            k2 = 0.61_wp * pt(k,j,i)
758                         ENDIF
759
760                         tend(k,j,i) = tend(k,j,i) + g /                       &
761                                    MERGE( vpt_reference, vpt(k,j,i),          &
762                                           use_single_reference_value ) *      &
763                                               ( k1 * surf_lsm_h%shf(m) +      &
764                                                 k2 * surf_lsm_h%qsws(m)       &
765                                               ) * drho_air_zw(k-1)
766                      ENDDO
767!
768!--                   Treat horizontal urban surfaces
769                      surf_s = surf_usm_h%start_index(j,i)
770                      surf_e = surf_usm_h%end_index(j,i)
771                      DO  m = surf_s, surf_e
772                         k = surf_lsm_h%k(m)
773
774                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
775                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
776                            k2 = 0.61_wp * pt(k,j,i)
777                         ELSE IF ( cloud_physics )  THEN
778                            IF ( ql(k,j,i) == 0.0_wp )  THEN
779                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
780                               k2 = 0.61_wp * pt(k,j,i)
781                            ELSE
782                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
783                               temp  = theta * t_d_pt(k)
784                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
785                                             ( q(k,j,i) - ql(k,j,i) ) *        &
786                                    ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /   &
787                                    ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *     &
788                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
789                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
790                            ENDIF
791                         ELSE IF ( cloud_droplets )  THEN
792                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
793                            k2 = 0.61_wp * pt(k,j,i)
794                         ENDIF
795
796                         tend(k,j,i) = tend(k,j,i) + g /                       &
797                                    MERGE( vpt_reference, vpt(k,j,i),          &
798                                           use_single_reference_value ) *      &
799                                               ( k1 * surf_usm_h%shf(m) +      &
800                                                 k2 * surf_usm_h%qsws(m)       &
801                                               ) * drho_air_zw(k-1)
802                      ENDDO
803
804                   ENDDO
805
806                ENDIF
807
808                IF ( use_top_fluxes )  THEN
809
810                   DO  j = nys, nyn
811
812                      surf_s = surf_def_h(2)%start_index(j,i)
813                      surf_e = surf_def_h(2)%end_index(j,i)
814                      DO  m = surf_s, surf_e
815                         k = surf_def_h(2)%k(m)
816
817                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
818                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
819                            k2 = 0.61_wp * pt(k,j,i)
820                         ELSE IF ( cloud_physics )  THEN
821                            IF ( ql(k,j,i) == 0.0_wp )  THEN
822                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
823                               k2 = 0.61_wp * pt(k,j,i)
824                            ELSE
825                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
826                               temp  = theta * t_d_pt(k)
827                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
828                                          ( q(k,j,i) - ql(k,j,i) ) *           &
829                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
830                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
831                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
832                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
833                            ENDIF
834                         ELSE IF ( cloud_droplets )  THEN
835                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
836                            k2 = 0.61_wp * pt(k,j,i)
837                         ENDIF
838
839                         tend(k,j,i) = tend(k,j,i) + g /                       &
840                                    MERGE( vpt_reference, vpt(k,j,i),          &
841                                           use_single_reference_value ) *      &
842                                               ( k1 * surf_def_h(2)%shf(m) +   &
843                                                 k2 * surf_def_h(2)%qsws(m)    &
844                                               ) * drho_air_zw(k)
845
846                      ENDDO
847
848                   ENDDO
849
850                ENDIF
851
852             ENDIF
853
854          ENDIF
855
856       ENDDO
857
858    END SUBROUTINE production_e
859
860
861!------------------------------------------------------------------------------!
862! Description:
863! ------------
864!> Call for grid point i,j
865!------------------------------------------------------------------------------!
866    SUBROUTINE production_e_ij( i, j )
867
868       USE arrays_3d,                                                          &
869           ONLY:  ddzw, dd2zu, drho_air_zw, kh, km, prho, pt, q, ql, tend, u,  &
870                  v, vpt, w
871
872       USE cloud_parameters,                                                   &
873           ONLY:  l_d_cp, l_d_r, pt_d_t, t_d_pt
874
875       USE control_parameters,                                                 &
876           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
877                  humidity, kappa, neutral, ocean, pt_reference,               &
878                  rho_reference, use_single_reference_value,                   &
879                  use_surface_fluxes, use_top_fluxes, vpt_reference
880
881       USE grid_variables,                                                     &
882           ONLY:  ddx, dx, ddy, dy
883
884       USE indices,                                                            &
885           ONLY:  nxl, nxr, nys, nyn, nzb, nzb, nzt, wall_flags_0
886
887       USE surface_mod,                                                        &
888           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
889                   surf_usm_v 
890
891       IMPLICIT NONE
892
893       INTEGER(iwp) ::  i       !< running index x-direction
894       INTEGER(iwp) ::  j       !< running index y-direction
895       INTEGER(iwp) ::  k       !< running index z-direction
896       INTEGER(iwp) ::  l       !< running index for different surface type orientation
897       INTEGER(iwp) ::  m       !< running index surface elements
898       INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
899       INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
900
901       REAL(wp)     ::  def         !<
902       REAL(wp)     ::  flag        !< flag to mask topography
903       REAL(wp)     ::  k1          !<
904       REAL(wp)     ::  k2          !<
905       REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
906       REAL(wp)     ::  theta       !<
907       REAL(wp)     ::  temp        !<
908       REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
909       REAL(wp)     ::  usvs        !< momentum flux u"v"
910       REAL(wp)     ::  vsus        !< momentum flux v"u"
911       REAL(wp)     ::  wsus        !< momentum flux w"u"
912       REAL(wp)     ::  wsvs        !< momentum flux w"v"
913
914
915       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudx        !< Gradient of u-component in x-direction
916       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudy        !< Gradient of u-component in y-direction
917       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudz        !< Gradient of u-component in z-direction
918       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdx        !< Gradient of v-component in x-direction
919       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdy        !< Gradient of v-component in y-direction
920       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdz        !< Gradient of v-component in z-direction
921       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdx        !< Gradient of w-component in x-direction
922       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdy        !< Gradient of w-component in y-direction
923       REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdz        !< Gradient of w-component in z-direction
924
925
926       IF ( constant_flux_layer )  THEN
927!
928!--       Calculate TKE production by shear. Calculate gradients at all grid
929!--       points first, gradients at surface-bounded grid points will be
930!--       overwritten further below.
931          DO  k = nzb+1, nzt
932
933             dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
934             dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
935                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 
936             dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
937                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 
938
939             dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
940                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx     
941             dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
942             dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
943                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 
944
945             dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
946                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     
947             dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
948                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy     
949             dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
950
951          ENDDO
952!
953!--       Compute gradients at north- and south-facing surfaces.
954!--       Note, no vertical natural surfaces so far.
955          DO  l = 0, 1
956!
957!--          Default surfaces
958             surf_s = surf_def_v(l)%start_index(j,i)
959             surf_e = surf_def_v(l)%end_index(j,i)
960             DO  m = surf_s, surf_e
961                k           = surf_def_v(l)%k(m)
962                usvs        = surf_def_v(l)%mom_flux_tke(0,m)
963                wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
964
965                km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp            &
966                                * 0.5_wp * dy
967!
968!--             -1.0 for right-facing wall, 1.0 for left-facing wall
969                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
970                                  BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
971                dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
972                dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
973             ENDDO   
974!
975!--          Natural surfaces
976             surf_s = surf_lsm_v(l)%start_index(j,i)
977             surf_e = surf_lsm_v(l)%end_index(j,i)
978             DO  m = surf_s, surf_e
979                k           = surf_lsm_v(l)%k(m)
980                usvs        = surf_lsm_v(l)%mom_flux_tke(0,m)
981                wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
982
983                km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp            &
984                                * 0.5_wp * dy
985!
986!--             -1.0 for right-facing wall, 1.0 for left-facing wall
987                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
988                                  BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
989                dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
990                dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
991             ENDDO 
992!
993!--          Urban surfaces
994             surf_s = surf_usm_v(l)%start_index(j,i)
995             surf_e = surf_usm_v(l)%end_index(j,i)
996             DO  m = surf_s, surf_e
997                k           = surf_usm_v(l)%k(m)
998                usvs        = surf_usm_v(l)%mom_flux_tke(0,m)
999                wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
1000
1001                km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp            &
1002                                * 0.5_wp * dy
1003!
1004!--             -1.0 for right-facing wall, 1.0 for left-facing wall
1005                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
1006                                  BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
1007                dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
1008                dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
1009             ENDDO 
1010          ENDDO
1011!
1012!--       Compute gradients at east- and west-facing walls
1013          DO  l = 2, 3
1014!
1015!--          Default surfaces
1016             surf_s = surf_def_v(l)%start_index(j,i)
1017             surf_e = surf_def_v(l)%end_index(j,i)
1018             DO  m = surf_s, surf_e
1019                k           = surf_def_v(l)%k(m)
1020                vsus        = surf_def_v(l)%mom_flux_tke(0,m)
1021                wsus        = surf_def_v(l)%mom_flux_tke(1,m)
1022
1023                km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp            &
1024                                   * 0.5_wp * dx
1025!
1026!--             -1.0 for right-facing wall, 1.0 for left-facing wall
1027                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
1028                                  BTEST( wall_flags_0(k,j,i-1), 0 ) )
1029                dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
1030                dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
1031             ENDDO 
1032!
1033!--          Natural surfaces
1034             surf_s = surf_lsm_v(l)%start_index(j,i)
1035             surf_e = surf_lsm_v(l)%end_index(j,i)
1036             DO  m = surf_s, surf_e
1037                k           = surf_lsm_v(l)%k(m)
1038                vsus        = surf_lsm_v(l)%mom_flux_tke(0,m)
1039                wsus        = surf_lsm_v(l)%mom_flux_tke(1,m)
1040
1041                km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp            &
1042                                   * 0.5_wp * dx
1043!
1044!--             -1.0 for right-facing wall, 1.0 for left-facing wall
1045                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
1046                                  BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
1047                dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
1048                dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
1049             ENDDO 
1050!
1051!--          Urban surfaces
1052             surf_s = surf_usm_v(l)%start_index(j,i)
1053             surf_e = surf_usm_v(l)%end_index(j,i)
1054             DO  m = surf_s, surf_e
1055                k           = surf_usm_v(l)%k(m)
1056                vsus        = surf_usm_v(l)%mom_flux_tke(0,m)
1057                wsus        = surf_usm_v(l)%mom_flux_tke(1,m)
1058
1059                km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp            &
1060                                   * 0.5_wp * dx
1061!
1062!--             -1.0 for right-facing wall, 1.0 for left-facing wall
1063                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
1064                                  BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
1065                dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
1066                dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
1067             ENDDO 
1068          ENDDO
1069!
1070!--       Compute gradients at upward-facing walls, first for
1071!--       non-natural default surfaces
1072          surf_s = surf_def_h(0)%start_index(j,i)
1073          surf_e = surf_def_h(0)%end_index(j,i)
1074          DO  m = surf_s, surf_e
1075             k = surf_def_h(0)%k(m)
1076!
1077!--          Please note, actually, an interpolation of u_0 and v_0
1078!--          onto the grid center would be required. However, this
1079!--          would require several data transfers between 2D-grid and
1080!--          wall type. The effect of this missing interpolation is
1081!--          negligible. (See also production_e_init).
1082             dudz(k)     = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)   
1083             dvdz(k)     = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
1084
1085          ENDDO
1086!
1087!--       Natural surfaces
1088          surf_s = surf_lsm_h%start_index(j,i)
1089          surf_e = surf_lsm_h%end_index(j,i)
1090          DO  m = surf_s, surf_e
1091             k = surf_lsm_h%k(m)
1092
1093             dudz(k)     = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
1094             dvdz(k)     = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
1095          ENDDO
1096!
1097!--       Urban surfaces
1098          surf_s = surf_usm_h%start_index(j,i)
1099          surf_e = surf_usm_h%end_index(j,i)
1100          DO  m = surf_s, surf_e
1101             k = surf_usm_h%k(m)
1102
1103             dudz(k)     = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
1104             dvdz(k)     = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
1105          ENDDO
1106!
1107!--       Compute gradients at downward-facing walls, only for
1108!--       non-natural default surfaces
1109          surf_s = surf_def_h(1)%start_index(j,i)
1110          surf_e = surf_def_h(1)%end_index(j,i)
1111          DO  m = surf_s, surf_e
1112             k = surf_def_h(1)%k(m)
1113
1114             dudz(k)     = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
1115             dvdz(k)     = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 
1116
1117          ENDDO
1118
1119          DO  k = nzb+1, nzt
1120
1121             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
1122                              dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
1123                              dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
1124                   2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
1125
1126             IF ( def < 0.0_wp )  def = 0.0_wp
1127
1128             flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
1129
1130             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
1131
1132          ENDDO
1133
1134       ELSE
1135!
1136!--       Calculate TKE production by shear. Here, no additional
1137!--       wall-bounded code is considered.
1138!--       Why?
1139          DO  k = nzb+1, nzt
1140
1141             dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
1142             dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
1143                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1144             dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
1145                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1146
1147             dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
1148                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1149             dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
1150             dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
1151                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1152
1153             dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
1154                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1155             dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
1156                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1157             dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
1158
1159             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
1160                              dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
1161                              dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
1162                   2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
1163
1164             IF ( def < 0.0_wp )  def = 0.0_wp
1165
1166             flag        = MERGE( 1.0_wp, 0.0_wp,                              &
1167                                  BTEST( wall_flags_0(k,j,i), 29 ) ) 
1168             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
1169
1170          ENDDO
1171
1172       ENDIF
1173
1174!
1175!--    If required, calculate TKE production by buoyancy
1176       IF ( .NOT. neutral )  THEN
1177
1178          IF ( .NOT. humidity )  THEN
1179
1180             IF ( ocean )  THEN
1181!
1182!--             So far in the ocean no special treatment of density flux in
1183!--             the bottom and top surface layer
1184                DO  k = nzb+1, nzt
1185 
1186                   tend(k,j,i) = tend(k,j,i) +                                 &
1187                                 kh(k,j,i) * g /                               &
1188                                 MERGE( rho_reference, prho(k,j,i),            &
1189                                        use_single_reference_value ) *         &
1190                                 ( prho(k+1,j,i) - prho(k-1,j,i) ) *           &
1191                                 dd2zu(k) *                                    &
1192                                   MERGE( 1.0_wp, 0.0_wp,                      &
1193                                          BTEST( wall_flags_0(k,j,i), 30 )     &
1194                                        ) *                                    &
1195                                   MERGE( 1.0_wp, 0.0_wp,                      &
1196                                          BTEST( wall_flags_0(k,j,i), 9 )      &
1197                                        )           
1198                ENDDO
1199
1200                IF ( use_surface_fluxes )  THEN
1201!
1202!--                Default surfaces, up- and downward-facing
1203                   DO  l = 0, 1
1204                      surf_s = surf_def_h(l)%start_index(j,i)
1205                      surf_e = surf_def_h(l)%end_index(j,i)
1206                      DO  m = surf_s, surf_e
1207                         k = surf_def_h(l)%k(m)
1208                         tend(k,j,i) = tend(k,j,i) + g /                       &
1209                                   MERGE( rho_reference, prho(k,j,i),          &
1210                                          use_single_reference_value ) *       &
1211                                   drho_air_zw(k-1) *                          &
1212                                   surf_def_h(l)%shf(m)
1213                      ENDDO
1214                   ENDDO
1215
1216                ENDIF
1217
1218                IF ( use_top_fluxes )  THEN
1219                   surf_s = surf_def_h(2)%start_index(j,i)
1220                   surf_e = surf_def_h(2)%end_index(j,i)
1221                   DO  m = surf_s, surf_e
1222                      k = surf_def_h(2)%k(m)
1223                      tend(k,j,i) = tend(k,j,i) + g /                          &
1224                                   MERGE( rho_reference, prho(k,j,i),          &
1225                                          use_single_reference_value ) *       &
1226                                   drho_air_zw(k) *                            &
1227                                   surf_def_h(2)%shf(m) 
1228                   ENDDO
1229                ENDIF
1230
1231             ELSE
1232
1233                DO  k = nzb+1, nzt
1234!
1235!--                Flag 9 is used to mask top fluxes, flag 30 to mask
1236!--                surface fluxes
1237                   tend(k,j,i) = tend(k,j,i) -                                 &
1238                                 kh(k,j,i) * g /                               &
1239                                   MERGE( pt_reference, pt(k,j,i),             &
1240                                          use_single_reference_value ) *       &
1241                                 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) *    &
1242                                   MERGE( 1.0_wp, 0.0_wp,                      &
1243                                          BTEST( wall_flags_0(k,j,i), 30 )     &
1244                                        ) *                                    &
1245                                   MERGE( 1.0_wp, 0.0_wp,                      &
1246                                          BTEST( wall_flags_0(k,j,i), 9 )      &
1247                                        )
1248
1249                ENDDO
1250
1251                IF ( use_surface_fluxes )  THEN
1252!
1253!--                Default surfaces, up- and downward-facing
1254                   DO  l = 0, 1
1255                      surf_s = surf_def_h(l)%start_index(j,i)
1256                      surf_e = surf_def_h(l)%end_index(j,i)
1257                      DO  m = surf_s, surf_e
1258                         k = surf_def_h(l)%k(m)
1259                         tend(k,j,i) = tend(k,j,i) + g /                       &
1260                                   MERGE( pt_reference, pt(k,j,i),             &
1261                                          use_single_reference_value ) *       &
1262                                   drho_air_zw(k-1) *                          &
1263                                   surf_def_h(l)%shf(m)
1264                      ENDDO
1265                   ENDDO
1266!
1267!--                Natural surfaces
1268                   surf_s = surf_lsm_h%start_index(j,i)
1269                   surf_e = surf_lsm_h%end_index(j,i)
1270                   DO  m = surf_s, surf_e
1271                      k = surf_lsm_h%k(m)
1272                      tend(k,j,i) = tend(k,j,i) + g /                          &
1273                                   MERGE( pt_reference, pt(k,j,i),             &
1274                                          use_single_reference_value ) *       &
1275                                   drho_air_zw(k-1) *                          &
1276                                   surf_lsm_h%shf(m) 
1277                   ENDDO
1278!
1279!--                Urban surfaces
1280                   surf_s = surf_usm_h%start_index(j,i)
1281                   surf_e = surf_usm_h%end_index(j,i)
1282                   DO  m = surf_s, surf_e
1283                      k = surf_usm_h%k(m)
1284                      tend(k,j,i) = tend(k,j,i) + g /                          &
1285                                   MERGE( pt_reference, pt(k,j,i),             &
1286                                          use_single_reference_value ) *       &
1287                                   drho_air_zw(k-1) *                          &
1288                                   surf_usm_h%shf(m) 
1289                   ENDDO
1290                ENDIF
1291
1292                IF ( use_top_fluxes )  THEN
1293                   surf_s = surf_def_h(2)%start_index(j,i)
1294                   surf_e = surf_def_h(2)%end_index(j,i)
1295                   DO  m = surf_s, surf_e
1296                      k = surf_def_h(2)%k(m)
1297                      tend(k,j,i) = tend(k,j,i) + g /                          &
1298                                   MERGE( pt_reference, pt(k,j,i),             &
1299                                          use_single_reference_value ) *       &
1300                                   drho_air_zw(k) *                            &
1301                                   surf_def_h(2)%shf(m) 
1302                   ENDDO
1303                ENDIF
1304
1305             ENDIF
1306
1307          ELSE
1308
1309             DO  k = nzb+1, nzt
1310!
1311!--             Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes
1312                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
1313                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1314                   k2 = 0.61_wp * pt(k,j,i)
1315                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
1316                                   MERGE( vpt_reference, vpt(k,j,i),           &
1317                                          use_single_reference_value ) *       &
1318                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
1319                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
1320                                         ) * dd2zu(k) *                        &
1321                                      MERGE( 1.0_wp, 0.0_wp,                   &
1322                                             BTEST( wall_flags_0(k,j,i), 30 )  &
1323                                           )          *                        &
1324                                      MERGE( 1.0_wp, 0.0_wp,                   &
1325                                             BTEST( wall_flags_0(k,j,i), 9 )   &
1326                                           )
1327
1328                ELSE IF ( cloud_physics )  THEN
1329                   IF ( ql(k,j,i) == 0.0_wp )  THEN
1330                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1331                      k2 = 0.61_wp * pt(k,j,i)
1332                   ELSE
1333                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1334                      temp  = theta * t_d_pt(k)
1335                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                     &
1336                                    ( q(k,j,i) - ql(k,j,i) ) *                 &
1337                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /            &
1338                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *              &
1339                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1340                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1341                   ENDIF
1342                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
1343                                   MERGE( vpt_reference, vpt(k,j,i),           &
1344                                          use_single_reference_value ) *       &
1345                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
1346                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
1347                                         ) * dd2zu(k) *                        &
1348                                      MERGE( 1.0_wp, 0.0_wp,                   &
1349                                             BTEST( wall_flags_0(k,j,i), 30 )  &
1350                                           )          *                        &
1351                                      MERGE( 1.0_wp, 0.0_wp,                   &
1352                                             BTEST( wall_flags_0(k,j,i), 9 )   &
1353                                           )
1354                ELSE IF ( cloud_droplets )  THEN
1355                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1356                   k2 = 0.61_wp * pt(k,j,i)
1357                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
1358                                   MERGE( vpt_reference, vpt(k,j,i),           &
1359                                          use_single_reference_value ) *       &
1360                                     ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +      &
1361                                       k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -      &
1362                                       pt(k,j,i) * ( ql(k+1,j,i) -             &
1363                                                     ql(k-1,j,i) ) ) * dd2zu(k)&
1364                                    * MERGE( 1.0_wp, 0.0_wp,                   &
1365                                             BTEST( wall_flags_0(k,j,i), 30 )  &
1366                                           )                                   &
1367                                    * MERGE( 1.0_wp, 0.0_wp,                   &
1368                                             BTEST( wall_flags_0(k,j,i), 9 )   &
1369                                           )   
1370                ENDIF
1371             ENDDO
1372
1373             IF ( use_surface_fluxes )  THEN
1374!
1375!--             Treat horizontal default surfaces, up- and downward-facing
1376                DO  l = 0, 1
1377                   surf_s = surf_def_h(l)%start_index(j,i)
1378                   surf_e = surf_def_h(l)%end_index(j,i)
1379                   DO  m = surf_s, surf_e
1380                      k = surf_def_h(l)%k(m)
1381
1382                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1383                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1384                         k2 = 0.61_wp * pt(k,j,i)
1385                      ELSE IF ( cloud_physics )  THEN
1386                         IF ( ql(k,j,i) == 0.0_wp )  THEN
1387                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1388                            k2 = 0.61_wp * pt(k,j,i)
1389                         ELSE
1390                           theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1391                           temp  = theta * t_d_pt(k)
1392                           k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                &
1393                                      ( q(k,j,i) - ql(k,j,i) ) *               &
1394                             ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /          &
1395                             ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *            & 
1396                             ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1397                           k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1398                         ENDIF
1399                      ELSE IF ( cloud_droplets )  THEN
1400                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1401                         k2 = 0.61_wp * pt(k,j,i)
1402                      ENDIF
1403
1404                      tend(k,j,i) = tend(k,j,i) + g /                          &
1405                                   MERGE( vpt_reference, vpt(k,j,i),           &
1406                                          use_single_reference_value ) *       &
1407                                         ( k1 * surf_def_h(l)%shf(m) +         &
1408                                           k2 * surf_def_h(l)%qsws(m)          &
1409                                         ) * drho_air_zw(k-1)
1410                   ENDDO
1411                ENDDO
1412!
1413!--             Treat horizontal natural surfaces
1414                surf_s = surf_lsm_h%start_index(j,i)
1415                surf_e = surf_lsm_h%end_index(j,i)
1416                DO  m = surf_s, surf_e
1417                   k = surf_lsm_h%k(m)
1418
1419                   IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1420                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1421                       k2 = 0.61_wp * pt(k,j,i)
1422                   ELSE IF ( cloud_physics )  THEN
1423                       IF ( ql(k,j,i) == 0.0_wp )  THEN
1424                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1425                          k2 = 0.61_wp * pt(k,j,i)
1426                       ELSE
1427                          theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1428                          temp  = theta * t_d_pt(k)
1429                          k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
1430                                        ( q(k,j,i) - ql(k,j,i) ) *             &
1431                               ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
1432                               ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
1433                               ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1434                          k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1435                      ENDIF
1436                   ELSE IF ( cloud_droplets )  THEN
1437                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1438                      k2 = 0.61_wp * pt(k,j,i)
1439                   ENDIF
1440
1441                   tend(k,j,i) = tend(k,j,i) + g /                             &
1442                                   MERGE( vpt_reference, vpt(k,j,i),           &
1443                                          use_single_reference_value ) *       &
1444                                            ( k1 * surf_lsm_h%shf(m) +         &
1445                                              k2 * surf_lsm_h%qsws(m)          &
1446                                            ) * drho_air_zw(k-1)
1447                ENDDO
1448!
1449!--             Treat horizontal urban surfaces
1450                surf_s = surf_usm_h%start_index(j,i)
1451                surf_e = surf_usm_h%end_index(j,i)
1452                DO  m = surf_s, surf_e
1453                   k = surf_usm_h%k(m)
1454
1455                   IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1456                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1457                       k2 = 0.61_wp * pt(k,j,i)
1458                   ELSE IF ( cloud_physics )  THEN
1459                       IF ( ql(k,j,i) == 0.0_wp )  THEN
1460                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1461                          k2 = 0.61_wp * pt(k,j,i)
1462                       ELSE
1463                          theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1464                          temp  = theta * t_d_pt(k)
1465                          k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
1466                                        ( q(k,j,i) - ql(k,j,i) ) *             &
1467                               ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
1468                               ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
1469                               ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1470                          k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1471                      ENDIF
1472                   ELSE IF ( cloud_droplets )  THEN
1473                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1474                      k2 = 0.61_wp * pt(k,j,i)
1475                   ENDIF
1476
1477                   tend(k,j,i) = tend(k,j,i) + g /                             &
1478                                   MERGE( vpt_reference, vpt(k,j,i),           &
1479                                          use_single_reference_value ) *       &
1480                                            ( k1 * surf_usm_h%shf(m) +         &
1481                                              k2 * surf_usm_h%qsws(m)          &
1482                                            ) * drho_air_zw(k-1)
1483                ENDDO
1484
1485             ENDIF
1486
1487             IF ( use_top_fluxes )  THEN
1488                surf_s = surf_def_h(2)%start_index(j,i)
1489                surf_e = surf_def_h(2)%end_index(j,i)
1490                DO  m = surf_s, surf_e
1491                   k = surf_def_h(2)%k(m)
1492
1493
1494
1495                   IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
1496                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1497                      k2 = 0.61_wp * pt(k,j,i)
1498                   ELSE IF ( cloud_physics )  THEN
1499                      IF ( ql(k,j,i) == 0.0_wp )  THEN
1500                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1501                         k2 = 0.61_wp * pt(k,j,i)
1502                      ELSE
1503                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1504                         temp  = theta * t_d_pt(k)
1505                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
1506                                    ( q(k,j,i) - ql(k,j,i) ) *                 &
1507                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /            &
1508                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *              &
1509                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1510                         k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1511                      ENDIF
1512                   ELSE IF ( cloud_droplets )  THEN
1513                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1514                      k2 = 0.61_wp * pt(k,j,i)
1515                   ENDIF
1516
1517                   tend(k,j,i) = tend(k,j,i) + g /                             &
1518                                   MERGE( vpt_reference, vpt(k,j,i),           &
1519                                          use_single_reference_value ) *       &
1520                               ( k1* surf_def_h(2)%shf(m) +                    &
1521                                 k2 * surf_def_h(2)%qsws(m)                    &
1522                               ) * drho_air_zw(k)
1523                ENDDO
1524
1525             ENDIF
1526
1527          ENDIF
1528
1529       ENDIF
1530
1531    END SUBROUTINE production_e_ij
1532
1533
1534!------------------------------------------------------------------------------!
1535! Description:
1536! ------------
1537!> @todo Missing subroutine description.
1538!------------------------------------------------------------------------------!
1539    SUBROUTINE production_e_init
1540
1541       USE arrays_3d,                                                          &
1542           ONLY:  kh, km, drho_air_zw, u, v, zu
1543
1544       USE control_parameters,                                                 &
1545           ONLY:  constant_flux_layer, kappa
1546
1547       USE indices,                                                            &
1548           ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_u_inner,     &
1549                  nzb_v_inner
1550
1551       USE surface_mod,                                                        &
1552           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h 
1553
1554       IMPLICIT NONE
1555
1556       INTEGER(iwp) ::  i   !< grid index x-direction
1557       INTEGER(iwp) ::  j   !< grid index y-direction
1558       INTEGER(iwp) ::  k   !< grid index z-direction
1559       INTEGER(iwp) ::  l   !< running index surface type (up- or downward-facing)
1560       INTEGER(iwp) ::  m   !< running index surface elements
1561
1562       IF ( constant_flux_layer )  THEN
1563!
1564!--       Calculate a virtual velocity at the surface in a way that the
1565!--       vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
1566!--       Prandtl law (-w'u'/km). This gradient is used in the TKE shear
1567!--       production term at k=1 (see production_e_ij).
1568!--       The velocity gradient has to be limited in case of too small km
1569!--       (otherwise the timestep may be significantly reduced by large
1570!--       surface winds).
1571!--       not available in case of non-cyclic boundary conditions.
1572!--       WARNING: the exact analytical solution would require the determination
1573!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
1574!--       Default surfaces, upward-facing
1575          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
1576          DO  m = 1, surf_def_h(0)%ns
1577
1578             i = surf_def_h(0)%i(m)           
1579             j = surf_def_h(0)%j(m)
1580             k = surf_def_h(0)%k(m)
1581!
1582!--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
1583!--          and km are not on the same grid. Actually, a further
1584!--          interpolation of km onto the u/v-grid is necessary. However, the
1585!--          effect of this error is negligible.
1586             surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *       &
1587                                        drho_air_zw(k-1) *                     &
1588                                        ( zu(k+1)    - zu(k-1)    )  /         &
1589                                        ( km(k,j,i)  + 1.0E-20_wp ) 
1590             surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *       &
1591                                        drho_air_zw(k-1) *                     &
1592                                        ( zu(k+1)    - zu(k-1)    )  /         &
1593                                        ( km(k,j,i)  + 1.0E-20_wp ) 
1594
1595             IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                  &
1596                  ABS( u(k+1,j,i) - u(k-1,j,i)           )                     &
1597                )  surf_def_h(0)%u_0(m) = u(k-1,j,i)
1598
1599             IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) )  >                  &
1600                  ABS( v(k+1,j,i) - v(k-1,j,i)           )                     &
1601                )  surf_def_h(0)%v_0(m) = v(k-1,j,i)
1602
1603          ENDDO
1604!
1605!--       Default surfaces, downward-facing surfaces
1606          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
1607          DO  m = 1, surf_def_h(1)%ns
1608
1609             i = surf_def_h(1)%i(m)           
1610             j = surf_def_h(1)%j(m)
1611             k = surf_def_h(1)%k(m)
1612
1613             surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *    &
1614                                        drho_air_zw(k-1) *                  &
1615                                        ( zu(k+1)    - zu(k-1)    )  /      &
1616                                        ( km(k,j,i)  + 1.0E-20_wp ) 
1617             surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *    &
1618                                        drho_air_zw(k-1) *                  &
1619                                        ( zu(k+1)    - zu(k-1)    )  /      &
1620                                        ( km(k,j,i)  + 1.0E-20_wp ) 
1621
1622             IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >               &
1623                  ABS( u(k+1,j,i)           - u(k-1,j,i) )                  &
1624                )  surf_def_h(1)%u_0(m) = u(k+1,j,i)
1625
1626             IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >               &
1627                  ABS( v(k+1,j,i)           - v(k-1,j,i) )                  &
1628                )  surf_def_h(1)%v_0(m) = v(k+1,j,i)
1629
1630          ENDDO
1631!
1632!--       Natural surfaces, upward-facing
1633          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
1634          DO  m = 1, surf_lsm_h%ns
1635
1636             i = surf_lsm_h%i(m)           
1637             j = surf_lsm_h%j(m)
1638             k = surf_lsm_h%k(m)
1639!
1640!--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
1641!--          and km are not on the same grid. Actually, a further
1642!--          interpolation of km onto the u/v-grid is necessary. However, the
1643!--          effect of this error is negligible.
1644             surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)      *        &
1645                                           drho_air_zw(k-1) *                  &
1646                                           ( zu(k+1)   - zu(k-1)    )  /       &
1647                                           ( km(k,j,i) + 1.0E-20_wp ) 
1648             surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)      *        &
1649                                           drho_air_zw(k-1) *                  &
1650                                           ( zu(k+1)   - zu(k-1)    )  /       &
1651                                           ( km(k,j,i) + 1.0E-20_wp )
1652
1653             IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                     &
1654                  ABS( u(k+1,j,i) - u(k-1,j,i)   )                             &
1655                )  surf_lsm_h%u_0(m) = u(k-1,j,i)
1656
1657             IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) )  >                     &
1658                  ABS( v(k+1,j,i) - v(k-1,j,i)   )                             &
1659                )  surf_lsm_h%v_0(m) = v(k-1,j,i)
1660
1661          ENDDO
1662!
1663!--       Urban surfaces, upward-facing
1664          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
1665          DO  m = 1, surf_usm_h%ns
1666
1667             i = surf_usm_h%i(m)           
1668             j = surf_usm_h%j(m)
1669             k = surf_usm_h%k(m)
1670!
1671!--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
1672!--          and km are not on the same grid. Actually, a further
1673!--          interpolation of km onto the u/v-grid is necessary. However, the
1674!--          effect of this error is negligible.
1675             surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)      *        &
1676                                           drho_air_zw(k-1) *                  &
1677                                           ( zu(k+1)   - zu(k-1)    )  /       &
1678                                           ( km(k,j,i) + 1.0E-20_wp ) 
1679             surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)      *        &
1680                                           drho_air_zw(k-1) *                  &
1681                                           ( zu(k+1)   - zu(k-1)    )  /       &
1682                                           ( km(k,j,i) + 1.0E-20_wp )
1683
1684             IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                     &
1685                  ABS( u(k+1,j,i) - u(k-1,j,i)   )                             &
1686                )  surf_usm_h%u_0(m) = u(k-1,j,i)
1687
1688             IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) )  >                     &
1689                  ABS( v(k+1,j,i) - v(k-1,j,i)   )                             &
1690                )  surf_usm_h%v_0(m) = v(k-1,j,i)
1691
1692          ENDDO
1693
1694       ENDIF
1695
1696    END SUBROUTINE production_e_init
1697
1698 END MODULE production_e_mod
Note: See TracBrowser for help on using the repository browser.