subroutine llintsp(f1, slon1, slat1, dlon1, dlat1, im1, jm1, i1, j1, fsp, &
rlon, rlat, izp, ierr)

! This version includes Jack Dostalek's pole correction (23 Jan 2006)

! This routine linearly interpolates f1 on an evenly spaced
! lat,lon grid to obtain fsp at the single point (rlon,rlat).
! The subroutine arguments are defined as follows:

! f1(im1,jm1):  Contains the function to be interpolated. The
!               first index represents longitude (increasing
!               eastward) and the second index represents
!               latitude (increasing northward).

! fsp:          The interpolated value of f1.

! slon1         The first longitude of f1 (deg E positive)
! slat1         The first latitude  of f1 (deg N positive)
! dlon1         The longitude increment of f1 (deg)
! dlat1         The latitude  increment of f1 (deg)
! rlon,rlat     The longitude,latitude of the point to interpolate f1.

! im1,jm1       The dimensions of f1

! i1,j1         The number of longitude,latitude points of f1
!               to use in the interpolation

! izp           Zonal Periodic flag:
!                   =1 when f1 is periodic in the zonal direction
!                      (normally used only when f1 spans 360 deg
!                      of longitude, starting at 0)
!                   =0 if not periodic

! ierr          Error flag: =0 for normal return
!                           =1 if fsp is outside of the f1 domain
!                              (non fatal)
!                           =2 if indices i1,j1 exceed
!                              dimension of f1 (fatal)
!                           =3 if dlon or dlat .le. 0.0 (fatal)

! Modifications:
!   4/26/21 (J. Dostalek) Started this Fortran 90 version
!-------------------------------------------------------------------------------

  implicit none

  real, intent(in), dimension(im1, jm1) :: f1
  real, intent(in) :: slon1
  real, intent(in) :: slat1
  real, intent(in) :: dlon1
  real, intent(in) :: dlat1
  real, intent(inout) :: rlon
  real, intent(inout) :: rlat
  integer, intent(in) :: im1
  integer, intent(in) :: jm1
  integer, intent(in) :: i1
  integer, intent(in) :: j1
  integer, intent(in) :: izp

  real, intent(out) :: fsp
  integer, intent(out) :: ierr

! Specify needed constants
  real, parameter :: pi = 3.14159265
  real, parameter :: dtr = pi/180.0
  real, parameter :: erad = 6371.0
  real, parameter :: adtr = erad*dtr
  real, parameter :: argmax = 89.5

  real :: rlonn1, rlonx1
  real :: rlatn1, rlatx1
  integer :: i00, j00
  real :: f00, f01, f10, f11
  real :: rlon00, rlat00
  real :: arg0, arg1, arg
  real :: dx0, dx1
  real :: dy
  real :: x, y
  real :: a, b, c, d

! Initialize error flag
  ierr = 0

! Check indices and dlat,dlon
  if (i1 > im1 .or. j1 > jm1) then
    ierr = 2
    return
  end if

  if (dlat1 <= 0.0 .or. dlon1 <= 0.0) then
    ierr = 3
    return
  end if

! Calculate min and max long,lat of f1 domain
  rlonn1 = slon1
  rlonx1 = slon1 + dlon1*float(i1 - 1)
  rlatn1 = slat1
  rlatx1 = slat1 + dlat1*float(j1 - 1)

! Check if fsp point is outside of f1 domain.
! If yes, move the fsp point to the nearest point in the
! f1 domain and set error flag.

  if (izp /= 1) then
!   Adjust fsp longitude for case without zonal periodicity
    if (rlon > rlonx1) then
      rlon = rlonx1
      ierr = 1
    end if

    if (rlon < rlonn1) then
      rlon = rlonn1
      ierr = 1
    end if
  else
!   Zonal periodic case
    if (rlon >= 360.0) rlon = rlon - 360.0
  end if

  if (rlat > rlatx1) then
    rlat = rlatx1
    ierr = 1
  end if

  if (rlat < rlatn1) then
    rlat = rlatn1
    ierr = 1
  end if

! Find the indices of the f1 point closest to,
! but with lon,lat less than the fsp point.
  i00 = 1 + ifix((rlon - rlonn1)/dlon1)
  j00 = 1 + ifix((rlat - rlatn1)/dlat1)
  if (i00 < 1) i00 = 1
  if (izp /= 1) then
    if (i00 > i1 - 1) i00 = i1 - 1
  end if
  if (j00 < 1) j00 = 1
  if (j00 > j1 - 1) j00 = j1 - 1

! Define the four f1 values to be used in the linear
! interpolation.
  f00 = f1(i00, j00)
  f01 = f1(i00, j00 + 1)

  if (izp == 1 .and. i00 == i1) then
    f10 = f1(1, j00)
    f11 = f1(1, j00 + 1)
  else
    f10 = f1(i00 + 1, j00)
    f11 = f1(i00 + 1, j00 + 1)
  end if

! Calculate the lon,lat of the point i00,j00
  rlon00 = rlonn1 + dlon1*float(i00 - 1)
  rlat00 = rlatn1 + dlat1*float(j00 - 1)

! Calculate the x,y distances between the four f1 points
! where x,y = 0,0 at i00,j00

  arg0 = abs(rlat00)
  arg1 = abs(rlat00 + dlat1)
  arg = abs(rlat)
  arg0 = min(arg0, argmax)
  arg1 = min(arg1, argmax)
  arg = min(arg, argmax)

  dx0 = dlon1*adtr*cos(dtr*(arg0))
  dx1 = dlon1*adtr*cos(dtr*(arg1))
  dy = dlat1*adtr

! Calculate the x,y coordinates of the current f2 point
  x = adtr*(rlon - rlon00)*cos(dtr*arg)
  y = adtr*(rlat - rlat00)

! Calculate the coefficients for the linear interpolation
  a = f00
  b = (f10 - f00)/dx0
  c = (f01 - f00)/dy
  d = (f11 - f00 - b*dx1 - c*dy)/(dx1*dy)

! Perform interpolation
  fsp = a + b*x + c*y + d*x*y

  return
end subroutine llintsp