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

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

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

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