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

Last change on this file since 2416 was 2329, checked in by knoop, 7 years ago

Bugfix for topography usage with anelastic approximation and boussinesq approximation with air density != 1

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