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

Last change on this file since 213 was 208, checked in by raasch, 15 years ago

bugfix in production_e

  • Property svn:keywords set to Id
File size: 45.6 KB
Line 
1 MODULE production_e_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Bugfix concerning the calculation of velocity gradients at vertical walls
7! in case of diabatic conditions
8!
9! Former revisions:
10! -----------------
11! $Id: production_e.f90 208 2008-10-20 06:02:59Z raasch $
12!
13! 187 2008-08-06 16:25:09Z letzel
14! Change: add 'minus' sign to fluxes obtained from subroutine wall_fluxes_e for
15! consistency with subroutine wall_fluxes
16!
17! 124 2007-10-19 15:47:46Z raasch
18! Bugfix: calculation of density flux in the ocean now starts from nzb+1
19!
20! 108 2007-08-24 15:10:38Z letzel
21! Bugfix: wrong sign removed from the buoyancy production term in the case
22! use_reference = .T.,
23! u_0 and v_0 are calculated for nxr+1, nyn+1 also (otherwise these values are
24! not available in case of non-cyclic boundary conditions)
25! Bugfix for ocean density flux at bottom
26!
27! 97 2007-06-21 08:23:15Z raasch
28! Energy production by density flux (in ocean) added
29! use_pt_reference renamed use_reference
30!
31! 75 2007-03-22 09:54:05Z raasch
32! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
33! reference temperature pt_reference can be used in buoyancy term,
34! moisture renamed humidity
35!
36! 37 2007-03-01 08:33:54Z raasch
37! Calculation extended for gridpoint nzt, extended for given temperature /
38! humidity fluxes at the top, wall-part is now executed in case that a
39! Prandtl-layer is switched on (instead of surfaces fluxes switched on)
40!
41! RCS Log replace by Id keyword, revision history cleaned up
42!
43! Revision 1.21  2006/04/26 12:45:35  raasch
44! OpenMP parallelization of production_e_init
45!
46! Revision 1.1  1997/09/19 07:45:35  raasch
47! Initial revision
48!
49!
50! Description:
51! ------------
52! Production terms (shear + buoyancy) of the TKE
53! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
54!          not considered well!
55!------------------------------------------------------------------------------!
56
57    USE wall_fluxes_mod
58
59    PRIVATE
60    PUBLIC production_e, production_e_init
61
62    LOGICAL, SAVE ::  first_call = .TRUE.
63
64    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  u_0, v_0
65
66    INTERFACE production_e
67       MODULE PROCEDURE production_e
68       MODULE PROCEDURE production_e_ij
69    END INTERFACE production_e
70   
71    INTERFACE production_e_init
72       MODULE PROCEDURE production_e_init
73    END INTERFACE production_e_init
74 
75 CONTAINS
76
77
78!------------------------------------------------------------------------------!
79! Call for all grid points
80!------------------------------------------------------------------------------!
81    SUBROUTINE production_e
82
83       USE arrays_3d
84       USE cloud_parameters
85       USE control_parameters
86       USE grid_variables
87       USE indices
88       USE statistics
89
90       IMPLICIT NONE
91
92       INTEGER ::  i, j, k
93
94       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
95                   k1, k2, km_neutral, theta, temp
96
97!       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
98       REAL, DIMENSION(nzb:nzt+1) ::   usvs, vsus, wsus, wsvs
99
100!
101!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
102!--    vertical walls, if neccessary
103!--    So far, results are slightly different from the ij-Version.
104!--    Therefore, ij-Version is called further below within the ij-loops.
105!       IF ( topography /= 'flat' )  THEN
106!          CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
107!          CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
108!          CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
109!          CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
110!       ENDIF
111
112!
113!--    Calculate TKE production by shear
114       DO  i = nxl, nxr
115
116          DO  j = nys, nyn
117             DO  k = nzb_diff_s_outer(j,i), nzt
118
119                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
120                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
121                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
122                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
123                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
124
125                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
126                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
127                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
128                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
129                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
130
131                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
132                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
133                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
134                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
135                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
136
137                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
138                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
139                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
140
141                IF ( def < 0.0 )  def = 0.0
142
143                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
144               
145             ENDDO
146          ENDDO
147
148          IF ( prandtl_layer )  THEN
149
150!
151!--          Position beneath wall
152!--          (2) - Will allways be executed.
153!--          'bottom and wall: use u_0,v_0 and wall functions'
154             DO  j = nys, nyn
155
156                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
157                THEN
158
159                   k = nzb_diff_s_inner(j,i) - 1
160                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
161                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
162                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
163                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
164                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
165                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
166                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
167
168                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
169!                     
170!--                   Inconsistency removed: as the thermal stratification is
171!--                   not taken into account for the evaluation of the wall
172!--                   fluxes at vertical walls, the eddy viscosity km must not
173!--                   be used for the evaluation of the velocity gradients dudy
174!--                   and dwdy
175!--                   Note: The validity of the new method has not yet been
176!--                         shown, as so far no suitable data for a validation
177!--                         has been available
178                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
179                                          usvs, 1.0, 0.0, 0.0, 0.0 )
180                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
181                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
182                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
183                                   0.5 * dy 
184                      dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
185                      dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
186                   ELSE
187                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
188                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
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                   ENDIF
192
193                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
194!                     
195!--                   Inconsistency removed: as the thermal stratification is
196!--                   not taken into account for the evaluation of the wall
197!--                   fluxes at vertical walls, the eddy viscosity km must not
198!--                   be used for the evaluation of the velocity gradients dvdx
199!--                   and dwdx
200!--                   Note: The validity of the new method has not yet been
201!--                         shown, as so far no suitable data for a validation
202!--                         has been available
203                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
204                                          vsus, 0.0, 1.0, 0.0, 0.0 )
205                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
206                                          wsus, 0.0, 0.0, 0.0, 1.0 )
207                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &
208                                   0.5 * dx
209                      dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
210                      dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
211                   ELSE
212                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
213                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
214                      dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
215                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
216                   ENDIF
217
218                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
219                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
220                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
221
222                   IF ( def < 0.0 )  def = 0.0
223
224                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
225
226
227!
228!--                (3) - will be executed only, if there is at least one level
229!--                between (2) and (4), i.e. the topography must have a
230!--                minimum height of 2 dz. Wall fluxes for this case have
231!--                already been calculated for (2).
232!--                'wall only: use wall functions'
233
234                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
235
236                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
237                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
238                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
239                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
240                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
241                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
242                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
243
244                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
245!                     
246!--                      Inconsistency removed: as the thermal stratification
247!--                      is not taken into account for the evaluation of the
248!--                      wall fluxes at vertical walls, the eddy viscosity km
249!--                      must not be used for the evaluation of the velocity
250!--                      gradients dudy and dwdy
251!--                      Note: The validity of the new method has not yet
252!--                            been shown, as so far no suitable data for a
253!--                            validation has been available
254                         km_neutral = kappa * ( usvs(k)**2 + & 
255                                                wsvs(k)**2 )**0.25 * 0.5 * dy
256                         dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
257                         dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
258                      ELSE
259                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
260                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
261                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
262                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
263                      ENDIF
264
265                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
266!                     
267!--                      Inconsistency removed: as the thermal stratification
268!--                      is not taken into account for the evaluation of the
269!--                      wall fluxes at vertical walls, the eddy viscosity km
270!--                      must not be used for the evaluation of the velocity
271!--                      gradients dvdx and dwdx
272!--                      Note: The validity of the new method has not yet
273!--                            been shown, as so far no suitable data for a
274!--                            validation has been available
275                         km_neutral = kappa * ( vsus(k)**2 + & 
276                                                wsus(k)**2 )**0.25 * 0.5 * dx
277                         dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
278                         dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
279                      ELSE
280                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
281                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
282                         dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
283                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
284                      ENDIF
285
286                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
287                           dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
288                           dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
289
290                      IF ( def < 0.0 )  def = 0.0
291
292                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
293
294                   ENDDO
295
296                ENDIF
297
298             ENDDO
299
300!
301!--          (4) - will allways be executed.
302!--          'special case: free atmosphere' (as for case (0))
303             DO  j = nys, nyn
304
305                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
306                THEN
307
308                   k = nzb_diff_s_outer(j,i)-1
309
310                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
311                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
312                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
313                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
314                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
315
316                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
317                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
318                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
319                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
320                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
321
322                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
323                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
324                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
325                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
326                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
327
328                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
329                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
330                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
331
332                   IF ( def < 0.0 )  def = 0.0
333
334                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
335
336                ENDIF
337
338             ENDDO
339
340!
341!--          Position without adjacent wall
342!--          (1) - will allways be executed.
343!--          'bottom only: use u_0,v_0'
344             DO  j = nys, nyn
345
346                IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
347                THEN
348
349                   k = nzb_diff_s_inner(j,i)-1
350
351                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
352                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
353                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
354                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
355                                    u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
356
357                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
358                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
359                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
360                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
361                                    v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
362
363                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
364                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
365                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
366                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
367                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
368
369                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
370                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
371                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
372
373                   IF ( def < 0.0 )  def = 0.0
374
375                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
376               
377                ENDIF
378
379             ENDDO
380
381          ELSEIF ( use_surface_fluxes )  THEN
382
383             DO  j = nys, nyn
384
385                k = nzb_diff_s_outer(j,i)-1
386
387                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
388                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
389                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
390                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
391                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
392
393                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
394                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
395                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
396                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
397                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
398
399                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
400                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
401                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
402                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
403                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
404
405                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
406                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
407                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
408
409                IF ( def < 0.0 )  def = 0.0
410
411                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
412
413             ENDDO
414
415          ENDIF
416
417!
418!--       Calculate TKE production by buoyancy
419          IF ( .NOT. humidity )  THEN
420
421             IF ( use_reference )  THEN
422
423                IF ( ocean )  THEN
424!
425!--                So far in the ocean no special treatment of density flux in
426!--                the bottom and top surface layer
427                   DO  j = nys, nyn
428                      DO  k = nzb_s_inner(j,i)+1, nzt
429                         tend(k,j,i) = tend(k,j,i) +                    &
430                                       kh(k,j,i) * g / prho_reference * &
431                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
432                      ENDDO
433                   ENDDO
434
435                ELSE
436
437                   DO  j = nys, nyn
438                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
439                         tend(k,j,i) = tend(k,j,i) -                  &
440                                       kh(k,j,i) * g / pt_reference * &
441                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
442                      ENDDO
443
444                      IF ( use_surface_fluxes )  THEN
445                         k = nzb_diff_s_inner(j,i)-1
446                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
447                      ENDIF
448
449                      IF ( use_top_fluxes )  THEN
450                         k = nzt
451                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
452                                                     tswst(j,i)
453                      ENDIF
454                   ENDDO
455
456                ENDIF
457
458             ELSE
459
460                IF ( ocean )  THEN
461!
462!--                So far in the ocean no special treatment of density flux in
463!--                the bottom and top surface layer
464                   DO  j = nys, nyn
465                      DO  k = nzb_s_inner(j,i)+1, nzt
466                         tend(k,j,i) = tend(k,j,i) -                &
467                                       kh(k,j,i) * g / rho(k,j,i) * &
468                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
469                      ENDDO
470                   ENDDO
471
472                ELSE
473
474                   DO  j = nys, nyn
475                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
476                         tend(k,j,i) = tend(k,j,i) -               &
477                                       kh(k,j,i) * g / pt(k,j,i) * &
478                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
479                      ENDDO
480
481                      IF ( use_surface_fluxes )  THEN
482                         k = nzb_diff_s_inner(j,i)-1
483                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
484                      ENDIF
485
486                      IF ( use_top_fluxes )  THEN
487                         k = nzt
488                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
489                      ENDIF
490                   ENDDO
491
492                ENDIF
493
494             ENDIF
495
496          ELSE
497
498             DO  j = nys, nyn
499
500                DO  k = nzb_diff_s_inner(j,i), nzt_diff
501
502                   IF ( .NOT. cloud_physics )  THEN
503                      k1 = 1.0 + 0.61 * q(k,j,i)
504                      k2 = 0.61 * pt(k,j,i)
505                   ELSE
506                      IF ( ql(k,j,i) == 0.0 )  THEN
507                         k1 = 1.0 + 0.61 * q(k,j,i)
508                         k2 = 0.61 * pt(k,j,i)
509                      ELSE
510                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
511                         temp  = theta * t_d_pt(k)
512                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
513                                    ( q(k,j,i) - ql(k,j,i) ) *          &
514                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
515                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
516                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
517                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
518                      ENDIF
519                   ENDIF
520
521                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
522                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
523                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
524                                      ) * dd2zu(k)
525                ENDDO
526
527             ENDDO
528
529             IF ( use_surface_fluxes )  THEN
530
531                DO  j = nys, nyn
532
533                   k = nzb_diff_s_inner(j,i)-1
534
535                   IF ( .NOT. cloud_physics )  THEN
536                      k1 = 1.0 + 0.61 * q(k,j,i)
537                      k2 = 0.61 * pt(k,j,i)
538                   ELSE
539                      IF ( ql(k,j,i) == 0.0 )  THEN
540                         k1 = 1.0 + 0.61 * q(k,j,i)
541                         k2 = 0.61 * pt(k,j,i)
542                      ELSE
543                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
544                         temp  = theta * t_d_pt(k)
545                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
546                                    ( q(k,j,i) - ql(k,j,i) ) *          &
547                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
548                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
549                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
550                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
551                      ENDIF
552                   ENDIF
553
554                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
555                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
556                ENDDO
557
558             ENDIF
559
560             IF ( use_top_fluxes )  THEN
561
562                DO  j = nys, nyn
563
564                   k = nzt
565
566                   IF ( .NOT. cloud_physics )  THEN
567                      k1 = 1.0 + 0.61 * q(k,j,i)
568                      k2 = 0.61 * pt(k,j,i)
569                   ELSE
570                      IF ( ql(k,j,i) == 0.0 )  THEN
571                         k1 = 1.0 + 0.61 * q(k,j,i)
572                         k2 = 0.61 * pt(k,j,i)
573                      ELSE
574                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
575                         temp  = theta * t_d_pt(k)
576                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
577                                    ( q(k,j,i) - ql(k,j,i) ) *          &
578                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
579                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
580                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
581                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
582                      ENDIF
583                   ENDIF
584
585                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
586                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
587                ENDDO
588
589             ENDIF
590
591          ENDIF
592
593       ENDDO
594
595    END SUBROUTINE production_e
596
597
598!------------------------------------------------------------------------------!
599! Call for grid point i,j
600!------------------------------------------------------------------------------!
601    SUBROUTINE production_e_ij( i, j )
602
603       USE arrays_3d
604       USE cloud_parameters
605       USE control_parameters
606       USE grid_variables
607       USE indices
608       USE statistics
609
610       IMPLICIT NONE
611
612       INTEGER ::  i, j, k
613
614       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
615                   k1, k2, km_neutral, theta, temp
616
617       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs
618
619!
620!--    Calculate TKE production by shear
621       DO  k = nzb_diff_s_outer(j,i), nzt
622
623          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
624          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
625                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
626          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
627                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
628
629          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
630                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
631          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
632          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
633                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
634
635          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
636                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
637          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
638                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
639          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
640
641          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
642                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
643                + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
644
645          IF ( def < 0.0 )  def = 0.0
646
647          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
648               
649       ENDDO
650
651       IF ( prandtl_layer )  THEN
652
653          IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
654
655!
656!--          Position beneath wall
657!--          (2) - Will allways be executed.
658!--          'bottom and wall: use u_0,v_0 and wall functions'
659             k = nzb_diff_s_inner(j,i)-1
660
661             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
662             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
663                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
664             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
665             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
666                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
667             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
668
669             IF ( wall_e_y(j,i) /= 0.0 )  THEN
670!                     
671!--             Inconsistency removed: as the thermal stratification
672!--             is not taken into account for the evaluation of the
673!--             wall fluxes at vertical walls, the eddy viscosity km
674!--             must not be used for the evaluation of the velocity
675!--             gradients dudy and dwdy
676!--             Note: The validity of the new method has not yet
677!--                   been shown, as so far no suitable data for a
678!--                   validation has been available
679                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
680                                    usvs, 1.0, 0.0, 0.0, 0.0 )
681                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
682                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
683                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
684                             0.5 * dy
685                dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
686                dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
687             ELSE
688                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
689                                u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
690                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
691                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
692             ENDIF
693
694             IF ( wall_e_x(j,i) /= 0.0 )  THEN
695!                     
696!--             Inconsistency removed: as the thermal stratification
697!--             is not taken into account for the evaluation of the
698!--             wall fluxes at vertical walls, the eddy viscosity km
699!--             must not be used for the evaluation of the velocity
700!--             gradients dvdx and dwdx
701!--             Note: The validity of the new method has not yet
702!--                   been shown, as so far no suitable data for a
703!--                   validation has been available
704                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
705                                    vsus, 0.0, 1.0, 0.0, 0.0 )
706                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
707                                    wsus, 0.0, 0.0, 0.0, 1.0 )
708                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * & 
709                             0.5 * dx
710                dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
711                dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
712             ELSE
713                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
714                                v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
715                dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
716                                w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
717             ENDIF
718
719             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
720                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
721                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
722
723             IF ( def < 0.0 )  def = 0.0
724
725             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
726
727!
728!--          (3) - will be executed only, if there is at least one level
729!--          between (2) and (4), i.e. the topography must have a
730!--          minimum height of 2 dz. Wall fluxes for this case have
731!--          already been calculated for (2).
732!--          'wall only: use wall functions'
733             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
734
735                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
736                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
737                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
738                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
739                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
740                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
741                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
742
743                IF ( wall_e_y(j,i) /= 0.0 )  THEN
744!                     
745!--                Inconsistency removed: as the thermal stratification
746!--                is not taken into account for the evaluation of the
747!--                wall fluxes at vertical walls, the eddy viscosity km
748!--                must not be used for the evaluation of the velocity
749!--                gradients dudy and dwdy
750!--                Note: The validity of the new method has not yet
751!--                      been shown, as so far no suitable data for a
752!--                      validation has been available
753                   km_neutral = kappa * ( usvs(k)**2 + & 
754                                          wsvs(k)**2 )**0.25 * 0.5 * dy
755                   dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
756                   dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
757                ELSE
758                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
759                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
760                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
761                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
762                ENDIF
763
764                IF ( wall_e_x(j,i) /= 0.0 )  THEN
765!                     
766!--                Inconsistency removed: as the thermal stratification
767!--                is not taken into account for the evaluation of the
768!--                wall fluxes at vertical walls, the eddy viscosity km
769!--                must not be used for the evaluation of the velocity
770!--                gradients dvdx and dwdx
771!--                Note: The validity of the new method has not yet
772!--                      been shown, as so far no suitable data for a
773!--                      validation has been available
774                   km_neutral = kappa * ( vsus(k)**2 + & 
775                                          wsus(k)**2 )**0.25 * 0.5 * dx
776                   dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
777                   dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
778                ELSE
779                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
780                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
781                   dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
782                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
783                ENDIF
784
785                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
786                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
787                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
788
789                IF ( def < 0.0 )  def = 0.0
790
791                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
792
793             ENDDO
794
795!
796!--          (4) - will allways be executed.
797!--          'special case: free atmosphere' (as for case (0))
798             k = nzb_diff_s_outer(j,i)-1
799
800             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
801             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
802                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
803             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
804                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
805
806             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
807                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
808             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
809             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
810                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
811
812             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
813                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
814             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
815                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
816             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
817
818             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
819                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
820                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
821
822             IF ( def < 0.0 )  def = 0.0
823
824             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
825
826          ELSE
827
828!
829!--          Position without adjacent wall
830!--          (1) - will allways be executed.
831!--          'bottom only: use u_0,v_0'
832             k = nzb_diff_s_inner(j,i)-1
833
834             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
835             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
836                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
837             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
838                              u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
839
840             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
841                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
842             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
843             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
844                              v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
845
846             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
847                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
848             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
849                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
850             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
851
852             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
853                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
854                   + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
855
856             IF ( def < 0.0 )  def = 0.0
857
858             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
859
860          ENDIF
861
862       ELSEIF ( use_surface_fluxes )  THEN
863
864          k = nzb_diff_s_outer(j,i)-1
865
866          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
867          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
868                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
869          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
870                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
871
872          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
873                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
874          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
875          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
876                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
877
878          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
879                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
880          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
881                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
882          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
883
884          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
885                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
886                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
887
888          IF ( def < 0.0 )  def = 0.0
889
890          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
891
892       ENDIF
893
894!
895!--    Calculate TKE production by buoyancy
896       IF ( .NOT. humidity )  THEN
897
898          IF ( use_reference )  THEN
899
900             IF ( ocean )  THEN
901!
902!--             So far in the ocean no special treatment of density flux in the
903!--             bottom and top surface layer
904                DO  k = nzb_s_inner(j,i)+1, nzt
905                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
906                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
907                ENDDO
908
909             ELSE
910
911                DO  k = nzb_diff_s_inner(j,i), nzt_diff
912                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
913                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
914                ENDDO
915
916                IF ( use_surface_fluxes )  THEN
917                   k = nzb_diff_s_inner(j,i)-1
918                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
919                ENDIF
920
921                IF ( use_top_fluxes )  THEN
922                   k = nzt
923                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
924                ENDIF
925
926             ENDIF
927
928          ELSE
929
930             IF ( ocean )  THEN
931!
932!--             So far in the ocean no special treatment of density flux in the
933!--             bottom and top surface layer
934                DO  k = nzb_s_inner(j,i)+1, nzt
935                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho(k,j,i) * &
936                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
937                ENDDO
938
939             ELSE
940
941                DO  k = nzb_diff_s_inner(j,i), nzt_diff
942                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
943                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
944                ENDDO
945
946                IF ( use_surface_fluxes )  THEN
947                   k = nzb_diff_s_inner(j,i)-1
948                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
949                ENDIF
950
951                IF ( use_top_fluxes )  THEN
952                   k = nzt
953                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
954                ENDIF
955
956             ENDIF
957
958          ENDIF
959
960       ELSE
961
962          DO  k = nzb_diff_s_inner(j,i), nzt_diff
963
964             IF ( .NOT. cloud_physics )  THEN
965                k1 = 1.0 + 0.61 * q(k,j,i)
966                k2 = 0.61 * pt(k,j,i)
967             ELSE
968                IF ( ql(k,j,i) == 0.0 )  THEN
969                   k1 = 1.0 + 0.61 * q(k,j,i)
970                   k2 = 0.61 * pt(k,j,i)
971                ELSE
972                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
973                   temp  = theta * t_d_pt(k)
974                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
975                              ( q(k,j,i) - ql(k,j,i) ) *          &
976                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
977                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
978                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
979                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
980                ENDIF
981             ENDIF
982
983             tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *      &
984                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
985                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
986                                      ) * dd2zu(k)
987          ENDDO
988
989          IF ( use_surface_fluxes )  THEN
990             k = nzb_diff_s_inner(j,i)-1
991
992             IF ( .NOT. cloud_physics ) THEN
993                k1 = 1.0 + 0.61 * q(k,j,i)
994                k2 = 0.61 * pt(k,j,i)
995             ELSE
996                IF ( ql(k,j,i) == 0.0 )  THEN
997                   k1 = 1.0 + 0.61 * q(k,j,i)
998                   k2 = 0.61 * pt(k,j,i)
999                ELSE
1000                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1001                   temp  = theta * t_d_pt(k)
1002                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
1003                              ( q(k,j,i) - ql(k,j,i) ) *          &
1004                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
1005                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
1006                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1007                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
1008                ENDIF
1009             ENDIF
1010
1011             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1012                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
1013          ENDIF
1014
1015          IF ( use_top_fluxes )  THEN
1016             k = nzt
1017
1018             IF ( .NOT. cloud_physics ) THEN
1019                k1 = 1.0 + 0.61 * q(k,j,i)
1020                k2 = 0.61 * pt(k,j,i)
1021             ELSE
1022                IF ( ql(k,j,i) == 0.0 )  THEN
1023                   k1 = 1.0 + 0.61 * q(k,j,i)
1024                   k2 = 0.61 * pt(k,j,i)
1025                ELSE
1026                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1027                   temp  = theta * t_d_pt(k)
1028                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
1029                              ( q(k,j,i) - ql(k,j,i) ) *          &
1030                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
1031                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
1032                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1033                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
1034                ENDIF
1035             ENDIF
1036
1037             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1038                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
1039          ENDIF
1040
1041       ENDIF
1042
1043    END SUBROUTINE production_e_ij
1044
1045
1046    SUBROUTINE production_e_init
1047
1048       USE arrays_3d
1049       USE control_parameters
1050       USE grid_variables
1051       USE indices
1052
1053       IMPLICIT NONE
1054
1055       INTEGER ::  i, j, ku, kv
1056
1057       IF ( prandtl_layer )  THEN
1058
1059          IF ( first_call )  THEN
1060             ALLOCATE( u_0(nys-1:nyn+1,nxl-1:nxr+1), &
1061                       v_0(nys-1:nyn+1,nxl-1:nxr+1) )
1062             first_call = .FALSE.
1063          ENDIF
1064
1065!
1066!--       Calculate a virtual velocity at the surface in a way that the
1067!--       vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
1068!--       Prandtl law (-w'u'/km). This gradient is used in the TKE shear
1069!--       production term at k=1 (see production_e_ij).
1070!--       The velocity gradient has to be limited in case of too small km
1071!--       (otherwise the timestep may be significantly reduced by large
1072!--       surface winds).
1073!--       Upper bounds are nxr+1 and nyn+1 because otherwise these values are
1074!--       not available in case of non-cyclic boundary conditions.
1075!--       WARNING: the exact analytical solution would require the determination
1076!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
1077          !$OMP PARALLEL DO PRIVATE( ku, kv )
1078          DO  i = nxl, nxr+1
1079             DO  j = nys, nyn+1
1080
1081                ku = nzb_u_inner(j,i)+1
1082                kv = nzb_v_inner(j,i)+1
1083
1084                u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
1085                                 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +       &
1086                                   1.0E-20 )
1087!                                  ( us(j,i) * kappa * zu(1) )
1088                v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
1089                                 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +       &
1090                                   1.0E-20 )
1091!                                  ( us(j,i) * kappa * zu(1) )
1092
1093                IF ( ABS( u(ku+1,j,i) - u_0(j,i) )  > &
1094                     ABS( u(ku+1,j,i) - u(ku-1,j,i) ) )  u_0(j,i) = u(ku-1,j,i)
1095                IF ( ABS( v(kv+1,j,i) - v_0(j,i) )  > &
1096                     ABS( v(kv+1,j,i) - v(kv-1,j,i) ) )  v_0(j,i) = v(kv-1,j,i)
1097
1098             ENDDO
1099          ENDDO
1100
1101          CALL exchange_horiz_2d( u_0 )
1102          CALL exchange_horiz_2d( v_0 )
1103
1104       ENDIF
1105
1106    END SUBROUTINE production_e_init
1107
1108 END MODULE production_e_mod
Note: See TracBrowser for help on using the repository browser.