source: palm/trunk/SOURCE/wall_fluxes.f90 @ 55

Last change on this file since 55 was 55, checked in by raasch, 17 years ago

small bugs removed

  • Property svn:keywords set to Id
File size: 8.8 KB
Line 
1 SUBROUTINE wall_fluxes( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
2!------------------------------------------------------------------------------!
3! Actual revisions:
4! -----------------
5!
6!
7! Former revisions:
8! -----------------
9! $Id: wall_fluxes.f90 55 2007-03-08 04:12:41Z raasch $
10! Initial version (2007/03/07)
11!
12! Description:
13! ------------
14! Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
15! similarity.
16! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
17!------------------------------------------------------------------------------!
18
19    USE arrays_3d
20    USE control_parameters
21    USE grid_variables
22    USE indices
23    USE statistics
24    USE user
25
26    IMPLICIT NONE
27
28    INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
29    REAL    ::  a, b, c1, c2, h1, h2, zp
30
31    REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
32
33    REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
34
35
36    zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
37    wall_flux  = 0.0
38    wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
39
40!
41!-- All subsequent variables are computed for the respective location where
42!-- the relevant variable is defined
43    DO  k = nzb_w, nzt_w
44
45!
46!--    (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
47       rifs  = rif_wall(k,j,i,wall_index)
48
49       u_i   = a * u(k,j,i) + c1 * 0.25 * &
50               ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
51
52       v_i   = b * v(k,j,i) + c2 * 0.25 * &
53               ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
54
55       ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                             &
56                   a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
57                 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
58                                               )
59       pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
60                       + ( c1 + c2 ) * pt(k+1,j,i) )
61
62       pts   = pt_i - hom(k,1,4,0)
63       wspts = ws * pts
64
65!
66!--    (2) Compute wall-parallel absolute velocity vel_total
67       vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
68
69!
70!--    (3) Compute wall friction velocity us_wall
71       IF ( rifs >= 0.0 )  THEN
72
73!
74!--       Stable stratification (and neutral)
75          us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +               &
76                                          5.0 * rifs * ( zp - z0(j,i) ) / zp  &
77                                        )
78       ELSE
79
80!
81!--       Unstable stratification
82          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
83          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
84
85!
86!--       If a borderline case occurs, the formula for stable stratification
87!--       must be used anyway, or else a zero division would occur in the
88!--       argument of the logarithm.
89          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
90             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
91                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
92                                           )
93          ELSE
94             us_wall = kappa * vel_total / (                                   &
95                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
96                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
97                                           )
98          ENDIF
99
100       ENDIF
101
102!
103!--    (4) Compute zp/L (corresponds to neutral Richardson flux number
104!--        rifs)
105       rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) )
106
107!
108!--    Limit the value range of the Richardson numbers.
109!--    This is necessary for very small velocities (u,w --> 0), because
110!--    the absolute value of rif can then become very large, which in
111!--    consequence would result in very large shear stresses and very
112!--    small momentum fluxes (both are generally unrealistic).
113       IF ( rifs < rif_min )  rifs = rif_min
114       IF ( rifs > rif_max )  rifs = rif_max
115
116!
117!--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
118       IF ( rifs >= 0.0 )  THEN
119
120!
121!--       Stable stratification (and neutral)
122          wall_flux(k) = kappa *                                          &
123                         ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
124                         (  LOG( zp / z0(j,i) ) +                         &
125                            5.0 * rifs * ( zp - z0(j,i) ) / zp            &
126                         )
127       ELSE
128
129!
130!--       Unstable stratification
131          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
132          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
133
134!
135!--       If a borderline case occurs, the formula for stable stratification
136!--       must be used anyway, or else a zero division would occur in the
137!--       argument of the logarithm.
138          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
139             wall_flux(k) = kappa *                                          &
140                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
141                            ( LOG( zp / z0(j,i) ) +                          &
142                              5.0 * rifs * ( zp - z0(j,i) ) / zp             &
143                            )
144          ELSE
145             wall_flux(k) = kappa *                                            &
146                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /   &
147                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
148                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
149                            )
150          ENDIF
151
152       ENDIF
153       wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
154
155!
156!--    store rifs for next time step
157       rif_wall(k,j,i,wall_index) = rifs
158
159    ENDDO
160
161 END SUBROUTINE wall_fluxes
162
163
164
165 SUBROUTINE wall_fluxes_e( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
166!------------------------------------------------------------------------------!
167! Actual revisions:
168! -----------------
169!
170!
171! Former revisions:
172! -----------------
173! Initial version (2007/03/07)
174!
175! Description:
176! ------------
177! Calculates momentum fluxes at vertical walls for routine production_e
178! assuming Monin-Obukhov similarity.
179! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
180!------------------------------------------------------------------------------!
181
182    USE arrays_3d
183    USE control_parameters
184    USE grid_variables
185    USE indices
186    USE statistics
187    USE user
188
189    IMPLICIT NONE
190
191    INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
192    REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
193
194    REAL ::  rifs
195
196    REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
197
198
199    zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
200    wall_flux  = 0.0
201    wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
202
203!
204!-- All subsequent variables are computed for the respective location where
205!-- the relevant variable is defined
206    DO  k = nzb_w, nzt_w
207
208!
209!--    (1) Compute rifs
210       IF ( k == nzb_w )  THEN
211          kk = nzb_w
212       ELSE
213          kk = k-1
214       ENDIF
215       rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
216                       a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
217                       c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
218                     )
219
220!
221!--    Skip (2) to (4) of wall_fluxes, because here rifs is already available
222!--    from (1)
223
224!
225!--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
226       vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
227                              b * ( v(k,j,i) + v(k,j+1,i) ) +  &
228                        (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
229                      )
230
231       IF ( rifs >= 0.0 )  THEN
232
233!
234!--       Stable stratification (and neutral)
235          wall_flux(k) = kappa *  vel_zp / &
236                         ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
237       ELSE
238
239!
240!--       Unstable stratification
241          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
242          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
243
244!
245!--       If a borderline case occurs, the formula for stable stratification
246!--       must be used anyway, or else a zero division would occur in the
247!--       argument of the logarithm.
248          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
249             wall_flux(k) = kappa * vel_zp /                                 &
250                            ( LOG( zp / z0(j,i) ) +                          &
251                              5.0 * rifs * ( zp - z0(j,i) ) / zp             &
252                            )
253          ELSE
254             wall_flux(k) = kappa * vel_zp /                                   &
255                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
256                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
257                            )
258          ENDIF
259
260       ENDIF
261       wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) )
262
263!
264!--    store rifs for next time step
265       rif_wall(k,j,i,wall_index) = rifs
266
267    ENDDO
268
269 END SUBROUTINE wall_fluxes_e
Note: See TracBrowser for help on using the repository browser.