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

Last change on this file since 785 was 760, checked in by raasch, 13 years ago

last commit documented

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