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

Last change on this file since 2113 was 2101, checked in by suehring, 7 years ago

last commit documented

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