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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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