Changeset 187 for palm/trunk/SOURCE/wall_fluxes.f90
 Timestamp:
 Aug 6, 2008 4:25:09 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/wall_fluxes.f90
r75 r187 3 3 ! Actual revisions: 4 4 !  5 ! 5 ! Bugfix: change definition of us_wall from 1D to 2D: 6 ! Modification of the evaluation of the vertical turbulent momentum 7 ! fluxes u'w' and v'w'; the first usws that is computed corresponds 8 ! to u'w'/u* and not as priorily assumed to (u'w')**0.5, the first 9 ! vsws that is computed corresponds to v'w'/u* and not as priorily 10 ! assumed to (v'w')**0.5. Therefore, the intermediate result for 11 ! usws has to be multiplied by u* instead by itself in order to 12 ! get u'w'. Accordingly, the intermediate result for vsws has to be 13 ! multiplied by u* instead by itself in order to get v'w'. 14 ! This requires the calculation of us_wall (and vel_total, u_i, v_i, ws) 15 ! also in wall_fluxes_e. 16 ! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed 17 ! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for 18 ! consistency with subroutine wall_fluxes 19 ! Change: Modification of the integrated version of the profile function for 20 ! momentum for unstable stratification 6 21 ! 7 22 ! Former revisions: … … 69 84 ! 70 85 ! All subsequent variables are computed for the respective 71 ! location where the re levant variable is defined86 ! location where the respective flux is defined. 72 87 DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i) 73 88 … … 109 124 ! 110 125 ! Unstable stratification 111 h1 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs ) ) 112 h2 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs / zp * z0(j,i) )) 113 114 ! 115 ! If a borderline case occurs, the formula for stable 116 ! stratification must be used anyway, or else a zero 117 ! division would occur in the argument of the logarithm. 118 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 119 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 120 5.0 * rifs * ( zp  z0(j,i) ) / zp & 121 ) 122 ELSE 123 us_wall = kappa * vel_total / ( & 124 LOG( (1.0+h2) / (1.0h2) * (1.0h1) / (1.0+h1) ) + & 125 2.0 * ( ATAN( h2 )  ATAN( h1 ) ) & 126 ) 127 ENDIF 128 126 h1 = SQRT( SQRT( 1.0  16.0 * rifs ) ) 127 h2 = SQRT( SQRT( 1.0  16.0 * rifs * z0(j,i) / zp ) ) 128 129 us_wall = kappa * vel_total / ( & 130 LOG( zp / z0(j,i) )  & 131 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 132 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 133 2.0 * ( ATAN( h1 )  ATAN( h2 ) ) & 134 ) 129 135 ENDIF 130 136 … … 160 166 ! 161 167 ! Unstable stratification 162 h1 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs ) ) 163 h2 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs / zp * z0(j,i) )) 164 165 ! 166 ! If a borderline case occurs, the formula for stable 167 ! stratification must be used anyway, or else a zero 168 ! division would occur in the argument of the logarithm. 169 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 170 wall_flux(k,j,i) = kappa * & 171 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 172 ( LOG( zp / z0(j,i) ) + & 173 5.0 * rifs * ( zp  z0(j,i) ) / zp & 174 ) 175 ELSE 176 wall_flux(k,j,i) = kappa * & 177 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 178 ( LOG( (1.0+h2) / (1.0h2) * (1.0h1) / (1.0+h1) )& 179 + 2.0 * ( ATAN( h2 )  ATAN( h1 ) ) & 180 ) 181 ENDIF 182 168 h1 = SQRT( SQRT( 1.0  16.0 * rifs ) ) 169 h2 = SQRT( SQRT( 1.0  16.0 * rifs * z0(j,i) / zp ) ) 170 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 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 175 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 176 2.0 * ( ATAN( h1 )  ATAN( h2 ) ) & 177 ) 183 178 ENDIF 184 wall_flux(k,j,i) = wall_flux(k,j,i) * ABS(wall_flux(k,j,i))179 wall_flux(k,j,i) = wall_flux(k,j,i) * us_wall 185 180 186 181 ! … … 226 221 ! 227 222 ! All subsequent variables are computed for the respective location where 228 ! the re levant variable is defined223 ! the respective flux is defined. 229 224 DO k = nzb_w, nzt_w 230 225 … … 266 261 ! 267 262 ! Unstable stratification 268 h1 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs ) ) 269 h2 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs / zp * z0(j,i) ) ) 270 271 ! 272 ! If a borderline case occurs, the formula for stable stratification 273 ! must be used anyway, or else a zero division would occur in the 274 ! argument of the logarithm. 275 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 276 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 277 5.0 * rifs * ( zp  z0(j,i) ) / zp & 278 ) 279 ELSE 280 us_wall = kappa * vel_total / ( & 281 LOG( (1.0+h2) / (1.0h2) * (1.0h1) / (1.0+h1) ) + & 282 2.0 * ( ATAN( h2 )  ATAN( h1 ) ) & 283 ) 284 ENDIF 285 263 h1 = SQRT( SQRT( 1.0  16.0 * rifs ) ) 264 h2 = SQRT( SQRT( 1.0  16.0 * rifs * z0(j,i) / zp ) ) 265 266 us_wall = kappa * vel_total / ( & 267 LOG( zp / z0(j,i) )  & 268 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 269 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 270 2.0 * ( ATAN( h1 )  ATAN( h2 ) ) & 271 ) 286 272 ENDIF 287 273 … … 315 301 ! 316 302 ! Unstable stratification 317 h1 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs ) ) 318 h2 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs / zp * z0(j,i) ) ) 319 320 ! 321 ! If a borderline case occurs, the formula for stable stratification 322 ! must be used anyway, or else a zero division would occur in the 323 ! argument of the logarithm. 324 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 325 wall_flux(k) = kappa * & 326 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 327 ( LOG( zp / z0(j,i) ) + & 328 5.0 * rifs * ( zp  z0(j,i) ) / zp & 329 ) 330 ELSE 331 wall_flux(k) = kappa * & 332 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 333 ( LOG( (1.0+h2) / (1.0h2) * (1.0h1) / (1.0+h1) )& 334 + 2.0 * ( ATAN( h2 )  ATAN( h1 ) ) & 335 ) 336 ENDIF 337 303 h1 = SQRT( SQRT( 1.0  16.0 * rifs ) ) 304 h2 = SQRT( SQRT( 1.0  16.0 * rifs * z0(j,i) / zp ) ) 305 306 wall_flux(k) = kappa * & 307 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / ( & 308 LOG( zp / z0(j,i) )  & 309 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 310 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 311 2.0 * ( ATAN( h1 )  ATAN( h2 ) ) & 312 ) 338 313 ENDIF 339 wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) )314 wall_flux(k) = wall_flux(k) * us_wall 340 315 341 316 ! … … 371 346 372 347 INTEGER :: i, j, k, kk, wall_index 373 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 348 REAL :: a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, & 349 ws, zp 374 350 375 351 REAL :: rifs … … 388 364 IF ( wall(j,i) /= 0.0 ) THEN 389 365 ! 390 ! All subsequent variables are computed for the respective 391 ! location where the relevant variable is defined 366 ! All subsequent variables are computed for scalar locations. 392 367 DO k = nzb_diff_s_inner(j,i)1, nzb_diff_s_outer(j,i)2 393 394 ! 395 ! (1) Compute rifs 368 ! 369 ! (1) Compute rifs, u_i, v_i, and ws 396 370 IF ( k == nzb_diff_s_inner(j,i)1 ) THEN 397 371 kk = nzb_diff_s_inner(j,i)1 … … 404 378 ) 405 379 406 ! 407 ! Skip (2) to (4) of wall_fluxes, because here rifs is 408 ! already available from (1) 409 380 u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) 381 v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) 382 ws = 0.5 * ( w(k,j,i) + w(k1,j,i) ) 383 ! 384 ! (2) Compute wallparallel absolute velocity vel_total and 385 ! interpolate appropriate velocity component vel_zp. 386 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 387 vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws ) 388 ! 389 ! (3) Compute wall friction velocity us_wall 390 IF ( rifs >= 0.0 ) THEN 391 392 ! 393 ! Stable stratification (and neutral) 394 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 395 5.0 * rifs * ( zp  z0(j,i) ) / zp & 396 ) 397 ELSE 398 399 ! 400 ! Unstable stratification 401 h1 = SQRT( SQRT( 1.0  16.0 * rifs ) ) 402 h2 = SQRT( SQRT( 1.0  16.0 * rifs * z0(j,i) / zp ) ) 403 404 us_wall = kappa * vel_total / ( & 405 LOG( zp / z0(j,i) )  & 406 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 407 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 408 2.0 * ( ATAN( h1 )  ATAN( h2 ) ) & 409 ) 410 ENDIF 411 412 ! 413 ! Skip step (4) of wall_fluxes, because here rifs is already 414 ! available from (1) 410 415 ! 411 416 ! (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 412 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + &413 b * ( v(k,j,i) + v(k,j+1,i) ) + &414 (c1+c2) * ( w(k,j,i) + w(k1,j,i) ) &415 )416 417 417 418 IF ( rifs >= 0.0 ) THEN … … 425 426 ! 426 427 ! Unstable stratification 427 h1 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs ) ) 428 h2 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs / zp * z0(j,i) )) 429 430 ! 431 ! If a borderline case occurs, the formula for stable 432 ! stratification must be used anyway, or else a zero 433 ! division would occur in the argument of the logarithm. 434 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 435 wall_flux(k,j,i) = kappa * vel_zp / & 436 ( LOG( zp / z0(j,i) ) + & 437 5.0 * rifs * ( zp  z0(j,i) ) / zp & 438 ) 439 ELSE 440 wall_flux(k,j,i) = kappa * vel_zp / & 441 ( LOG( (1.0+h2) / (1.0h2) * (1.0h1) / (1.0+h1) ) & 442 + 2.0 * ( ATAN( h2 )  ATAN( h1 ) ) & 443 ) 444 ENDIF 445 428 h1 = SQRT( SQRT( 1.0  16.0 * rifs ) ) 429 h2 = SQRT( SQRT( 1.0  16.0 * rifs * z0(j,i) / zp ) ) 430 431 wall_flux(k,j,i) = kappa * vel_zp / ( & 432 LOG( zp / z0(j,i) )  & 433 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 434 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 435 2.0 * ( ATAN( h1 )  ATAN( h2 ) ) & 436 ) 446 437 ENDIF 447 wall_flux(k,j,i) = wall_flux(k,j,i) * ABS( wall_flux(k,j,i) ) 448 449 ! 450 ! Store rifs for next time step 451 rif_wall(k,j,i,wall_index) = rifs 438 wall_flux(k,j,i) =  wall_flux(k,j,i) * us_wall 452 439 453 440 ENDDO … … 476 463 477 464 INTEGER :: i, j, k, kk, nzb_w, nzt_w, wall_index 478 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 465 REAL :: a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, & 466 ws, zp 479 467 480 468 REAL :: rifs … … 488 476 489 477 ! 490 ! All subsequent variables are computed for the respective location where 491 ! the relevant variable is defined 478 ! All subsequent variables are computed for scalar locations. 492 479 DO k = nzb_w, nzt_w 493 480 494 481 ! 495 ! (1) Compute rifs 482 ! (1) Compute rifs, u_i, v_i, and ws 496 483 IF ( k == nzb_w ) THEN 497 484 kk = nzb_w … … 504 491 ) 505 492 506 ! 507 ! Skip (2) to (4) of wall_fluxes, because here rifs is already available 508 ! from (1) 509 493 u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) 494 v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) 495 ws = 0.5 * ( w(k,j,i) + w(k1,j,i) ) 496 ! 497 ! (2) Compute wallparallel absolute velocity vel_total and 498 ! interpolate appropriate velocity component vel_zp. 499 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 500 vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws ) 501 ! 502 ! (3) Compute wall friction velocity us_wall 503 IF ( rifs >= 0.0 ) THEN 504 505 ! 506 ! Stable stratification (and neutral) 507 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 508 5.0 * rifs * ( zp  z0(j,i) ) / zp & 509 ) 510 ELSE 511 512 ! 513 ! Unstable stratification 514 h1 = SQRT( SQRT( 1.0  16.0 * rifs ) ) 515 h2 = SQRT( SQRT( 1.0  16.0 * rifs * z0(j,i) / zp ) ) 516 517 us_wall = kappa * vel_total / ( & 518 LOG( zp / z0(j,i) )  & 519 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 520 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 521 2.0 * ( ATAN( h1 )  ATAN( h2 ) ) & 522 ) 523 ENDIF 524 525 ! 526 ! Skip step (4) of wall_fluxes, because here rifs is already 527 ! available from (1) 510 528 ! 511 529 ! (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 530 ! First interpolate the velocity (this is different from 531 ! subroutine wall_fluxes because fluxes in subroutine 532 ! wall_fluxes_e are defined at scalar locations). 512 533 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & 513 534 b * ( v(k,j,i) + v(k,j+1,i) ) + & … … 525 546 ! 526 547 ! Unstable stratification 527 h1 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs ) ) 528 h2 = 1.0 / SQRT( SQRT( 1.0  16.0 * rifs / zp * z0(j,i) ) ) 529 530 ! 531 ! If a borderline case occurs, the formula for stable stratification 532 ! must be used anyway, or else a zero division would occur in the 533 ! argument of the logarithm. 534 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 535 wall_flux(k) = kappa * vel_zp / & 536 ( LOG( zp / z0(j,i) ) + & 537 5.0 * rifs * ( zp  z0(j,i) ) / zp & 538 ) 539 ELSE 540 wall_flux(k) = kappa * vel_zp / & 541 ( LOG( (1.0+h2) / (1.0h2) * (1.0h1) / (1.0+h1) ) & 542 + 2.0 * ( ATAN( h2 )  ATAN( h1 ) ) & 543 ) 544 ENDIF 545 548 h1 = SQRT( SQRT( 1.0  16.0 * rifs ) ) 549 h2 = SQRT( SQRT( 1.0  16.0 * rifs * z0(j,i) / zp ) ) 550 551 wall_flux(k) = kappa * vel_zp / ( & 552 LOG( zp / z0(j,i) )  & 553 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 554 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 555 2.0 * ( ATAN( h1 )  ATAN( h2 ) ) & 556 ) 546 557 ENDIF 547 wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) ) 548 549 ! 550 ! Store rifs for next time step 551 rif_wall(k,j,i,wall_index) = rifs 558 wall_flux(k) =  wall_flux(k) * us_wall 552 559 553 560 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.