subroutine hinterp(fieldin,fieldint,hlon,hlat,nxin,nyin,nzin,im,jm,nm)
!---------------------------------------------------------------------------c
!          horizontal interpolate to get the initial data at h points       c
!                                                                           c
!          --> fieldin    input data field                                  c
!          <--fieldint    output interpolated dataset                       c
!          --> zmax       number of vertical levels in input data           c
!          --> gds        gds of input data                                 c
!                                                                           c
!---------------------------------------------------------------------------c

        implicit none

        integer, intent(in):: im,jm,nm
        
        integer::nxin,nyin,nzin,l,i,j,n
        real,dimension(nxin,nyin,nzin):: fieldin
        real,dimension(0:im+1,0:jm+1,nm,nzin):: fieldint
        REAL,dimension(3,0:im+1,0:jm+1,nm)::coh,cov
        integer,dimension(4,0:im+1,0:jm+1,nm)::inh,jnh,inv,jnv
        real,dimension(nxin,nyin)::tmp
        real,dimension(0:im+1,0:jm+1,nm)::tmpout
        real,dimension(0:im+1,0:jm+1,nm)::hlon,hlat
        
        call gtll(coh,inh,jnh,hlon,hlat,nxin,nyin,im,jm,nm)
!        call gtll2(coh,inh,jnh,hlon,hlat,nxin,nyin,im,jm,nm)

!
! *** Horizontally interpolate input data (ht, mr, u, and v) to ETA grid 
!        using bilinear interpolation.
!
        DO L=1,nzin
          tmp=fieldin(:,:,l)
          call bilinb(coh,inh,jnh,tmp,tmpout,1,nxin,nyin,im,jm,nm) 
!          call polinb(coh,inh,jnh,tmp,tmpout,1,nxin,nyin,im,jm,nm) 
          fieldint(:,:,:,l)=tmpout
        ENDDO

      END SUBROUTINE hinterp