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

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

last commit documented

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