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

Last change on this file since 550 was 484, checked in by raasch, 14 years ago

typo in file headers removed

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