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

Last change on this file since 1874 was 1874, checked in by maronga, 8 years ago

last commit documented

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