subroutine simpin1(wgts,wgtsx1,wgtsx11,iwgts, & iflag,x1in,nin,iord,lbig,x1grid,n1grid,nf, & ndx1,ndx11,ixi,tl,alocal,blocal) !$$$ subprogram documentation block ! . . . . ! subprogram: simpin1 ! prgmmr: parrish org: np22 date: 2000-01-01 ! ! abstract: This routine computes interpolation weights for 1-d ! interpolation ! ! program history log: ! 2000-01-01 parrish ! 2004-06-22 treadon - update documentation ! 2008-04-12 safford - rm unused vars ! ! input argument list: ! x1in - coordinates of interpolatees ! nin - number of interpolatees ! iord - order of interpolation (1=linear, 2=quadratic, etc) ! lbig - number of interpolating points (=iord+1) ! x1grid - coordinates of interpolator grid ! (can be monotonic increasing or decreasing) ! n1grid - dimension of interpolator grid ! nf - =1, then return interpolation weights for function f ! ndx1 - =1, then return interpolation weights for df/dx1 ! ndx11 - =1, then return interpolation weights for d2f/(dx1*dx1) ! ! output argument list: ! wgts - interpolation weights for function ( wgts(nin,lbig) ) ! wgtsx1 - interpolation weights for df/dx1 ! wgtsx11 - interpolation weights for d2f/(dx1*dx1) ! note: if any of these 3 are not asked for, they ! can be dummy arguments ! iwgts - absolute grid addresses ( iwgts(nin,lbig) ) ! iflag - flag for each interpolatee (iflag(nin)) ! =0, then point too close to edge of domain, but weights ! are computed anyway, in case it is considered ! OK to extrapolate. ! =1, then weights computed ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ end documentation block use kinds, only: r_kind,i_kind use constants, only: izero, ione, zero, half, one implicit none integer(i_kind),intent(in ) :: nin,iord,lbig,n1grid,nf,ndx1,ndx11 real(r_kind) ,intent(in ) :: x1in(nin),x1grid(n1grid) real(r_kind) ,intent( out) :: wgts(nin,lbig) integer(i_kind),intent( out) :: iwgts(nin,lbig) real(r_kind) ,intent( out) :: wgtsx1(nin,lbig),wgtsx11(nin,lbig) integer(i_kind),intent(in ) :: ixi(0:iord) real(r_kind) ,intent(in ) :: tl(lbig,lbig,n1grid,2),alocal(n1grid,2),blocal(n1grid,2) real(r_kind) dx1gridi(-3:n1grid+3_i_kind) integer(i_kind) i1ref(nin) integer(i_kind) ix1sign(nin) real(r_kind) x1p(nin) integer(i_kind) ix1temp(nin) integer(i_kind) iflag(nin) real(r_kind) z0(max(64,nin/lbig),lbig) integer(i_kind), allocatable, dimension(:) :: ix1grid real(r_kind) dxa,dxb,dxc,dxd,dxe,dxf,dxg,x1ref,xboundleft,xboundright real(r_kind) dxfine,dxminmin,dxthis,dxfinei,dxmax,dxmin integer(i_kind) i1ref0,i1fine,l,iflip,nn,j,k,ip,nend,nstart,lp,& nthis,nmaxright,n,nminleft,iximx,ntl integer(i_kind) ii,n1fine,iximax,iximn,i ntl=max(64_i_kind,nin/lbig) ! Find and mark all points that are outside interval of interpolation iximx=maxval(ixi) ; iximn=minval(ixi) iximax=max(abs(iximx),abs(iximn)) ! write(6,*)' iximx,mn=',iximx,iximn ! write(6,*)' iximax=',iximax ! write(6,*)' ixi=',ixi nminleft=abs(iximn)+ione ; nmaxright=n1grid-iximax if(iord==ione) then xboundleft=x1grid(1) xboundright=x1grid(n1grid) else xboundleft=x1grid(ione+iximax)-.49999_r_kind*(x1grid(ione+iximax)-x1grid(iximax)) xboundright=x1grid(n1grid-iximax)+ & .49999_r_kind*(x1grid(n1grid-iximax+ione)-x1grid(n1grid-iximax)) end if ! write(6,*) ' xboundleft=',xboundleft ! write(6,*) ' xboundright=',xboundright ! do i=1,iximax+ione ! write(6,*) ' x1grid(',i,')=',x1grid(i) ! end do ! do i=n1grid,n1grid-iximax,-1 ! write(6,*)' x1grid(',i,')=',x1grid(i) ! end do if(x1grid(n1grid)>x1grid(1)) then iflag=ione do n=1,nin if(x1in(n)<=xboundleft.or.x1in(n)>=xboundright) iflag(n)=izero end do end if if(x1grid(n1grid)=xboundleft.or.x1in(n)<=xboundright) iflag(n)=izero end do end if ! Set up uniform fine grid to use in finding interpolation coordinates dxmax=-huge(dxmax) ; dxmin=huge(dxmin) do i=1,n1grid-ione dxthis=x1grid(i+ione)-x1grid(i) dx1gridi(i)=one/dxthis dxmax=max(dxthis,dxmax) ; dxmin=min(dxthis,dxmin) end do dx1gridi(-3:0)=dx1gridi(1) dx1gridi(n1grid:n1grid+3_i_kind)=dx1gridi(n1grid-ione) if(dxmax*dxmin<=zero) then write(6,*)' INTERPOLATION GRID NOT MONOTONIC IN SIMPIN1' stop end if dxminmin=min(abs(dxmax),abs(dxmin)) dxfine=sign(dxminmin,dxmax) dxfinei=one/dxfine n1fine=ione+ceiling((x1grid(n1grid)-x1grid(1))/dxfine) ! write(6,*)' in simpin1, n1grid,n1fine=',n1grid,n1fine allocate (ix1grid(n1fine)) ii=1 do i=1,n1fine x1ref=x1grid(1)+(i-ione)*dxfine if(dxfinei*(x1ref-x1grid(ii+ione))>=-.001_r_kind) ii=min(ii+ione,n1grid-ione) ix1grid(i)=ii ! write(6,*)' x1fine,x1grid=',x1ref,x1grid(ix1grid(i)) end do ! Now get i1ref, index of interval containing point i1ref=-ione do n=1,nin i1fine=ione+nint((x1in(n)-x1grid(1))*dxfinei) i1fine=max(ione,min(i1fine,n1fine)) i1ref0=max(4_i_kind,min(ix1grid(i1fine),n1grid-4_i_kind)) dxa=(x1in(n)-x1grid(i1ref0-3_i_kind))*dx1gridi(i1ref0-3_i_kind) dxb=(x1in(n)-x1grid(i1ref0-2_i_kind))*dx1gridi(i1ref0-2_i_kind) dxc=(x1in(n)-x1grid(i1ref0-ione ))*dx1gridi(i1ref0-ione ) dxd=(x1in(n)-x1grid(i1ref0 ))*dx1gridi(i1ref0 ) dxe=(x1in(n)-x1grid(i1ref0+ione ))*dx1gridi(i1ref0+ione ) dxf=(x1in(n)-x1grid(i1ref0+2_i_kind))*dx1gridi(i1ref0+2_i_kind) dxg=(x1in(n)-x1grid(i1ref0+3_i_kind))*dx1gridi(i1ref0+3_i_kind) if(dxa<=one) i1ref(n)=i1ref0-3_i_kind if(dxb>=zero.and.dxb<=one) i1ref(n)=i1ref0-2_i_kind if(dxc>=zero.and.dxc<=one) i1ref(n)=i1ref0-ione if(dxd>=zero.and.dxd<=one) i1ref(n)=i1ref0 if(dxe>=zero.and.dxe<=one) i1ref(n)=i1ref0+ione if(dxf>=zero.and.dxf<=one) i1ref(n)=i1ref0+2_i_kind if(dxg>=zero) i1ref(n)=i1ref0+3_i_kind i1ref(n)=max(nminleft,min(i1ref(n),nmaxright)) end do deallocate (ix1grid) ! write(6,*)' max,min(i1ref)=',maxval(i1ref),minval(i1ref) ! i1ref now has index of interval containing point ! adjust to be closest one to interpolating point ix1sign=ione do n=1,nin dxa=(x1in(n)-x1grid(i1ref(n)))*dx1gridi(i1ref(n)) if(dxa>half.and.i1ref(n)= 1 and <= ',n1grid ! check that we got all points ! dx1max=-huge(dx1max) ! dx1min=huge(dx1min) ! numinside=izero ! numgtp5=izero ! numoutside=izero ! halfpeps=half+epsilon(x) ! do n=1,nin ! if(iflag(n)/=izero) then ! numinside=numinside+ione ! dx1top=x1in(n)-x1grid(i1ref(n)) ! if(dx1top*dx1gridi(i1ref(n))halfpeps) then ! numgtp5=numgtp5+ione ! end if ! dx1max=max(dx1this,dx1max) ! dx1min=min(dx1this,dx1min) ! else ! numoutside=numoutside+ione ! end if ! end do ! write(6,*)' numinside,outside=',numinside,numoutside ! write(6,*)' numgtp5=',numgtp5 ! write(6,*)' dx1min,max=',dx1min,dx1max ! Get taylor matrices and invert nstart=ione nend=min(nin,ntl) do while (nstart<=nend) nthis=nend-nstart+ione if(nf==ione) then ! get weights for interpolating function ! get taylor vector(s) lp=izero x1p(nstart:nend)=one do ip=0,iord lp=lp+ione z0(1:nthis,lp)=x1p(nstart:nend) do n=1,nthis nn=n+nstart-ione iflip=ione if(ix1sign(nn)izero) then do n=1,nthis nn=n+nstart-ione iflip=ione if(ix1sign(nn)ione) then do n=1,nthis nn=n+nstart-ione iflip=ione if(ix1sign(nn)ajj(n))then ipiv(n,j)=i ajj(n)=aa end if end do end do ! swap rows, recording changed sign of determinant do n=1,ninv if(ipiv(n,j)/=j) d(n)=-d(n) end do do k=1,m do n=1,ninv s(n)=a(n,j,k) a(n,j,k)=a(n,ipiv(n,j),k) a(n,ipiv(n,j),k)=s(n) end do end do do n=1,ninv ajj(n)=a(n,j,j) if(ajj(n)==zero)then jm=j-ione write(6,100)jm 100 format(' failure in lufact:',/,' matrix singular, rank=',i3) endif end do ajj=one/ajj do i=jp,m a(1:ninv,i,j)=ajj(1:ninv)*a(1:ninv,i,j) do k=jp,m a(1:ninv,i,k)=a(1:ninv,i,k)-a(1:ninv,i,j)*a(1:ninv,j,k) end do end do end do return end subroutine vlufm