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

Last change on this file since 325 was 226, checked in by raasch, 16 years ago

preparations for the next release

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