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

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

last commit documented

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