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

Last change on this file since 1165 was 1132, checked in by raasch, 12 years ago

r1028 documented

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