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

Last change on this file since 722 was 668, checked in by suehring, 13 years ago

last commit documented

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