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

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

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

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