source: palm/tags/release-4.0/SOURCE/production_e.f90 @ 1613

Last change on this file since 1613 was 1375, checked in by raasch, 10 years ago

last commit documented

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