source: palm/trunk/SOURCE/production_e_mod.f90 @ 1859

Last change on this file since 1859 was 1851, checked in by maronga, 9 years ago

last commit documented

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