Changeset 2118 for palm/trunk/SOURCE/wall_fluxes.f90
- Timestamp:
- Jan 17, 2017 4:38:49 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wall_fluxes.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wall_fluxes.f90
r2101 r2118 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! OpenACC versions of subroutines removed 23 23 ! 24 24 ! Former revisions: … … 90 90 91 91 PRIVATE 92 PUBLIC wall_fluxes, wall_fluxes_ acc, wall_fluxes_e, wall_fluxes_e_acc92 PUBLIC wall_fluxes, wall_fluxes_e 93 93 94 94 INTERFACE wall_fluxes … … 97 97 END INTERFACE wall_fluxes 98 98 99 INTERFACE wall_fluxes_acc100 MODULE PROCEDURE wall_fluxes_acc101 END INTERFACE wall_fluxes_acc102 103 99 INTERFACE wall_fluxes_e 104 100 MODULE PROCEDURE wall_fluxes_e … … 106 102 END INTERFACE wall_fluxes_e 107 103 108 INTERFACE wall_fluxes_e_acc109 MODULE PROCEDURE wall_fluxes_e_acc110 END INTERFACE wall_fluxes_e_acc111 112 104 CONTAINS 113 105 … … 299 291 ! Description: 300 292 ! ------------ 301 !> Call for all grid points - accelerator version 302 !------------------------------------------------------------------------------! 303 SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, & 304 nzb_uvw_outer, wall ) 293 !> Call for all grid point i,j 294 !------------------------------------------------------------------------------! 295 SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) 305 296 306 297 USE arrays_3d, & … … 314 305 315 306 USE indices, & 316 ONLY: i_left, i_right, j_north, j_south, nxl, nxlg, nxr, nxrg, & 317 nyn, nyng, nys, nysg, nzb, nzt 307 ONLY: nzb, nzt 318 308 319 309 USE kinds … … 327 317 INTEGER(iwp) :: j !< 328 318 INTEGER(iwp) :: k !< 329 INTEGER(iwp) :: max_outer!<330 INTEGER(iwp) :: min_inner!<319 INTEGER(iwp) :: nzb_w !< 320 INTEGER(iwp) :: nzt_w !< 331 321 INTEGER(iwp) :: wall_index !< 332 333 INTEGER(iwp), &334 DIMENSION(nysg:nyng,nxlg:nxrg) :: &335 nzb_uvw_inner !<336 INTEGER(iwp), &337 DIMENSION(nysg:nyng,nxlg:nxrg) :: &338 nzb_uvw_outer !<339 322 340 323 REAL(wp) :: a !< … … 355 338 REAL(wp) :: wspts !< 356 339 357 REAL(wp), &358 DIMENSION(nysg:nyng,nxlg:nxrg) :: &359 wall !<360 361 REAL(wp), &362 DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &363 wall_flux !<364 365 366 zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )367 wall_flux = 0.0_wp368 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )369 370 min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1371 max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) )372 373 !$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) &374 !$acc present( u, v, w, wall, wall_flux, z0 )375 !$acc loop independent376 DO i = i_left, i_right377 DO j = j_south, j_north378 379 IF ( wall(j,i) /= 0.0_wp ) THEN380 !381 !-- All subsequent variables are computed for the respective382 !-- location where the respective flux is defined.383 !$acc loop independent384 DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)385 386 !387 !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'388 rifs = rif_wall(k,j,i,wall_index)389 390 u_i = a * u(k,j,i) + c1 * 0.25_wp * &391 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )392 393 v_i = b * v(k,j,i) + c2 * 0.25_wp * &394 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )395 396 ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &397 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &398 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &399 )400 pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + &401 b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )402 403 pts = pt_i - hom(k,1,4,0)404 wspts = ws * pts405 406 !407 !-- (2) Compute wall-parallel absolute velocity vel_total408 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )409 410 !411 !-- (3) Compute wall friction velocity us_wall412 IF ( rifs >= 0.0_wp ) THEN413 414 !415 !-- Stable stratification (and neutral)416 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &417 5.0_wp * rifs * ( zp - z0(j,i) ) / zp &418 )419 ELSE420 421 !422 !-- Unstable stratification423 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )424 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )425 426 us_wall = kappa * vel_total / ( &427 LOG( zp / z0(j,i) ) - &428 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( &429 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +&430 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &431 )432 ENDIF433 434 !435 !-- (4) Compute zp/L (corresponds to neutral Richardson flux436 !-- number rifs)437 rifs = -1.0_wp * zp * kappa * g * wspts / &438 ( pt_i * ( us_wall**3 + 1E-30 ) )439 440 !441 !-- Limit the value range of the Richardson numbers.442 !-- This is necessary for very small velocities (u,w --> 0),443 !-- because the absolute value of rif can then become very444 !-- large, which in consequence would result in very large445 !-- shear stresses and very small momentum fluxes (both are446 !-- generally unrealistic).447 IF ( rifs < zeta_min ) rifs = zeta_min448 IF ( rifs > zeta_max ) rifs = zeta_max449 450 !451 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')452 IF ( rifs >= 0.0_wp ) THEN453 454 !455 !-- Stable stratification (and neutral)456 wall_flux(k,j,i) = kappa * &457 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &458 ( LOG( zp / z0(j,i) ) + &459 5.0_wp * rifs * ( zp - z0(j,i) ) / zp &460 )461 ELSE462 463 !464 !-- Unstable stratification465 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )466 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )467 468 wall_flux(k,j,i) = kappa * &469 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &470 LOG( zp / z0(j,i) ) - &471 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( &472 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +&473 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &474 )475 ENDIF476 wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall477 478 !479 !-- store rifs for next time step480 rif_wall(k,j,i,wall_index) = rifs481 482 ENDDO483 484 ENDIF485 486 ENDDO487 ENDDO488 !$acc end kernels489 490 END SUBROUTINE wall_fluxes_acc491 492 493 !------------------------------------------------------------------------------!494 ! Description:495 ! ------------496 !> Call for all grid point i,j497 !------------------------------------------------------------------------------!498 SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )499 500 USE arrays_3d, &501 ONLY: rif_wall, pt, u, v, w, z0502 503 USE control_parameters, &504 ONLY: g, kappa, zeta_max, zeta_min505 506 USE grid_variables, &507 ONLY: dx, dy508 509 USE indices, &510 ONLY: nzb, nzt511 512 USE kinds513 514 USE statistics, &515 ONLY: hom516 517 IMPLICIT NONE518 519 INTEGER(iwp) :: i !<520 INTEGER(iwp) :: j !<521 INTEGER(iwp) :: k !<522 INTEGER(iwp) :: nzb_w !<523 INTEGER(iwp) :: nzt_w !<524 INTEGER(iwp) :: wall_index !<525 526 REAL(wp) :: a !<527 REAL(wp) :: b !<528 REAL(wp) :: c1 !<529 REAL(wp) :: c2 !<530 REAL(wp) :: h1 !<531 REAL(wp) :: h2 !<532 REAL(wp) :: zp !<533 REAL(wp) :: pts !<534 REAL(wp) :: pt_i !<535 REAL(wp) :: rifs !<536 REAL(wp) :: u_i !<537 REAL(wp) :: v_i !<538 REAL(wp) :: us_wall !<539 REAL(wp) :: vel_total !<540 REAL(wp) :: ws !<541 REAL(wp) :: wspts !<542 543 340 REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !< 544 341 … … 811 608 ! Description: 812 609 ! ------------ 813 !> Call for all grid points - accelerator version 814 !> Calculates momentum fluxes at vertical walls for routine production_e 815 !> assuming Monin-Obukhov similarity. 816 !> Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0). 817 !------------------------------------------------------------------------------! 818 SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall ) 819 610 !> Call for grid point i,j 611 !------------------------------------------------------------------------------! 612 SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) 820 613 821 614 USE arrays_3d, & … … 829 622 830 623 USE indices, & 831 ONLY: i_left, i_right, j_north, j_south, nxl, nxlg, nxr, nxrg, & 832 nyn, nyng, nys, nysg, nzb, nzb_diff_s_inner, & 833 nzb_diff_s_outer, nzt 624 ONLY: nzb, nzt 834 625 835 626 USE kinds … … 841 632 INTEGER(iwp) :: k !< 842 633 INTEGER(iwp) :: kk !< 843 INTEGER(iwp) :: max_outer!<844 INTEGER(iwp) :: min_inner!<634 INTEGER(iwp) :: nzb_w !< 635 INTEGER(iwp) :: nzt_w !< 845 636 INTEGER(iwp) :: wall_index !< 846 637 … … 860 651 REAL(wp) :: rifs !< 861 652 862 REAL(wp), &863 DIMENSION(nysg:nyng,nxlg:nxrg) :: &864 wall !<865 866 REAL(wp), &867 DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &868 wall_flux !<869 870 871 zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )872 wall_flux = 0.0_wp873 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )874 875 min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1876 max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2877 878 !$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, rif_wall ) &879 !$acc present( u, v, w, wall, wall_flux, z0 )880 DO i = i_left, i_right881 DO j = j_south, j_north882 DO k = min_inner, max_outer883 !884 !-- All subsequent variables are computed for scalar locations885 IF ( k >= nzb_diff_s_inner(j,i)-1 .AND. &886 k <= nzb_diff_s_outer(j,i)-2 .AND. &887 wall(j,i) /= 0.0_wp ) THEN888 !889 !-- (1) Compute rifs, u_i, v_i, and ws890 IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN891 kk = nzb_diff_s_inner(j,i)-1892 ELSE893 kk = k-1894 ENDIF895 rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + &896 a * rif_wall(k,j,i+1,1) + &897 b * rif_wall(k,j+1,i,2) + &898 c1 * rif_wall(kk,j,i,3) + &899 c2 * rif_wall(kk,j,i,4) &900 )901 902 u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )903 v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )904 ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )905 !906 !-- (2) Compute wall-parallel absolute velocity vel_total and907 !-- interpolate appropriate velocity component vel_zp.908 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )909 vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )910 !911 !-- (3) Compute wall friction velocity us_wall912 IF ( rifs >= 0.0_wp ) THEN913 914 !915 !-- Stable stratification (and neutral)916 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &917 5.0_wp * rifs * ( zp - z0(j,i) ) / zp &918 )919 ELSE920 921 !922 !-- Unstable stratification923 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )924 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )925 926 us_wall = kappa * vel_total / ( &927 LOG( zp / z0(j,i) ) - &928 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( &929 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +&930 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &931 )932 ENDIF933 934 !935 !-- Skip step (4) of wall_fluxes, because here rifs is already936 !-- available from (1)937 !938 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')939 940 IF ( rifs >= 0.0_wp ) THEN941 942 !943 !-- Stable stratification (and neutral)944 wall_flux(k,j,i) = kappa * vel_zp / ( &945 LOG( zp/z0(j,i) ) + &946 5.0_wp * rifs * ( zp-z0(j,i) ) / zp &947 )948 ELSE949 950 !951 !-- Unstable stratification952 h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )953 h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )954 955 wall_flux(k,j,i) = kappa * vel_zp / ( &956 LOG( zp / z0(j,i) ) - &957 LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / ( &958 ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 ) ) ) +&959 2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &960 )961 ENDIF962 wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall963 964 ENDIF965 966 ENDDO967 ENDDO968 ENDDO969 !$acc end kernels970 971 END SUBROUTINE wall_fluxes_e_acc972 973 974 !------------------------------------------------------------------------------!975 ! Description:976 ! ------------977 !> Call for grid point i,j978 !------------------------------------------------------------------------------!979 SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )980 981 USE arrays_3d, &982 ONLY: rif_wall, u, v, w, z0983 984 USE control_parameters, &985 ONLY: kappa986 987 USE grid_variables, &988 ONLY: dx, dy989 990 USE indices, &991 ONLY: nzb, nzt992 993 USE kinds994 995 IMPLICIT NONE996 997 INTEGER(iwp) :: i !<998 INTEGER(iwp) :: j !<999 INTEGER(iwp) :: k !<1000 INTEGER(iwp) :: kk !<1001 INTEGER(iwp) :: nzb_w !<1002 INTEGER(iwp) :: nzt_w !<1003 INTEGER(iwp) :: wall_index !<1004 1005 REAL(wp) :: a !<1006 REAL(wp) :: b !<1007 REAL(wp) :: c1 !<1008 REAL(wp) :: c2 !<1009 REAL(wp) :: h1 !<1010 REAL(wp) :: h2 !<1011 REAL(wp) :: u_i !<1012 REAL(wp) :: v_i !<1013 REAL(wp) :: us_wall !<1014 REAL(wp) :: vel_total !<1015 REAL(wp) :: vel_zp !<1016 REAL(wp) :: ws !<1017 REAL(wp) :: zp !<1018 REAL(wp) :: rifs !<1019 1020 653 REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !< 1021 654
Note: See TracChangeset
for help on using the changeset viewer.
