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

Last change on this file since 52 was 52, checked in by raasch, 15 years ago

preliminary version of wall_fluxes.f90 added

File size: 5.7 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$
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
29    REAL    ::  a, b, c1, c2, h1, h2, delta_p
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    delta_p   = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
37    wall_flux = 0.0
38
39!
40!-- All subsequent variables are computed for the respective location where
41!-- the relevant variable is defined
42    DO  k = nzb_w, nzt_w
43!
44!--    (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
45       rifs    = rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2))
46       u_i     = a * u(k,j,i) +  &
47            c1 * 0.25 * ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
48       v_i     = b * v(k,j,i) +  &
49            c2 * 0.25 * ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
50       ws      = ( c1 + c2 ) * w(k,j,i) +  &
51            a * 0.25 * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) + &
52            b * 0.25 * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) )
53       pt_i    = 0.5 * ( pt(k,j,i) +  &
54            a *  pt(k,j,i-1) + b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
55       pts     = pt_i - hom(k,1,4,0)
56       wspts   = ws * pts
57!
58!--    (2) Compute wall-parallel absolute velocity vel_total
59       vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
60!
61!--    (3) Compute wall friction velocity us_wall
62       IF ( rifs >= 0.0 )  THEN
63!
64!--       Stable stratification (and neutral)
65          us_wall = kappa * vel_total / (                         &
66                    LOG( delta_p / z0(j,i) ) +                    &
67                    5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p  &
68                                        )
69       ELSE
70!
71!--       Unstable stratification
72          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
73          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / delta_p * z0(j,i) ) )
74!
75!--       If a borderline case occurs, the formula for stable stratification
76!--       must be used anyway, or else a zero division would occur in the
77!--       argument of the logarithm.
78          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
79             us_wall = kappa * vel_total / (                           &
80                       LOG( delta_p / z0(j,i) ) +                         &
81                       5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p &
82                                                   )
83          ELSE
84             us_wall = kappa * vel_total / (                           &
85                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
86                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
87                                                   )
88          ENDIF
89       ENDIF
90!
91!--    (4) Compute delta_p/L (corresponds to neutral Richardson flux number
92!--        rifs)
93       rifs = -1.0 * delta_p * kappa * g * wspts / &
94                       ( pt_i * ( us_wall**3 + 1E-30 ) )
95!
96!--    Limit the value range of the Richardson numbers.
97!--    This is necessary for very small velocities (u,w --> 0), because
98!--    the absolute value of rif can then become very large, which in
99!--    consequence would result in very large shear stresses and very
100!--    small momentum fluxes (both are generally unrealistic).
101       IF ( rifs < rif_min )  rifs = rif_min
102       IF ( rifs > rif_max )  rifs = rif_max
103!
104!--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
105       IF ( rifs >= 0.0 )  THEN
106!
107!--       Stable stratification (and neutral)
108          wall_flux(k) = kappa *  &
109               ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / (  &
110               LOG( delta_p / z0(j,i) ) +                                &
111               5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        &
112                                                                         )
113       ELSE
114!
115!--       Unstable stratification
116          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
117          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / delta_p * z0(j,i) ) )
118!
119!--       If a borderline case occurs, the formula for stable stratification
120!--       must be used anyway, or else a zero division would occur in the
121!--       argument of the logarithm.
122          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
123             wall_flux(k) = kappa *  &
124                  ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / (  &
125                  LOG( delta_p / z0(j,i) ) +                                &
126                  5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        &
127                                                                            )
128          ELSE
129             wall_flux(k) = kappa *  &
130                  ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / (  &
131                  LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) +           &
132                  2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                            &
133                                                                            )
134          ENDIF
135       ENDIF
136       wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
137
138!
139!--    store rifs for next time step
140       rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2)) = rifs
141
142    ENDDO
143 END SUBROUTINE wall_fluxes
Note: See TracBrowser for help on using the repository browser.