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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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