Changeset 1015 for palm/trunk/SOURCE/wall_fluxes.f90
- Timestamp:
- Sep 27, 2012 9:23:24 AM (12 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wall_fluxes.f90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wall_fluxes.f90
r667 r1015 3 3 ! Current revisions: 4 4 ! ----------------- 5 ! 5 ! accelerator version (*_acc) added 6 6 ! 7 7 ! Former revisions: … … 32 32 !------------------------------------------------------------------------------! 33 33 PRIVATE 34 PUBLIC wall_fluxes, wall_fluxes_ e34 PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc 35 35 36 36 INTERFACE wall_fluxes … … 39 39 END INTERFACE wall_fluxes 40 40 41 INTERFACE wall_fluxes_acc 42 MODULE PROCEDURE wall_fluxes_acc 43 END INTERFACE wall_fluxes_acc 44 41 45 INTERFACE wall_fluxes_e 42 46 MODULE PROCEDURE wall_fluxes_e … … 44 48 END INTERFACE wall_fluxes_e 45 49 50 INTERFACE wall_fluxes_e_acc 51 MODULE PROCEDURE wall_fluxes_e_acc 52 END INTERFACE wall_fluxes_e_acc 53 46 54 CONTAINS 47 55 … … 189 197 END SUBROUTINE wall_fluxes 190 198 199 200 !------------------------------------------------------------------------------! 201 ! Call for all grid points - accelerator version 202 !------------------------------------------------------------------------------! 203 SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, & 204 nzb_uvw_outer, wall ) 205 206 USE arrays_3d 207 USE control_parameters 208 USE grid_variables 209 USE indices 210 USE statistics 211 212 IMPLICIT NONE 213 214 INTEGER :: i, j, k, max_outer, min_inner, wall_index 215 216 INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) :: nzb_uvw_inner, & 217 nzb_uvw_outer 218 REAL :: a, b, c1, c2, h1, h2, zp 219 REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts 220 221 REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: wall 222 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wall_flux 223 224 225 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 226 wall_flux = 0.0 227 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 228 229 min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1 230 max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) ) 231 232 !$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) & 233 !$acc present( u, v, w, wall, wall_flux, z0 ) 234 !$acc loop 235 DO i = nxl, nxr 236 DO j = nys, nyn 237 !$acc loop vector( 32 ) 238 DO k = min_inner, max_outer 239 ! 240 !-- All subsequent variables are computed for the respective 241 !-- location where the respective flux is defined. 242 IF ( k >= nzb_uvw_inner(j,i)+1 .AND. & 243 k <= nzb_uvw_outer(j,i) .AND. wall(j,i) /= 0.0 ) THEN 244 ! 245 !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' 246 rifs = rif_wall(k,j,i,wall_index) 247 248 u_i = a * u(k,j,i) + c1 * 0.25 * & 249 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 250 251 v_i = b * v(k,j,i) + c2 * 0.25 * & 252 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 253 254 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( & 255 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 256 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 257 ) 258 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + & 259 b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) ) 260 261 pts = pt_i - hom(k,1,4,0) 262 wspts = ws * pts 263 264 ! 265 !-- (2) Compute wall-parallel absolute velocity vel_total 266 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 267 268 ! 269 !-- (3) Compute wall friction velocity us_wall 270 IF ( rifs >= 0.0 ) THEN 271 272 ! 273 !-- Stable stratification (and neutral) 274 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 275 5.0 * rifs * ( zp - z0(j,i) ) / zp & 276 ) 277 ELSE 278 279 ! 280 !-- Unstable stratification 281 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 282 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 283 284 us_wall = kappa * vel_total / ( & 285 LOG( zp / z0(j,i) ) - & 286 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 287 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 288 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 289 ) 290 ENDIF 291 292 ! 293 !-- (4) Compute zp/L (corresponds to neutral Richardson flux 294 !-- number rifs) 295 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * & 296 ( us_wall**3 + 1E-30 ) ) 297 298 ! 299 !-- Limit the value range of the Richardson numbers. 300 !-- This is necessary for very small velocities (u,w --> 0), 301 !-- because the absolute value of rif can then become very 302 !-- large, which in consequence would result in very large 303 !-- shear stresses and very small momentum fluxes (both are 304 !-- generally unrealistic). 305 IF ( rifs < rif_min ) rifs = rif_min 306 IF ( rifs > rif_max ) rifs = rif_max 307 308 ! 309 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 310 IF ( rifs >= 0.0 ) THEN 311 312 ! 313 !-- Stable stratification (and neutral) 314 wall_flux(k,j,i) = kappa * & 315 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 316 ( LOG( zp / z0(j,i) ) + & 317 5.0 * rifs * ( zp - z0(j,i) ) / zp & 318 ) 319 ELSE 320 321 ! 322 !-- Unstable stratification 323 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 324 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 325 326 wall_flux(k,j,i) = kappa * & 327 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / ( & 328 LOG( zp / z0(j,i) ) - & 329 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 330 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 331 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 332 ) 333 ENDIF 334 wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall 335 336 ! 337 !-- store rifs for next time step 338 rif_wall(k,j,i,wall_index) = rifs 339 340 ENDIF 341 342 ENDDO 343 ENDDO 344 ENDDO 345 !$acc end kernels 346 347 END SUBROUTINE wall_fluxes_acc 191 348 192 349 … … 444 601 END SUBROUTINE wall_fluxes_e 445 602 603 604 !------------------------------------------------------------------------------! 605 ! Call for all grid points - accelerator version 606 !------------------------------------------------------------------------------! 607 SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall ) 608 609 !------------------------------------------------------------------------------! 610 ! Description: 611 ! ------------ 612 ! Calculates momentum fluxes at vertical walls for routine production_e 613 ! assuming Monin-Obukhov similarity. 614 ! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0). 615 !------------------------------------------------------------------------------! 616 617 USE arrays_3d 618 USE control_parameters 619 USE grid_variables 620 USE indices 621 USE statistics 622 623 IMPLICIT NONE 624 625 INTEGER :: i, j, k, kk, max_outer, min_inner, wall_index 626 REAL :: a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, & 627 ws, zp 628 629 REAL :: rifs 630 631 REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: wall 632 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wall_flux 633 634 635 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 636 wall_flux = 0.0 637 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 638 639 min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1 640 max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2 641 642 !$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, pt, rif_wall ) & 643 !$acc present( u, v, w, wall, wall_flux, z0 ) 644 !$acc loop 645 DO i = nxl, nxr 646 DO j = nys, nyn 647 !$acc loop vector(32) 648 DO k = min_inner, max_outer 649 ! 650 !-- All subsequent variables are computed for scalar locations 651 IF ( k >= nzb_diff_s_inner(j,i)-1 .AND. & 652 k <= nzb_diff_s_outer(j,i)-2 .AND. wall(j,i) /= 0.0 ) THEN 653 ! 654 !-- (1) Compute rifs, u_i, v_i, and ws 655 IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN 656 kk = nzb_diff_s_inner(j,i)-1 657 ELSE 658 kk = k-1 659 ENDIF 660 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 661 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 662 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 663 ) 664 665 u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) 666 v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) 667 ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) ) 668 ! 669 !-- (2) Compute wall-parallel absolute velocity vel_total and 670 !-- interpolate appropriate velocity component vel_zp. 671 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 672 vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws ) 673 ! 674 !-- (3) Compute wall friction velocity us_wall 675 IF ( rifs >= 0.0 ) THEN 676 677 ! 678 !-- Stable stratification (and neutral) 679 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 680 5.0 * rifs * ( zp - z0(j,i) ) / zp & 681 ) 682 ELSE 683 684 ! 685 !-- Unstable stratification 686 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 687 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 688 689 us_wall = kappa * vel_total / ( & 690 LOG( zp / z0(j,i) ) - & 691 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 692 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 693 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 694 ) 695 ENDIF 696 697 ! 698 !-- Skip step (4) of wall_fluxes, because here rifs is already 699 !-- available from (1) 700 ! 701 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 702 703 IF ( rifs >= 0.0 ) THEN 704 705 ! 706 !-- Stable stratification (and neutral) 707 wall_flux(k,j,i) = kappa * vel_zp / & 708 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp ) 709 ELSE 710 711 ! 712 !-- Unstable stratification 713 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 714 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 715 716 wall_flux(k,j,i) = kappa * vel_zp / ( & 717 LOG( zp / z0(j,i) ) - & 718 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 719 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 720 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 721 ) 722 ENDIF 723 wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall 724 725 ENDIF 726 727 ENDDO 728 ENDDO 729 ENDDO 730 !$acc end kernels 731 732 END SUBROUTINE wall_fluxes_e_acc 446 733 447 734
Note: See TracChangeset
for help on using the changeset viewer.
