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

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

last commit documented

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