Changeset 56 for palm/trunk/SOURCE/wall_fluxes.f90
- Timestamp:
- Mar 8, 2007 1:57:07 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wall_fluxes.f90
r55 r56 1 SUBROUTINE wall_fluxes( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )1 MODULE wall_fluxes_mod 2 2 !------------------------------------------------------------------------------! 3 3 ! Actual revisions: … … 15 15 ! similarity. 16 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 17 ! The all-gridpoint version of wall_fluxes_e is not used so far, because 18 ! it gives slightly different results from the ij-version for some unknown 19 ! reason. 20 !------------------------------------------------------------------------------! 21 PRIVATE 22 PUBLIC wall_fluxes, wall_fluxes_e 23 24 INTERFACE wall_fluxes 25 MODULE PROCEDURE wall_fluxes 26 MODULE PROCEDURE wall_fluxes_ij 27 END INTERFACE wall_fluxes 28 29 INTERFACE wall_fluxes_e 30 MODULE PROCEDURE wall_fluxes_e 31 MODULE PROCEDURE wall_fluxes_e_ij 32 END INTERFACE wall_fluxes_e 33 34 CONTAINS 35 36 !------------------------------------------------------------------------------! 37 ! Call for all grid points 38 !------------------------------------------------------------------------------! 39 SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, ixp, jyp, nzb_uvw_inner, & 40 nzb_uvw_outer, wall ) 41 42 USE arrays_3d 43 USE control_parameters 44 USE grid_variables 45 USE indices 46 USE statistics 47 USE user 48 49 IMPLICIT NONE 50 51 INTEGER :: i, ixp, j, jyp, k, wall_index 52 53 INTEGER, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: nzb_uvw_inner, & 54 nzb_uvw_outer 55 REAL :: a, b, c1, c2, h1, h2, zp 56 REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts 57 58 REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: wall 59 REAL, DIMENSION(nzb:nzt+1,nys:nyn+jyp,nxl:nxr+ixp) :: wall_flux 60 61 62 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 63 wall_flux = 0.0 64 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 65 66 DO i = nxl, nxr+ixp 67 DO j = nys, nyn+jyp 68 69 IF ( wall(j,i) /= 0.0 ) THEN 70 ! 71 !-- All subsequent variables are computed for the respective 72 !-- location where the relevant variable is defined 73 DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i) 74 75 ! 76 !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' 77 rifs = rif_wall(k,j,i,wall_index) 78 79 u_i = a * u(k,j,i) + c1 * 0.25 * & 80 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 81 82 v_i = b * v(k,j,i) + c2 * 0.25 * & 83 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 84 85 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( & 86 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 87 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 88 ) 89 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + & 90 b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) ) 91 92 pts = pt_i - hom(k,1,4,0) 93 wspts = ws * pts 94 95 ! 96 !-- (2) Compute wall-parallel absolute velocity vel_total 97 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 98 99 ! 100 !-- (3) Compute wall friction velocity us_wall 101 IF ( rifs >= 0.0 ) THEN 102 103 ! 104 !-- Stable stratification (and neutral) 105 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 106 5.0 * rifs * ( zp - z0(j,i) ) / zp & 107 ) 108 ELSE 109 110 ! 111 !-- Unstable stratification 112 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 113 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 114 115 ! 116 !-- If a borderline case occurs, the formula for stable 117 !-- stratification must be used anyway, or else a zero 118 !-- division would occur in the argument of the logarithm. 119 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 120 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 121 5.0 * rifs * ( zp - z0(j,i) ) / zp & 122 ) 123 ELSE 124 us_wall = kappa * vel_total / ( & 125 LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & 126 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 127 ) 128 ENDIF 129 130 ENDIF 131 132 ! 133 !-- (4) Compute zp/L (corresponds to neutral Richardson flux 134 !-- number rifs) 135 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * & 136 ( us_wall**3 + 1E-30 ) ) 137 138 ! 139 !-- Limit the value range of the Richardson numbers. 140 !-- This is necessary for very small velocities (u,w --> 0), 141 !-- because the absolute value of rif can then become very 142 !-- large, which in consequence would result in very large 143 !-- shear stresses and very small momentum fluxes (both are 144 !-- generally unrealistic). 145 IF ( rifs < rif_min ) rifs = rif_min 146 IF ( rifs > rif_max ) rifs = rif_max 147 148 ! 149 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 150 IF ( rifs >= 0.0 ) THEN 151 152 ! 153 !-- Stable stratification (and neutral) 154 wall_flux(k,j,i) = kappa * & 155 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 156 ( LOG( zp / z0(j,i) ) + & 157 5.0 * rifs * ( zp - z0(j,i) ) / zp & 158 ) 159 ELSE 160 161 ! 162 !-- Unstable stratification 163 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 164 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 165 166 ! 167 !-- If a borderline case occurs, the formula for stable 168 !-- stratification must be used anyway, or else a zero 169 !-- division would occur in the argument of the logarithm. 170 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 171 wall_flux(k,j,i) = kappa * & 172 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 173 ( LOG( zp / z0(j,i) ) + & 174 5.0 * rifs * ( zp - z0(j,i) ) / zp & 175 ) 176 ELSE 177 wall_flux(k,j,i) = kappa * & 178 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 179 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )& 180 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 181 ) 182 ENDIF 183 184 ENDIF 185 wall_flux(k,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i)) 186 187 ! 188 !-- store rifs for next time step 189 rif_wall(k,j,i,wall_index) = rifs 190 191 ENDDO 192 193 ENDIF 194 195 ENDDO 196 ENDDO 197 198 END SUBROUTINE wall_fluxes 199 200 201 202 !------------------------------------------------------------------------------! 203 ! Call for all grid point i,j 204 !------------------------------------------------------------------------------! 205 SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) 206 207 USE arrays_3d 208 USE control_parameters 209 USE grid_variables 210 USE indices 211 USE statistics 212 USE user 213 214 IMPLICIT NONE 215 216 INTEGER :: i, j, k, nzb_w, nzt_w, wall_index 217 REAL :: a, b, c1, c2, h1, h2, zp 218 219 REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts 220 221 REAL, DIMENSION(nzb:nzt+1) :: wall_flux 222 223 224 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 225 wall_flux = 0.0 226 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 227 228 ! 229 !-- All subsequent variables are computed for the respective location where 230 !-- the relevant variable is defined 231 DO k = nzb_w, nzt_w 232 233 ! 234 !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' 235 rifs = rif_wall(k,j,i,wall_index) 236 237 u_i = a * u(k,j,i) + c1 * 0.25 * & 238 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 239 240 v_i = b * v(k,j,i) + c2 * 0.25 * & 241 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 242 243 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( & 244 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 245 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 246 ) 247 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) & 248 + ( c1 + c2 ) * pt(k+1,j,i) ) 249 250 pts = pt_i - hom(k,1,4,0) 251 wspts = ws * pts 252 253 ! 254 !-- (2) Compute wall-parallel absolute velocity vel_total 255 vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 ) 256 257 ! 258 !-- (3) Compute wall friction velocity us_wall 259 IF ( rifs >= 0.0 ) THEN 260 261 ! 262 !-- Stable stratification (and neutral) 90 263 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 91 264 5.0 * rifs * ( zp - z0(j,i) ) / zp & 92 265 ) 93 266 ELSE 94 us_wall = kappa * vel_total / ( & 267 268 ! 269 !-- Unstable stratification 270 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 271 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 272 273 ! 274 !-- If a borderline case occurs, the formula for stable stratification 275 !-- must be used anyway, or else a zero division would occur in the 276 !-- argument of the logarithm. 277 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 278 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 279 5.0 * rifs * ( zp - z0(j,i) ) / zp & 280 ) 281 ELSE 282 us_wall = kappa * vel_total / ( & 95 283 LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & 96 284 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 97 ) 285 ) 286 ENDIF 287 98 288 ENDIF 99 289 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 290 ! 291 !-- (4) Compute zp/L (corresponds to neutral Richardson flux number 292 !-- rifs) 293 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) ) 294 295 ! 296 !-- Limit the value range of the Richardson numbers. 297 !-- This is necessary for very small velocities (u,w --> 0), because 298 !-- the absolute value of rif can then become very large, which in 299 !-- consequence would result in very large shear stresses and very 300 !-- small momentum fluxes (both are generally unrealistic). 301 IF ( rifs < rif_min ) rifs = rif_min 302 IF ( rifs > rif_max ) rifs = rif_max 303 304 ! 305 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 306 IF ( rifs >= 0.0 ) THEN 307 308 ! 309 !-- Stable stratification (and neutral) 139 310 wall_flux(k) = kappa * & 140 311 ( 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&312 ( LOG( zp / z0(j,i) ) + & 313 5.0 * rifs * ( zp - z0(j,i) ) / zp & 143 314 ) 144 315 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 ) 316 317 ! 318 !-- Unstable stratification 319 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 320 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 321 322 ! 323 !-- If a borderline case occurs, the formula for stable stratification 324 !-- must be used anyway, or else a zero division would occur in the 325 !-- argument of the logarithm. 326 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 327 wall_flux(k) = kappa * & 328 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 329 ( LOG( zp / z0(j,i) ) + & 330 5.0 * rifs * ( zp - z0(j,i) ) / zp & 331 ) 332 ELSE 333 wall_flux(k) = kappa * & 334 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 335 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )& 336 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 337 ) 338 ENDIF 339 150 340 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 ! 341 wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) ) 342 343 ! 344 !-- store rifs for next time step 345 rif_wall(k,j,i,wall_index) = rifs 346 347 ENDDO 348 349 END SUBROUTINE wall_fluxes_ij 350 351 352 353 !------------------------------------------------------------------------------! 354 ! Call for all grid points 355 !------------------------------------------------------------------------------! 356 SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall ) 357 358 !------------------------------------------------------------------------------! 175 359 ! Description: 176 360 ! ------------ … … 180 364 !------------------------------------------------------------------------------! 181 365 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 / & 366 USE arrays_3d 367 USE control_parameters 368 USE grid_variables 369 USE indices 370 USE statistics 371 USE user 372 373 IMPLICIT NONE 374 375 INTEGER :: i, j, k, kk, wall_index 376 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 377 378 REAL :: rifs 379 380 REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: wall 381 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wall_flux 382 383 384 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 385 wall_flux = 0.0 386 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 387 388 DO i = nxl, nxr 389 DO j = nys, nyn 390 391 IF ( wall(j,i) /= 0.0 ) THEN 392 ! 393 !-- All subsequent variables are computed for the respective 394 !-- location where the relevant variable is defined 395 DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2 396 397 ! 398 !-- (1) Compute rifs 399 IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN 400 kk = nzb_diff_s_inner(j,i)-1 401 ELSE 402 kk = k-1 403 ENDIF 404 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 405 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 406 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 407 ) 408 409 ! 410 !-- Skip (2) to (4) of wall_fluxes, because here rifs is 411 !-- already available from (1) 412 413 ! 414 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 415 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & 416 b * ( v(k,j,i) + v(k,j+1,i) ) + & 417 (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & 418 ) 419 420 IF ( rifs >= 0.0 ) THEN 421 422 ! 423 !-- Stable stratification (and neutral) 424 wall_flux(k,j,i) = kappa * vel_zp / & 425 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp ) 426 ELSE 427 428 ! 429 !-- Unstable stratification 430 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 431 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 432 433 ! 434 !-- If a borderline case occurs, the formula for stable 435 !-- stratification must be used anyway, or else a zero 436 !-- division would occur in the argument of the logarithm. 437 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 438 wall_flux(k,j,i) = kappa * vel_zp / & 439 ( LOG( zp / z0(j,i) ) + & 440 5.0 * rifs * ( zp - z0(j,i) ) / zp & 441 ) 442 ELSE 443 wall_flux(k,j,i) = kappa * vel_zp / & 255 444 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & 256 445 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 257 446 ) 447 ENDIF 448 449 ENDIF 450 wall_flux(k,j,i) = wall_flux(k,j,i) * ABS( wall_flux(k,j,i) ) 451 452 ! 453 !-- Store rifs for next time step 454 rif_wall(k,j,i,wall_index) = rifs 455 456 ENDDO 457 458 ENDIF 459 460 ENDDO 461 ENDDO 462 463 END SUBROUTINE wall_fluxes_e 464 465 466 467 !------------------------------------------------------------------------------! 468 ! Call for grid point i,j 469 !------------------------------------------------------------------------------! 470 SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) 471 472 USE arrays_3d 473 USE control_parameters 474 USE grid_variables 475 USE indices 476 USE statistics 477 USE user 478 479 IMPLICIT NONE 480 481 INTEGER :: i, j, k, kk, nzb_w, nzt_w, wall_index 482 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 483 484 REAL :: rifs 485 486 REAL, DIMENSION(nzb:nzt+1) :: wall_flux 487 488 489 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 490 wall_flux = 0.0 491 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 492 493 ! 494 !-- All subsequent variables are computed for the respective location where 495 !-- the relevant variable is defined 496 DO k = nzb_w, nzt_w 497 498 ! 499 !-- (1) Compute rifs 500 IF ( k == nzb_w ) THEN 501 kk = nzb_w 502 ELSE 503 kk = k-1 258 504 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 505 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 506 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 507 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 508 ) 509 510 ! 511 !-- Skip (2) to (4) of wall_fluxes, because here rifs is already available 512 !-- from (1) 513 514 ! 515 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 516 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & 517 b * ( v(k,j,i) + v(k,j+1,i) ) + & 518 (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & 519 ) 520 521 IF ( rifs >= 0.0 ) THEN 522 523 ! 524 !-- Stable stratification (and neutral) 525 wall_flux(k) = kappa * vel_zp / & 526 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp ) 527 ELSE 528 529 ! 530 !-- Unstable stratification 531 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 532 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 533 534 ! 535 !-- If a borderline case occurs, the formula for stable stratification 536 !-- must be used anyway, or else a zero division would occur in the 537 !-- argument of the logarithm. 538 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 539 wall_flux(k) = kappa * vel_zp / & 540 ( LOG( zp / z0(j,i) ) + & 541 5.0 * rifs * ( zp - z0(j,i) ) / zp & 542 ) 543 ELSE 544 wall_flux(k) = kappa * vel_zp / & 545 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & 546 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 547 ) 548 ENDIF 549 550 ENDIF 551 wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) ) 552 553 ! 554 !-- Store rifs for next time step 555 rif_wall(k,j,i,wall_index) = rifs 556 557 ENDDO 558 559 END SUBROUTINE wall_fluxes_e_ij 560 561 END MODULE wall_fluxes_mod
Note: See TracChangeset
for help on using the changeset viewer.