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

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