source: palm/tags/release-3.10/SOURCE/production_e.f90 @ 1614

Last change on this file since 1614 was 1258, checked in by raasch, 10 years ago

last commit documented

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