c  INTP.F - NSIDC 
c 
c  INTERPOLATES ICE CONCENTRATION DATA ONTO ARCTIC MODEL GRID 
c 

        subroutine intp (g,xg,yg,ig,jg, 
     >                       v,xr,yr,ir,jr, 
     >  phi,theta,psi,kmtg,kmtr,m,nyr,method,radius) 

        real  g(ig,jg,m,nyr),xg(ig,jg),yg(ig,jg) 
        real  xr(ir),yr(jr) 
        real  phi, theta, psi, bad,radius 
        real  gln,glt, wt, w(ir,jr), v(ir,jr,m,nyr) 
        integer i, j, ig, jg, ir, jr, m, kmtg(ig,jg),kmtr(ir,jr) 
        integer method 

        if (method.eq.1) write(*,*)'GAUSS INTERP' 
        if (method.eq.2) write(*,*)'WITCH INTERP' 
        if (method.eq.3) write(*,*)'NHBR INTERP' 
c 
        do 10 i=1,ir 
        do 10 j=1,jr 
           if (method.eq. 1 .or. method .eq. 2) then 
              do n=1,nyr 
                 do k=1,m 
                   v(i,j,k,n) = 0. ! sum of data*wt 
                 enddo 
              enddo 
                   w(i,j) = 0.  ! sum of weights 
           else 
              do n=1,nyr 
                 do k=1,m 
                   v(i,j,k,n)=0. 
                 enddo 
              enddo 
           endif 
c 
           length=9999. ! initializing length, for nhbr interp 
c 
           if (kmtr(i,j).gt.0) then 
             call rotate(yr(j),xr(i),-psi,-theta,-phi,glt,gln) 
c 
             do 20 ii=1,ig 
             do 20 jj=1,jg 
                   if (kmtg(ii,jj).gt.0) then 
                     call dist(glt,gln,yg(ii,jj),xg(ii,jj),dst) 
c 
c the following IF statement is for determining 
c method of interpolation 
                     if(method.eq.3) then 
                        if(dst .lt. length) then 
                           lat=jj    ! finding index of nrst 
                           lng=ii    ! geog. coords 
                           length=dst 
                        endif 
                     else 
c 
                        if (method .eq. 1) 
     >                       wt = exp(-(dst/radius)**2) 
                        if(method .eq. 2) 
     >                       wt = (1./(1.+(dst/radius)**2))**2 
                        w(i,j) = wt + w(i,j) 
                        do n=1,nyr 
                           do k=1,m 
                           v(i,j,k,n) = v(i,j,k,n) +g(ii,jj,k,n)*wt 
                           enddo 
                        enddo 
                    endif ! method loop 
                  endif   ! kmtg loop 
20           continue 
             do n=1,nyr 
                do k=1,m 
                   if (method.eq.3)  then 
                       v(i,j,k,n)=g(lng,lat,k,n) 
                   else 
                       v(i,j,k,n) = v(i,j,k,n)/w(i,j) 
                   endif 
                enddo 
             enddo 
         endif            ! kmtr loop 
10      continue 
c 
        return 
        end