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

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

last commit documented

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