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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

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