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

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

revised renaming of modules

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