cubically interpolate in t
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer | :: | nx | ||||
integer | :: | ny | ||||
integer | :: | nz | ||||
integer | :: | nt | ||||
real(kind=wp) | :: | p(nx,ny,nz,nt) | ||||
integer | :: | mx | ||||
integer | :: | my | ||||
integer | :: | mz | ||||
integer | :: | mt | ||||
integer | :: | mxmy | ||||
integer | :: | mxmymz | ||||
real(kind=wp) | :: | q(mxmymz,mt) | ||||
integer | :: | intpol(4) | ||||
integer | :: | lt(mt) | ||||
real(kind=wp) | :: | dtm(mt) | ||||
real(kind=wp) | :: | dt(mt) | ||||
real(kind=wp) | :: | dtp(mt) | ||||
real(kind=wp) | :: | dtpp(mt) | ||||
real(kind=wp) | :: | ptm(mxmymz) | ||||
real(kind=wp) | :: | pt(mxmymz) | ||||
real(kind=wp) | :: | ptp(mxmymz) | ||||
real(kind=wp) | :: | ptpp(mxmymz) | ||||
integer | :: | kz(mz) | ||||
real(kind=wp) | :: | dzm(mz) | ||||
real(kind=wp) | :: | dz(mz) | ||||
real(kind=wp) | :: | dzp(mz) | ||||
real(kind=wp) | :: | dzpp(mz) | ||||
real(kind=wp) | :: | pkm(mxmy) | ||||
real(kind=wp) | :: | pk(mxmy) | ||||
real(kind=wp) | :: | pkp(mxmy) | ||||
real(kind=wp) | :: | pkpp(mxmy) | ||||
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) |
subroutine cubt4(nx,ny,nz,nt,p,mx,my,mz,mt,mxmy,mxmymz,q,intpol,& lt,dtm,dt,dtp,dtpp,ptm,pt,ptp,ptpp,& kz,dzm,dz,dzp,dzpp,pkm,pk,pkp,pkpp,& jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) implicit none integer :: nx,ny,nz,nt,mx,my,mz,mt,mxmy,mxmymz,lsave,ll,l,iijjkk integer :: lt(mt),kz(mz),jy(my),ix(mx),intpol(4) real(wp) :: p(nx,ny,nz,nt),q(mxmymz,mt) real(wp) :: dtm(mt),dt(mt),dtp(mt),dtpp(mt) real(wp) :: ptm(mxmymz),pt(mxmymz),ptp(mxmymz),ptpp(mxmymz) real(wp) :: dzm(mz),dz(mz),dzp(mz),dzpp(mz) real(wp) :: pkm(mxmy),pk(mxmy),pkp(mxmy),pkpp(mxmy) 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(3) == 1) then ! linear in z lsave = -3 do ll=1,mt l = lt(ll) if (l==lsave) then ! l pointer has not moved since last pass (no updates or interpolation) else if (l==lsave+1) then ! update l-1,l,l+1 and interpolate l+2 do iijjkk=1,mxmymz ptm(iijjkk) = pt(iijjkk) pt(iijjkk) = ptp(iijjkk) ptp(iijjkk) = ptpp(iijjkk) end do call lint3(nx,ny,nz,p(1,1,1,l+2),mx,my,mxmy,mz,ptpp,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) else if (l==lsave+2) then ! update l-1,l and interpolate l+1,l+2 do iijjkk=1,mxmymz ptm(iijjkk) = ptp(iijjkk) pt(iijjkk) = ptpp(iijjkk) end do call lint3(nx,ny,nz,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) call lint3(nx,ny,nz,p(1,1,1,l+2),mx,my,mxmy,mz,ptpp,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) else if (l==lsave+3) then ! update l-1 and interpolate l,l+1,l+2 do iijjkk=1,mxmymz ptm(iijjkk) = ptpp(iijjkk) end do call lint3(nx,ny,nz,p(1,1,1,l),mx,my,mxmy,mz,pt,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) call lint3(nx,ny,nz,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) call lint3(nx,ny,nz,p(1,1,1,l+2),mx,my,mxmy,mz,ptpp,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) else ! interpolate all four l-1,l,l+1,l+2 call lint3(nx,ny,nz,p(1,1,1,l-1),mx,my,mxmy,mz,ptm,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) call lint3(nx,ny,nz,p(1,1,1,l),mx,my,mxmy,mz,pt,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) call lint3(nx,ny,nz,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) call lint3(nx,ny,nz,p(1,1,1,l+2),mx,my,mxmy,mz,ptpp,intpol,kz,dz,& pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) end if ! save l pointer for next pass lsave = l ! cubically interpolate q(ii,jj,kk,ll) from ptm,pt,ptp,ptpp in t direction do iijjkk=1,mxmymz q(iijjkk,ll) = dtm(ll)*ptm(iijjkk) + dt(ll)*pt(iijjkk) + dtp(ll)*ptp(iijjkk) + dtpp(ll)*ptpp(iijjkk) end do end do else ! cubic in z lsave = -3 do ll=1,mt l = lt(ll) if (l==lsave) then ! l pointer has not moved since last pass (no updates or interpolation) else if (l==lsave+1) then ! update l-1,l,l+1 and interpolate l+2 do iijjkk=1,mxmymz ptm(iijjkk) = pt(iijjkk) pt(iijjkk) = ptp(iijjkk) ptp(iijjkk) = ptpp(iijjkk) end do call cubt3(nx,ny,nz,p(1,1,1,l+2),mx,my,mxmy,mz,ptpp,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) else if (l==lsave+2) then ! update l-1,l and interpolate l+1,l+2 do iijjkk=1,mxmymz ptm(iijjkk) = ptp(iijjkk) pt(iijjkk) = ptpp(iijjkk) end do call cubt3(nx,ny,nz,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) call cubt3(nx,ny,nz,p(1,1,1,l+2),mx,my,mxmy,mz,ptpp,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) else if (l==lsave+3) then ! update l-1 and interpolate l,l+1,l+2 do iijjkk=1,mxmymz ptm(iijjkk) = ptpp(iijjkk) end do call cubt3(nx,ny,nz,p(1,1,1,l),mx,my,mxmy,mz,pt,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) call cubt3(nx,ny,nz,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) call cubt3(nx,ny,nz,p(1,1,1,l+2),mx,my,mxmy,mz,ptpp,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) else ! interpolate all four l-1,l,l+1,l+2 call cubt3(nx,ny,nz,p(1,1,1,l-1),mx,my,mxmy,mz,ptm,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) call cubt3(nx,ny,nz,p(1,1,1,l),mx,my,mxmy,mz,pt,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) call cubt3(nx,ny,nz,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) call cubt3(nx,ny,nz,p(1,1,1,l+2),mx,my,mxmy,mz,ptpp,intpol,kz,dzm,& dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,& ix,dxm,dx,dxp,dxpp) end if ! save l pointer for next pass lsave = l ! cubically interpolate q(ii,jj,kk,ll) from ptm,pt,ptp,ptpp in t direction do iijjkk=1,mxmymz q(iijjkk,ll) = dtm(ll)*ptm(iijjkk) + dt(ll)*pt(iijjkk) + dtp(ll)*ptp(iijjkk) + dtpp(ll)*ptpp(iijjkk) end do end do end if end subroutine cubt4