subroutine gdland_table(lon_in, lat_in, dtl)

      ! VERSION 1.0.1
      ! Last updated: 02 Mar 2020
      ! Modified Feb 2023 (KM) for NCO specifications

      ! This program calculates the distance to land (dtl in km)

      ! Input: 
      !   lon_in,lat_in - (deg W neg, and deg N pos)

      ! Output: 
      !   dtl - Distance to nearest landmass in km
  
      ! Created 10 Feb 2020 (v1.0.0) (SS) 
      ! ++ Modified from gfland.f - WSP version using lookup table for dtl (not fractional land)

      ! Updated 02 Mar 2020 (v1.0.1) (KM, AB)
      ! ++ IMPLICIT NONE and freeform format
      ! ++ adjusted ioper to include versions for SHIPS and WSP (uncomment needed version)
      ! ++ dtl returns missing (9999.) in case of error instead of 1000.
      ! ++ removed checks in subroutine for lsdiag length; moved to shipst7

      IMPLICIT NONE
      ! Calling variables
      real, intent(in) :: lon_in, lat_in
      real, intent(out) :: dtl

      ! Local variables
      integer, parameter :: mx=3601, my=1801
      integer, parameter :: iperiodic = 1

      real,dimension(mx,my) :: gdtab
      character(len=200) :: include_path, fntab
      character(len=10), parameter :: gdtab_fmt='(10(f9.1))'

      integer :: iread           !has table file been read
      integer :: ierr = 1        !initialize ierr to 1
      integer :: lutab           !file number
      integer :: i, j, ii, jj    !counters 
      ! Read in from table file header
      real :: slat, elat, dlat 
      integer :: nlat 
      real :: slon, elon, dlon 
      integer :: nlon 

      ! Temporary calculations 
      real :: x, xlon, y, ylat, xy, rlon, rlat
      real :: w00, w01, w10, w11

      ! Save table to only read file once during program execution
      save gdtab 
      save slat, elat, dlat, nlat, slon, elon, dlon, nlon, iread
      data iread /1/
      integer :: ioper 
      common /oper/ ioper

      ! If running SHIPS then use this version of ioper flag
      !if     (ioper .eq. 1) then 
        call getenv("SHIPS_COEF", include_path )      ! operational SHIPS
      !else
      !  include_path="."                              ! local
      !endif 

      ! If running WSP then use this version of ioper flag 
      !if     (ioper .eq. 1) then
      !   call getenv ( "ATCFPROBLTY", include_path )      ! ATCF
      !elseif (ioper .eq. 2) then
      !   call getenv ( "PRBLTYINC", include_path )        ! IBM
      !else
      !   include_path = "include"                         ! local
      !endif

      ! Set defaults
      dtl      = 9999.!-999.

      if (iread .eq. 1) then
         ! Read the data tables

         lutab = 27
         fntab = trim(include_path)//"/gdland_table.dat"
         open(file=trim(fntab),unit=lutab,form='formatted', &
              status='old',err=927)

         read(lutab,*,err=927,end=927) slon,elon,dlon,nlon, &
                                       slat,elat,dlat,nlat

         do j=1,nlat
            read(lutab,gdtab_fmt,err=927,end=927)  (gdtab(i,j),i=1,nlon)
         enddo

         ! Add a column at lon=360 matching lon=0
         if (iperiodic .eq. 1) then
            nlon=nlon+1
            elon=elon+dlon
            do j=1,nlat
               gdtab(nlon,j)=gdtab(1,j)
            enddo
         endif

         iread=0
      endif

      ! normalize input longitudes to 0-360 regardless of input 
      rlon = mod(lon_in+360., 360.)
      rlat = lat_in

      ! Check input values
      if (rlat .lt. slat .or. rlat .gt. elat) then
         ierr=2
         call gdland_err_handling(1,6,ierr, &
                                  slon,elon,slat,elat,rlon,rlat)
      endif

      if (rlon .lt. slon .or. rlon .gt. elon) then
         ierr=3
         call gdland_err_handling(1,6,ierr, &
                                  slon,elon,slat,elat,rlon,rlat)
      endif


      ! Calculate indices of table value closest but to the lower left
      ! of the requested lat/lon
      ii = 1 + ifix( (rlon-slon)/dlon )
      jj = 1 + ifix( (rlat-slat)/dlat )

      if (ii .eq. nlon) ii = nlon-1
      if (jj .eq. nlat) jj = nlat-1

      ! Calculate normalized x,y values relative to lower left table point
      xlon = slon + dlon*float(ii-1)
      x    = (rlon-xlon)/dlon

      ylat = slat + dlat*float(jj-1)
      y    = (rlat-ylat)/dlat

      xy = x*y

      ! Interpolate table values to find dtl and afraction
      w00 = 1.0+xy-(x+y)
      w10 = x-xy
      w01 = y-xy
      w11 = xy

      dtl = w00*gdtab(ii  ,jj  ) + &
            w10*gdtab(ii+1,jj  ) + &
            w01*gdtab(ii  ,jj+1) + &
            w11*gdtab(ii+1,jj+1)

      return

  927 continue
      call gdland_err_handling(1,6,ierr, &
                               slon,elon,slat,elat,rlon,rlat)
      return 

      end subroutine gdland_table

      subroutine gdland_err_handling(ierrtype,luerr,ierr, &
                                     slon,elon,slat,elat,rlon,rlat)
      ! This routine handles error processing for tcliper.
      ! Types of errors handled (ierrtype=):
      !    -1: error in gdland_table
      ! Input
      !    ierrtype - error type to process
      !    luerr    - unit number to write to (assumes file is already open)
      !    ierr     - error code to report

      IMPLICIT NONE
 
      !++ Passed variables
      integer, intent(in) :: ierrtype, luerr, ierr
      real, intent(in) :: slon,elon,slat,elat,rlon,rlat

      if (ierrtype .eq. 1) then
         write(luerr,*) 'Error in subroutine gdland_table, ierr=',ierr
         write(luerr,*) 'slon,elon,slat,elat,rlon,rlat: ', &
                         slon,elon,slat,elat,rlon,rlat 
      else
         write(luerr,*) 'Unrecognized error in gdland_table'
         stop
      endif

      return
      end subroutine gdland_err_handling