subroutine llintp(f1, slon1, slat1, dlon1, dlat1, im1, jm1, i1, j1, &
f2, slon2, slat2, dlon2, dlat2, im2, jm2, i2, j2, izp, ierr)

! This routine linearly interpolates f1 on an evenly spaced
! lat,lon grid to obtain f2 on a different lat,lon grid.
! 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).

! f2(im2,jm2):  The interpolated function with indices defined
!               the same as for 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)
! slon2         The first longitude of f2 (deg E positive)
! slat2         The first latitude  of f2 (deg N positive)
! dlon2         The longitude increment of f2 (deg)
! dlat2         The latitude  increment of f2 (deg)

! im1,jm1       The dimensions of f1
! im2,jm2       The dimensions of f2

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

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

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

! Modifications:
!   5/6/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
  integer, intent(in) :: im1
  integer, intent(in) :: jm1
  integer, intent(in) :: i1
  integer, intent(in) :: j1
  real, intent(in) :: slon2
  real, intent(in) :: slat2
  real, intent(in) :: dlon2
  real, intent(in) :: dlat2
  integer, intent(in) :: im2
  integer, intent(in) :: jm2
  integer, intent(in) :: i2
  integer, intent(in) :: j2
  integer, intent(in) :: izp

  real, intent(out), dimension(im2, jm2) :: f2
  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 :: rlonn1, rlonx1
  real :: rlatn1, rlatx1
  integer :: i, j
  integer :: i00, j00
  real :: f00, f01, f10, f11
  real :: rlon, rlat
  real :: rlon00, rlat00
  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 .or. i2 > im2 .or. j2 > jm2) then
    ierr = 2
    return
  end if

  if (dlat1 <= 0.0 .or. dlon1 <= 0.0 .or. dlat2 <= 0.0 .or. dlon2 <= 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)

! Start loop for f2 points
  do j = 1, j2
    do i = 1, i2
      rlon = slon2 + dlon2*float(i - 1)
      rlat = slat2 + dlat2*float(j - 1)

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

      if (izp /= 1) then
!       Adjust f2 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 current f2 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
      dx0 = dlon1*adtr*cos(dtr*(rlat00))
      dx1 = dlon1*adtr*cos(dtr*(rlat00 + dlat1))
      dy = dlat1*adtr

!     Calculate the x,y coordinates of the current f2 point
      x = adtr*(rlon - rlon00)*cos(dtr*rlat)
      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 and then go to the next f2 point
      f2(i, j) = a + b*x + c*y + d*x*y
    end do
  end do

  return
end subroutine llintp