subroutine mapirr2reg
      
! mapping of K_qlam, K_qphi, J_mr from the irregular TIEGCM grid
! to the regular grid by cubic polynomials
! assume same longitudinal points
      
      use param_module,only: mlat,rmlat,mlon
      use eqvcur_module,only: vallat,reglat,
     |    kqlam_org,kqlam,
     |    kqphi_org,kqphi

      implicit none
      integer, parameter :: nint=4
      integer :: i,j,ns,ns_pre,jx,ishift
      real :: yint,dy,dif,dift
      
      do j = 1,rmlat	! loop over all regular latitudes
	 
! search for interval so that reglat in [vallat(-2),vallat(+2)]
         ns     = 1
         ns_pre = 1
         dif = abs(reglat(j)-vallat(ns_pre))	! start from previous x value
         do jx = ns_pre,mlat
	   dift = abs(reglat(j)-vallat(jx))
	   if(dift.lt.dif) then
	      ns = jx
	      dif = dift
	   endif
         enddo
	 
	 ns_pre = ns
  	
! check on which side reglat(ns) is lying and take corrsp. vallat interval        
	 ishift = -1
	 if(reglat(j).lt.0.and.reglat(j)-vallat(ns).gt.0) ishift=-2	! SH
	 if(reglat(j).ge.0.and.reglat(j)-vallat(ns).lt.0) ishift=-2	! NH
	 
	 ! am 4/2011 change (nint-1) to (nint) otherwise out of bound
	 ! with interpolation
	 ! hope gives reasonable results
	 ns = ns+ishift	! starting point for x,y in polint
	 if(ns.lt.1) then	   ! shift if at the beginning
	    ns = 1
	 elseif(ns+(nint).gt.mlat) then  ! shift if at the end
	    ns = mlat - (nint)
	 endif
	 
!         write(6,'(a7,i4,2f15.8)') 'xint2: ',ns,reglat(j),vallat(ns)
         
	 do i = 1,mlon	! longitude loop
! inter-/ extrapolate 
           call polint(vallat(ns:ns+4),kqlam_org(i,ns:ns+4),nint,
     |         reglat(j),kqlam(i,j),dy) 
           call polint(vallat(ns:ns+4),kqphi_org(i,ns:ns+4),nint,
     |         reglat(j),kqphi(i,j),dy) 
	 enddo
	 
      enddo
       
      return
      
      end subroutine mapirr2reg
!---------------------------------------------------------------
      subroutine polint(xa,ya,n,x,y,dy)
      
! polynomial interpolation/ extrapolation
! Numerical Recipes in Fortran77 Vol 1 pp. 103   
!   
      implicit none
      integer, intent(in) :: n	! n-1 polynomial order
      real, intent(in) :: x,	! interpolation value
     |   xa(n),			! x values
     |   ya(n)			! y values
      real, intent(out) :: y,	! interpolated y-value
     |   dy 			! error of y
      
      integer :: i,m,ns
      real :: den,dif,dift,ho,hp,w,c(n),d(n)
      
      ns = 1
      dif = abs(x-xa(1))
      
! find closest entry in table      
      do i = 1,n
        dift= abs(x-xa(i))
	if(dift.lt.dif) then
	   ns = i
	   dif = dift
	endif
	c(i) = ya(i)	! initialize the tableau of c and d
	d(i) = ya(i)
      enddo
      
      y = ya(ns)	! initial approximation to y
      ns = ns-1
      
      do m = 1,n-1	! for each column of the tableau
         do i = 1,n-m	! loop over the current c, d & update them
	   ho = xa(i) - x
	   hp = xa(i+m) - x
	   w = c(i+1)-d(i)
	   den = ho-hp
	   if(den == 0.) then	! if two input xa is identical
	     write(6,'(a)') 'failure in interpolation'
	     stop
	   endif
	   den = w/den
	   d(i) = hp*den
	   c(i) = ho*den
	 enddo
	 
	 if(2*ns.lt.n-m) then	! decide which correction,c d, to use
	    dy = c(ns+1)	! correction
	 else
	    dy = d(ns)	! correction
	    ns = ns -1
	 endif
	 y = y + dy	! update
       enddo
	 
      return
      
      end subroutine polint