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

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