Changeset 1353 for palm/trunk/SOURCE/wall_fluxes.f90
- Timestamp:
- Apr 8, 2014 3:21:23 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wall_fluxes.f90 (modified) (41 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wall_fluxes.f90
r1321 r1353 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants provided with KIND-attribute 23 23 ! 24 24 ! Former revisions: … … 148 148 149 149 150 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )151 wall_flux = 0.0 150 zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx ) 151 wall_flux = 0.0_wp 152 152 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 153 153 … … 155 155 DO j = nys, nyn 156 156 157 IF ( wall(j,i) /= 0.0 ) THEN157 IF ( wall(j,i) /= 0.0_wp ) THEN 158 158 ! 159 159 !-- All subsequent variables are computed for the respective … … 165 165 rifs = rif_wall(k,j,i,wall_index) 166 166 167 u_i = a * u(k,j,i) + c1 * 0.25 *&167 u_i = a * u(k,j,i) + c1 * 0.25_wp * & 168 168 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 169 169 170 v_i = b * v(k,j,i) + c2 * 0.25 *&170 v_i = b * v(k,j,i) + c2 * 0.25_wp * & 171 171 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 172 172 173 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * (&173 ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( & 174 174 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 175 175 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 176 )177 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) +&176 ) 177 pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + & 178 178 b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) ) 179 179 … … 187 187 ! 188 188 !-- (3) Compute wall friction velocity us_wall 189 IF ( rifs >= 0.0 ) THEN189 IF ( rifs >= 0.0_wp ) THEN 190 190 191 191 ! 192 192 !-- Stable stratification (and neutral) 193 193 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 194 5.0* rifs * ( zp - z0(j,i) ) / zp &194 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 195 195 ) 196 196 ELSE … … 198 198 ! 199 199 !-- Unstable stratification 200 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )201 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )200 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 201 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 202 202 203 203 us_wall = kappa * vel_total / ( & 204 204 LOG( zp / z0(j,i) ) - & 205 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&206 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&207 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&205 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 206 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +& 207 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 208 208 ) 209 209 ENDIF … … 212 212 !-- (4) Compute zp/L (corresponds to neutral Richardson flux 213 213 !-- number rifs) 214 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i *&215 214 rifs = -1.0_wp * zp * kappa * g * wspts / & 215 ( pt_i * ( us_wall**3 + 1E-30 ) ) 216 216 217 217 ! … … 227 227 ! 228 228 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 229 IF ( rifs >= 0.0 ) THEN229 IF ( rifs >= 0.0_wp ) THEN 230 230 231 231 ! … … 234 234 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 235 235 ( LOG( zp / z0(j,i) ) + & 236 5.0 * rifs * ( zp - z0(j,i) ) / zp&236 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 237 237 ) 238 238 ELSE … … 240 240 ! 241 241 !-- Unstable stratification 242 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )243 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )242 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 243 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 244 244 245 245 wall_flux(k,j,i) = kappa * & 246 246 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / ( & 247 247 LOG( zp / z0(j,i) ) - & 248 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&249 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&250 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&248 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 249 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +& 250 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 251 251 ) 252 252 ENDIF … … 333 333 334 334 335 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )336 wall_flux = 0.0 335 zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx ) 336 wall_flux = 0.0_wp 337 337 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 338 338 … … 346 346 DO j = j_south, j_north 347 347 348 IF ( wall(j,i) /= 0.0 ) THEN348 IF ( wall(j,i) /= 0.0_wp ) THEN 349 349 ! 350 350 !-- All subsequent variables are computed for the respective … … 357 357 rifs = rif_wall(k,j,i,wall_index) 358 358 359 u_i = a * u(k,j,i) + c1 * 0.25 *&359 u_i = a * u(k,j,i) + c1 * 0.25_wp * & 360 360 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 361 361 362 v_i = b * v(k,j,i) + c2 * 0.25 *&362 v_i = b * v(k,j,i) + c2 * 0.25_wp * & 363 363 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 364 364 365 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * (&365 ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( & 366 366 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 367 367 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 368 )369 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) +&368 ) 369 pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + & 370 370 b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) ) 371 371 … … 379 379 ! 380 380 !-- (3) Compute wall friction velocity us_wall 381 IF ( rifs >= 0.0 ) THEN381 IF ( rifs >= 0.0_wp ) THEN 382 382 383 383 ! 384 384 !-- Stable stratification (and neutral) 385 385 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 386 5.0* rifs * ( zp - z0(j,i) ) / zp &386 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 387 387 ) 388 388 ELSE … … 390 390 ! 391 391 !-- Unstable stratification 392 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )393 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )392 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 393 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 394 394 395 395 us_wall = kappa * vel_total / ( & 396 396 LOG( zp / z0(j,i) ) - & 397 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&398 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&399 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&397 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 398 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +& 399 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 400 400 ) 401 401 ENDIF … … 404 404 !-- (4) Compute zp/L (corresponds to neutral Richardson flux 405 405 !-- number rifs) 406 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i *&407 406 rifs = -1.0_wp * zp * kappa * g * wspts / & 407 ( pt_i * ( us_wall**3 + 1E-30 ) ) 408 408 409 409 ! … … 419 419 ! 420 420 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 421 IF ( rifs >= 0.0 ) THEN421 IF ( rifs >= 0.0_wp ) THEN 422 422 423 423 ! … … 426 426 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 427 427 ( LOG( zp / z0(j,i) ) + & 428 5.0 * rifs * ( zp - z0(j,i) ) / zp&428 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 429 429 ) 430 430 ELSE … … 432 432 ! 433 433 !-- Unstable stratification 434 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )435 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )434 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 435 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 436 436 437 437 wall_flux(k,j,i) = kappa * & 438 438 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / ( & 439 439 LOG( zp / z0(j,i) ) - & 440 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&441 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&442 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&440 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 441 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +& 442 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 443 443 ) 444 444 ENDIF … … 511 511 512 512 513 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )514 wall_flux = 0.0 513 zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx ) 514 wall_flux = 0.0_wp 515 515 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 516 516 … … 524 524 rifs = rif_wall(k,j,i,wall_index) 525 525 526 u_i = a * u(k,j,i) + c1 * 0.25 *&526 u_i = a * u(k,j,i) + c1 * 0.25_wp * & 527 527 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 528 528 529 v_i = b * v(k,j,i) + c2 * 0.25 *&529 v_i = b * v(k,j,i) + c2 * 0.25_wp * & 530 530 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 531 531 532 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * (&532 ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( & 533 533 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 534 534 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 535 )536 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i)&535 ) 536 pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) & 537 537 + ( c1 + c2 ) * pt(k+1,j,i) ) 538 538 … … 546 546 ! 547 547 !-- (3) Compute wall friction velocity us_wall 548 IF ( rifs >= 0.0 ) THEN548 IF ( rifs >= 0.0_wp ) THEN 549 549 550 550 ! 551 551 !-- Stable stratification (and neutral) 552 552 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 553 5.0* rifs * ( zp - z0(j,i) ) / zp &553 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 554 554 ) 555 555 ELSE … … 557 557 ! 558 558 !-- Unstable stratification 559 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )560 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )559 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 560 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 561 561 562 562 us_wall = kappa * vel_total / ( & 563 563 LOG( zp / z0(j,i) ) - & 564 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&565 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&566 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&564 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 565 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) + & 566 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 567 567 ) 568 568 ENDIF … … 571 571 !-- (4) Compute zp/L (corresponds to neutral Richardson flux number 572 572 !-- rifs) 573 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) ) 573 rifs = -1.0_wp * zp * kappa * g * wspts / & 574 ( pt_i * (us_wall**3 + 1E-30) ) 574 575 575 576 ! … … 584 585 ! 585 586 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 586 IF ( rifs >= 0.0 ) THEN587 IF ( rifs >= 0.0_wp ) THEN 587 588 588 589 ! … … 591 592 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 592 593 ( LOG( zp / z0(j,i) ) + & 593 5.0 * rifs * ( zp - z0(j,i) ) / zp&594 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 594 595 ) 595 596 ELSE … … 597 598 ! 598 599 !-- Unstable stratification 599 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )600 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )600 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 601 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 601 602 602 603 wall_flux(k) = kappa * & 603 604 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / ( & 604 605 LOG( zp / z0(j,i) ) - & 605 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&606 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&607 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&606 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 607 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) + & 608 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 608 609 ) 609 610 ENDIF … … 680 681 681 682 682 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )683 wall_flux = 0.0 683 zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx ) 684 wall_flux = 0.0_wp 684 685 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 685 686 … … 687 688 DO j = nys, nyn 688 689 689 IF ( wall(j,i) /= 0.0 ) THEN690 IF ( wall(j,i) /= 0.0_wp ) THEN 690 691 ! 691 692 !-- All subsequent variables are computed for scalar locations. … … 698 699 kk = k-1 699 700 ENDIF 700 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 701 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 702 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 703 ) 704 705 u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) 706 v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) 707 ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) ) 701 rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + & 702 a * rif_wall(k,j,i+1,1) + & 703 b * rif_wall(k,j+1,i,2) + & 704 c1 * rif_wall(kk,j,i,3) + & 705 c2 * rif_wall(kk,j,i,4) & 706 ) 707 708 u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 709 v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 710 ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 708 711 ! 709 712 !-- (2) Compute wall-parallel absolute velocity vel_total and 710 713 !-- interpolate appropriate velocity component vel_zp. 711 714 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 712 vel_zp = 0.5* ( a * u_i + b * v_i + (c1+c2) * ws )715 vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws ) 713 716 ! 714 717 !-- (3) Compute wall friction velocity us_wall 715 IF ( rifs >= 0.0 ) THEN718 IF ( rifs >= 0.0_wp ) THEN 716 719 717 720 ! 718 721 !-- Stable stratification (and neutral) 719 722 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 720 5.0* rifs * ( zp - z0(j,i) ) / zp &723 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 721 724 ) 722 725 ELSE … … 724 727 ! 725 728 !-- Unstable stratification 726 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )727 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )729 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 730 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 728 731 729 732 us_wall = kappa * vel_total / ( & 730 733 LOG( zp / z0(j,i) ) - & 731 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&732 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&733 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&734 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 735 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +& 736 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 734 737 ) 735 738 ENDIF … … 741 744 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 742 745 743 IF ( rifs >= 0.0 ) THEN746 IF ( rifs >= 0.0_wp ) THEN 744 747 745 748 ! 746 749 !-- Stable stratification (and neutral) 747 wall_flux(k,j,i) = kappa * vel_zp /&748 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )750 wall_flux(k,j,i) = kappa * vel_zp / ( LOG( zp/z0(j,i) ) +& 751 5.0_wp * rifs * ( zp-z0(j,i) ) / zp ) 749 752 ELSE 750 753 751 754 ! 752 755 !-- Unstable stratification 753 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )754 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )756 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 757 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 755 758 756 759 wall_flux(k,j,i) = kappa * vel_zp / ( & 757 760 LOG( zp / z0(j,i) ) - & 758 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&759 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&760 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&761 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 762 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +& 763 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 761 764 ) 762 765 ENDIF … … 836 839 837 840 838 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )839 wall_flux = 0.0 841 zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx ) 842 wall_flux = 0.0_wp 840 843 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 841 844 … … 851 854 !-- All subsequent variables are computed for scalar locations 852 855 IF ( k >= nzb_diff_s_inner(j,i)-1 .AND. & 853 k <= nzb_diff_s_outer(j,i)-2 .AND. wall(j,i) /= 0.0 ) THEN 856 k <= nzb_diff_s_outer(j,i)-2 .AND. & 857 wall(j,i) /= 0.0_wp ) THEN 854 858 ! 855 859 !-- (1) Compute rifs, u_i, v_i, and ws … … 859 863 kk = k-1 860 864 ENDIF 861 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 862 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 863 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 864 ) 865 866 u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) 867 v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) 868 ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) ) 865 rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + & 866 a * rif_wall(k,j,i+1,1) + & 867 b * rif_wall(k,j+1,i,2) + & 868 c1 * rif_wall(kk,j,i,3) + & 869 c2 * rif_wall(kk,j,i,4) & 870 ) 871 872 u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 873 v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 874 ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 869 875 ! 870 876 !-- (2) Compute wall-parallel absolute velocity vel_total and 871 877 !-- interpolate appropriate velocity component vel_zp. 872 878 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 873 vel_zp = 0.5* ( a * u_i + b * v_i + (c1+c2) * ws )879 vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws ) 874 880 ! 875 881 !-- (3) Compute wall friction velocity us_wall 876 IF ( rifs >= 0.0 ) THEN882 IF ( rifs >= 0.0_wp ) THEN 877 883 878 884 ! 879 885 !-- Stable stratification (and neutral) 880 886 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 881 5.0* rifs * ( zp - z0(j,i) ) / zp &887 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 882 888 ) 883 889 ELSE … … 885 891 ! 886 892 !-- Unstable stratification 887 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )888 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )893 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 894 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 889 895 890 896 us_wall = kappa * vel_total / ( & 891 897 LOG( zp / z0(j,i) ) - & 892 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&893 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&894 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&898 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 899 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +& 900 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 895 901 ) 896 902 ENDIF … … 902 908 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 903 909 904 IF ( rifs >= 0.0 ) THEN910 IF ( rifs >= 0.0_wp ) THEN 905 911 906 912 ! 907 913 !-- Stable stratification (and neutral) 908 wall_flux(k,j,i) = kappa * vel_zp / & 909 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp ) 914 wall_flux(k,j,i) = kappa * vel_zp / ( & 915 LOG( zp/z0(j,i) ) + & 916 5.0_wp * rifs * ( zp-z0(j,i) ) / zp & 917 ) 910 918 ELSE 911 919 912 920 ! 913 921 !-- Unstable stratification 914 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )915 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )922 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 923 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 916 924 917 925 wall_flux(k,j,i) = kappa * vel_zp / ( & 918 926 LOG( zp / z0(j,i) ) - & 919 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&920 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&921 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&927 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 928 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +& 929 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 922 930 ) 923 931 ENDIF … … 981 989 982 990 983 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )984 wall_flux = 0.0 991 zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx ) 992 wall_flux = 0.0_wp 985 993 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 986 994 … … 996 1004 kk = k-1 997 1005 ENDIF 998 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 999 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 1000 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 1001 ) 1002 1003 u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) 1004 v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) 1005 ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) ) 1006 rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + & 1007 a * rif_wall(k,j,i+1,1) + & 1008 b * rif_wall(k,j+1,i,2) + & 1009 c1 * rif_wall(kk,j,i,3) + & 1010 c2 * rif_wall(kk,j,i,4) & 1011 ) 1012 1013 u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 1014 v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 1015 ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 1006 1016 ! 1007 1017 !-- (2) Compute wall-parallel absolute velocity vel_total and 1008 1018 !-- interpolate appropriate velocity component vel_zp. 1009 1019 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 1010 vel_zp = 0.5* ( a * u_i + b * v_i + (c1+c2) * ws )1020 vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws ) 1011 1021 ! 1012 1022 !-- (3) Compute wall friction velocity us_wall 1013 IF ( rifs >= 0.0 ) THEN1023 IF ( rifs >= 0.0_wp ) THEN 1014 1024 1015 1025 ! 1016 1026 !-- Stable stratification (and neutral) 1017 1027 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 1018 5.0* rifs * ( zp - z0(j,i) ) / zp &1028 5.0_wp * rifs * ( zp - z0(j,i) ) / zp & 1019 1029 ) 1020 1030 ELSE … … 1022 1032 ! 1023 1033 !-- Unstable stratification 1024 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )1025 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )1034 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 1035 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 1026 1036 1027 1037 us_wall = kappa * vel_total / ( & 1028 1038 LOG( zp / z0(j,i) ) - & 1029 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&1030 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&1031 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&1039 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 1040 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) + & 1041 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 1032 1042 ) 1033 1043 ENDIF … … 1041 1051 !-- subroutine wall_fluxes because fluxes in subroutine 1042 1052 !-- wall_fluxes_e are defined at scalar locations). 1043 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) +&1044 b * ( v(k,j,i) + v(k,j+1,i) ) +&1045 (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )&1046 )1047 1048 IF ( rifs >= 0.0 ) THEN1053 vel_zp = 0.5_wp * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & 1054 b * ( v(k,j,i) + v(k,j+1,i) ) + & 1055 (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & 1056 ) 1057 1058 IF ( rifs >= 0.0_wp ) THEN 1049 1059 1050 1060 ! 1051 1061 !-- Stable stratification (and neutral) 1052 1062 wall_flux(k) = kappa * vel_zp / & 1053 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )1063 ( LOG( zp/z0(j,i) ) + 5.0_wp * rifs * ( zp-z0(j,i) ) / zp ) 1054 1064 ELSE 1055 1065 1056 1066 ! 1057 1067 !-- Unstable stratification 1058 h1 = SQRT( SQRT( 1.0 - 16.0* rifs ) )1059 h2 = SQRT( SQRT( 1.0 - 16.0* rifs * z0(j,i) / zp ) )1068 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) ) 1069 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) ) 1060 1070 1061 1071 wall_flux(k) = kappa * vel_zp / ( & 1062 1072 LOG( zp / z0(j,i) ) - & 1063 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (&1064 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) +&1065 2.0 * ( ATAN( h1 ) - ATAN( h2 ) )&1066 1073 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( & 1074 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) + & 1075 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) & 1076 ) 1067 1077 ENDIF 1068 1078 wall_flux(k) = - wall_flux(k) * us_wall
Note: See TracChangeset
for help on using the changeset viewer.
