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

Last change on this file since 388 was 388, checked in by raasch, 12 years ago

in-situ AND potential density are calculated and used in the ocean version

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