subroutine gtll(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,nxinp
      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

       if(nxin.eq.145)then
         nxinp=144
       else
         nxinp=nxin
       endif

       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*nxinp/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)
           inh(2,i,j,n)=inh(1,i,j,n)+1
           inh(3,i,j,n)=inh(1,i,j,n)
           inh(4,i,j,n)=inh(1,i,j,n)+1

           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)
           jnh(2,i,j,n)=jnh(1,i,j,n)
           jnh(3,i,j,n)=jnh(1,i,j,n)+1
           jnh(4,i,j,n)=jnh(1,i,j,n)+1

           if (jnh(1,i,j,n).eq.nyin) then
              print *,i,j,n, y, int(y), aphd
              jnh(3,i,j,n)=jnh(1,i,j,n)
              jnh(4,i,j,n)=jnh(1,i,j,n)
           end if

       enddo
       enddo
       enddo
      end