Changeset 53 for palm/trunk/SOURCE/wall_fluxes.f90
- Timestamp:
- Mar 7, 2007 12:33:47 PM (17 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wall_fluxes.f90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wall_fluxes.f90
r52 r53 26 26 IMPLICIT NONE 27 27 28 INTEGER :: i, j, k, nzb_w, nzt_w 29 REAL :: a, b, c1, c2, h1, h2, delta_p28 INTEGER :: i, j, k, nzb_w, nzt_w, wall_index 29 REAL :: a, b, c1, c2, h1, h2, zp 30 30 31 31 REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts … … 34 34 35 35 36 delta_p = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 37 wall_flux = 0.0 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 ) 38 39 39 40 ! … … 41 42 !-- the relevant variable is defined 42 43 DO k = nzb_w, nzt_w 44 43 45 ! 44 46 !-- (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 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 57 65 ! 58 66 !-- (2) Compute wall-parallel absolute velocity vel_total 59 67 vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 ) 68 60 69 ! 61 70 !-- (3) Compute wall friction velocity us_wall 62 71 IF ( rifs >= 0.0 ) THEN 72 63 73 ! 64 74 !-- 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 & 75 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 76 5.0 * rifs * ( zp - z0(j,i) ) / zp & 68 77 ) 69 78 ELSE 79 70 80 ! 71 81 !-- Unstable stratification 72 82 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) ) ) 83 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 84 74 85 ! 75 86 !-- If a borderline case occurs, the formula for stable stratification … … 77 88 !-- argument of the logarithm. 78 89 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 ) 90 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 91 5.0 * rifs * ( zp - z0(j,i) ) / zp & 92 ) 83 93 ELSE 84 us_wall = kappa * vel_total / ( &94 us_wall = kappa * vel_total / ( & 85 95 LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & 86 96 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 87 97 ) 88 98 ENDIF 89 ENDIF 90 ! 91 !-- (4) Compute delta_p/L (corresponds to neutral Richardson flux number 99 100 ENDIF 101 102 ! 103 !-- (4) Compute zp/L (corresponds to neutral Richardson flux number 92 104 !-- rifs) 93 rifs = -1.0 * delta_p * kappa * g * wspts / &94 ( pt_i * ( us_wall**3 + 1E-30 ) ) 105 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) ) 106 95 107 ! 96 108 !-- Limit the value range of the Richardson numbers. … … 101 113 IF ( rifs < rif_min ) rifs = rif_min 102 114 IF ( rifs > rif_max ) rifs = rif_max 115 103 116 ! 104 117 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 105 118 IF ( rifs >= 0.0 ) THEN 119 106 120 ! 107 121 !-- 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 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 114 129 ! 115 130 !-- Unstable stratification 116 131 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) ) ) 132 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 133 118 134 ! 119 135 !-- If a borderline case occurs, the formula for stable stratification … … 121 137 !-- argument of the logarithm. 122 138 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 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 ) 128 144 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 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 ) 134 150 ENDIF 151 135 152 ENDIF 136 153 wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) ) … … 138 155 ! 139 156 !-- store rifs for next time step 140 rif_wall(k,j,i, NINT(a+2*b+3*c1+4*c2)) = rifs157 rif_wall(k,j,i,wall_index) = rifs 141 158 142 159 ENDDO 160 143 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, 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 IF ( rifs >= 0.0 ) THEN 227 228 ! 229 !-- Stable stratification (and neutral) 230 wall_flux(k) = kappa * & 231 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 232 ( LOG( zp / z0(j,i) ) + & 233 5.0 * rifs * ( zp - z0(j,i) ) / zp & 234 ) 235 ELSE 236 237 ! 238 !-- Unstable stratification 239 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 240 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 241 242 ! 243 !-- If a borderline case occurs, the formula for stable stratification 244 !-- must be used anyway, or else a zero division would occur in the 245 !-- argument of the logarithm. 246 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 247 wall_flux(k) = kappa * & 248 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 249 ( LOG( zp / z0(j,i) ) + & 250 5.0 * rifs * ( zp - z0(j,i) ) / zp & 251 ) 252 ELSE 253 wall_flux(k) = kappa * & 254 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 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 TracChangeset
for help on using the changeset viewer.
