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

Last change on this file since 2008 was 2001, checked in by knoop, 8 years ago

last commit documented

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