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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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