subroutine gtll2(coh,inh,jnh,hlon,hlat,nxin,nyin,im,jm,nm)
!==========================================================================c
!                                                                          c
!    calculate the position of each grid point within the projection       c
!    system                                                                c
!                                                                          c
!    <-- coh         relative positions of the h grid point                c
!    <-- inh,jnh     nearest four points indices to be used in             c
!                    the interpolation of h point                          c
!                                                                          c
!--------------------------------------------------------------------------c

      implicit none

      integer::nxin,nyin,im,jm,nm,i,j,n
      real,parameter::pi=3.141592654
      real,dimension(3,0:im+1,0:jm+1,nm)::coh
      integer,dimension(4,0:im+1,0:jm+1,nm)::inh,jnh
      real,dimension(0:im+1,0:jm+1,nm):: hlat,hlon
      real::aphd,almd,x,y


       nxin=nxin/2
       nyin=(nyin-1)/2+1

       do n=1,nm
       do i=1,im
       do j=1,jm
           aphd=hlat(i,j,n)
           almd=hlon(i,j,n)
	   if(almd.ge.2.*pi)then
	      almd=almd-2.*pi
           else if(almd.lt.0)then
	      almd=almd+2.*pi
           endif
	   x=almd*0.5*nxin/pi+1.
	   y=aphd*(nyin-1)/pi+(nyin+1)*0.5

           coh(2,i,j,n)=x-float(int(x))
           coh(3,i,j,n)=y-float(int(y))
           coh(1,i,j,n)=coh(2,i,j,n)*coh(3,i,j,n)

           inh(1,i,j,n)=int(x)*2-1
           inh(2,i,j,n)=inh(1,i,j,n)+2
           inh(3,i,j,n)=inh(1,i,j,n)
           inh(4,i,j,n)=inh(1,i,j,n)+2

	   if(inh(1,i,j,n).eq.nxin)then
	     inh(2,i,j,n)=1
	     inh(4,i,j,n)=1
           endif

           jnh(1,i,j,n)=int(y)*2-1
           jnh(2,i,j,n)=jnh(1,i,j,n)
           jnh(3,i,j,n)=jnh(1,i,j,n)+2
           jnh(4,i,j,n)=jnh(1,i,j,n)+2
	   if(i.eq.10.and.j.eq.11.and.n.eq.5)then
!	      print *,i,j,n,inh(1,i,j,n),jnh(1,i,j,n)
           endif
       enddo
       enddo
       enddo

       nxin=nxin*2
       nyin=(nyin-1)*2+1

    
      end