cubically interpolate p onto q in y
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer | :: | nx | ||||
integer | :: | ny | ||||
real(kind=wp) | :: | p(nx,ny) | ||||
integer | :: | mx | ||||
integer | :: | my | ||||
real(kind=wp) | :: | q(mx,my) | ||||
integer | :: | intpol(2) | ||||
integer | :: | jy(my) | ||||
real(kind=wp) | :: | dym(my) | ||||
real(kind=wp) | :: | dy(my) | ||||
real(kind=wp) | :: | dyp(my) | ||||
real(kind=wp) | :: | dypp(my) | ||||
real(kind=wp) | :: | pjm(mx) | ||||
real(kind=wp) | :: | pj(mx) | ||||
real(kind=wp) | :: | pjp(mx) | ||||
real(kind=wp) | :: | pjpp(mx) | ||||
integer | :: | ix(mx) | ||||
real(kind=wp) | :: | dxm(mx) | ||||
real(kind=wp) | :: | dx(mx) | ||||
real(kind=wp) | :: | dxp(mx) | ||||
real(kind=wp) | :: | dxpp(mx) | ||||
integer | :: | inmx | ||||
integer | :: | jnmy | ||||
integer | :: | isubx | ||||
integer | :: | jsuby |
subroutine cubt2u(nx,ny,p,mx,my,q,intpol,jy,dym,dy,dyp,dypp,& pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp,inmx,jnmy,& isubx,jsuby) implicit none integer :: nx,ny,mx,my,intpol(2),jy(my),ix(mx),inmx,jnmy,isubx,jsuby integer :: j,jj,ii,jsave real(wp) :: p(nx,ny),q(mx,my) real(wp) :: dym(my),dy(my),dyp(my),dypp(my) real(wp) :: pjm(mx),pj(mx),pjp(mx),pjpp(mx) real(wp) :: dxm(mx),dx(mx),dxp(mx),dxpp(mx) if (intpol(1) == 1) then ! linear in x if (jsuby == 1) then ! my grid is subset of ny grid do jj=1,my j = jnmy*(jj-1)+1 call lint1u(nx,p(1,j),mx,q(1,jj),ix,dx,inmx,isubx) end do return end if jsave = -3 do jj=1,my j = jy(jj) ! load pjm,pj,pjp,pjpp if (j==jsave) then ! no updates or x interpolation necessary else if (j==jsave+1) then do ii=1,mx pjm(ii) = pj(ii) pj(ii) = pjp(ii) pjp(ii) = pjpp(ii) end do call lint1u(nx,p(1,j+2),mx,pjpp,ix,dx,inmx,isubx) else if (j==jsave+2) then do ii=1,mx pjm(ii) = pjp(ii) pj(ii) = pjpp(ii) end do call lint1u(nx,p(1,j+1),mx,pjp,ix,dx,inmx,isubx) call lint1u(nx,p(1,j+2),mx,pjpp,ix,dx,inmx,isubx) else if (j==jsave+3) then do ii=1,mx pjm(ii) = pjpp(ii) end do call lint1u(nx,p(1,j),mx,pj,ix,dx,inmx,isubx) call lint1u(nx,p(1,j+1),mx,pjp,ix,dx,inmx,isubx) call lint1u(nx,p(1,j+2),mx,pjpp,ix,dx,inmx,isubx) else ! load all four (no updates) call lint1u(nx,p(1,j-1),mx,pjm,ix,dx,inmx,isubx) call lint1u(nx,p(1,j),mx,pj,ix,dx,inmx,isubx) call lint1u(nx,p(1,j+1),mx,pjp,ix,dx,inmx,isubx) call lint1u(nx,p(1,j+2),mx,pjpp,ix,dx,inmx,isubx) end if ! update pointer jsave = j do ii=1,mx q(ii,jj) = dym(jj)*pjm(ii) + dy(jj)*pj(ii) + dyp(jj)*pjp(ii) + dypp(jj)*pjpp(ii) end do end do else ! cubic in x if (jsuby == 1) then ! my grid is subset of ny grid do jj=1,my j = jnmy*(jj-1)+1 call cubt1u(nx,p(1,j),mx,q(1,jj),ix,dxm,dx,dxp,dxpp,inmx,isubx) end do return end if jsave = -3 do jj=1,my j = jy(jj) ! load pjm,pj,pjp,pjpp if (j==jsave) then ! no updates or x interpolation necessary else if (j==jsave+1) then do ii=1,mx pjm(ii) = pj(ii) pj(ii) = pjp(ii) pjp(ii) = pjpp(ii) end do call cubt1u(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp,inmx,isubx) else if (j==jsave+2) then do ii=1,mx pjm(ii) = pjp(ii) pj(ii) = pjpp(ii) end do call cubt1u(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp,inmx,isubx) call cubt1u(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp,inmx,isubx) else if (j==jsave+3) then do ii=1,mx pjm(ii) = pjpp(ii) end do call cubt1u(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp,inmx,isubx) call cubt1u(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp,inmx,isubx) call cubt1u(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp,inmx,isubx) else ! load all four (no updates) call cubt1u(nx,p(1,j-1),mx,pjm,ix,dxm,dx,dxp,dxpp,inmx,isubx) call cubt1u(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp,inmx,isubx) call cubt1u(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp,inmx,isubx) call cubt1u(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp,inmx,isubx) end if ! update pointer jsave = j do ii=1,mx q(ii,jj) = dym(jj)*pjm(ii) + dy(jj)*pj(ii) + dyp(jj)*pjp(ii) + dypp(jj)*pjpp(ii) end do end do return end if end subroutine cubt2u