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

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

added _mod string to several filenames to meet the naming convection for modules

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