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

Last change on this file since 1874 was 1874, checked in by maronga, 5 years ago

last commit documented

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