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

Last change on this file since 2492 was 2478, checked in by suehring, 7 years ago

Bugfixes concerning top fluxes and TKE production

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