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

Last change on this file since 1052 was 1037, checked in by raasch, 12 years ago

last commit documented

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