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

Last change on this file since 2053 was 2032, checked in by knoop, 8 years ago

last commit documented

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