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

Last change on this file since 1207 was 1182, checked in by raasch, 11 years ago

last commit documented, rc-file for example run updated

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