subroutine tspdcal(tlat, tlon, ttime, mft, rmiss, cx, cy, cmag)
! This routine calculates the components of the storm motion (cx,cy) in
! knots, given the lat/lon as a function of time. Centered or one-sided
! differences are used as a appropriate.
! Missing lat/lon points (equal to rmiss) are accounted for.

! This routine assumes lat is in degree N
!                      lon is in degree W positive

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

  implicit none

  real, intent(in), dimension(-2:mft) :: tlat
  real, intent(in), dimension(-2:mft) :: tlon
  real, intent(in), dimension(-2:mft) :: ttime

  integer, intent(in) :: mft
  real, intent(in) :: rmiss

  real, intent(out), dimension(-2:mft) :: cx
  real, intent(out), dimension(-2:mft) :: cy
  real, intent(out), dimension(-2:mft) :: cmag

! Local variables
  integer, parameter :: itmx = 500
  real, dimension(-3:itmx) :: ttlat
  real, dimension(-3:itmx) :: ttlon
  real, parameter :: dtnmi = 60.0
  real, parameter :: dtr = 3.14159/180.0
  integer :: k
  integer :: im, i0, ip, icode
  real :: dlat, dlon, dt, alat, cfac

! Check to make sure mft is not too big
  if (mft > itmx) then
    write (6, *) 'mft too large in routine tspdcal, increase itmx'
    stop
  end if

! Initialize variables to missing
  do k = -2, mft
    cx(k) = rmiss
    cy(k) = rmiss
    cmag(k) = rmiss
  end do

  do k = -3, itmx
    ttlat(k) = rmiss
    ttlon(k) = rmiss
  end do

  do k = -2, mft
    ttlat(k) = tlat(k)
    ttlon(k) = -tlon(k)
  end do

  do k = -2, mft
    im = 0
    i0 = 0
    ip = 0
    if (ttlat(k - 1) < rmiss) im = 1
    if (ttlat(k) < rmiss) i0 = 1
    if (ttlat(k + 1) < rmiss) ip = 1

    icode = ip + 10*i0 + 100*im

    if (icode == 101 .or. icode == 111) then
!     Use centered differences
      dlat = (ttlat(k + 1) - ttlat(k - 1))
      dlon = (ttlon(k + 1) - ttlon(k - 1))
      dt = (ttime(k + 1) - ttime(k - 1))
      alat = 0.5*(ttlat(k + 1) + ttlat(k - 1))
      cfac = cos(dtr*alat)

      cy(k) = dtnmi*dlat/dt
      cx(k) = cfac*dtnmi*dlon/dt
      cmag(k) = sqrt(cx(k)**2 + cy(k)**2)
    elseif (icode == 011) then
!     Use forward difference
      dlat = (ttlat(k + 1) - ttlat(k))
      dlon = (ttlon(k + 1) - ttlon(k))
      dt = (ttime(k + 1) - ttime(k))
      alat = 0.5*(ttlat(k + 1) + ttlat(k))
      cfac = cos(dtr*alat)

      cy(k) = dtnmi*dlat/dt
      cx(k) = cfac*dtnmi*dlon/dt
      cmag(k) = sqrt(cx(k)**2 + cy(k)**2)
    elseif (icode == 110) then
!     Use backward difference
      dlat = (ttlat(k) - ttlat(k - 1))
      dlon = (ttlon(k) - ttlon(k - 1))
      dt = (ttime(k) - ttime(k - 1))
      alat = 0.5*(ttlat(k) + ttlat(k - 1))
      cfac = cos(dtr*alat)

      cy(k) = dtnmi*dlat/dt
      cx(k) = cfac*dtnmi*dlon/dt
      cmag(k) = sqrt(cx(k)**2 + cy(k)**2)
    else
!     Two of three points are missing. Can not calculate speed components
      cx(k) = rmiss
      cy(k) = rmiss
      cmag(k) = rmiss
    end if
  end do

  return
end subroutine tspdcal

subroutine tspdcale(tlat, tlon, ttime, mft, rmiss, cx, cy, cmag)
! This routine calculates the components of the storm motion (cx,cy) in
! knots, given the lat/lon as a function of time. Centered or one-sided
! differences are used as a appropriate.
! Missing lat/lon points (equal to rmiss) are accounted for.

! This routine assumes lat is in degree N
!                      lon is in degree E positive, deg W negative

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

  implicit none

  real, intent(in), dimension(-2:mft) :: tlat
  real, intent(in), dimension(-2:mft) :: tlon
  real, intent(in), dimension(-2:mft) :: ttime

  integer, intent(in) :: mft
  real, intent(in) :: rmiss

  real, intent(out), dimension(-2:mft) :: cx
  real, intent(out), dimension(-2:mft) :: cy
  real, intent(out), dimension(-2:mft) :: cmag

! Local variables
  integer, parameter :: itmx = 500
  real, dimension(-3:itmx) :: ttlat
  real, dimension(-3:itmx) :: ttlon
  real, parameter :: dtnmi = 60.0
  real, parameter :: dtr = 3.14159/180.0
  integer :: k
  integer :: im, i0, ip, icode
  real :: dlat, dlon, dt, alat, cfac

! Check to make sure mft is not too big
  if (mft > itmx) then
    write (6, *) 'mft too large in routine tspdcal, increase itmx'
    stop
  end if

! Initialize variables to missing
  do k = -2, mft
    cx(k) = rmiss
    cy(k) = rmiss
    cmag(k) = rmiss
  end do

  do k = -3, itmx
    ttlat(k) = rmiss
    ttlon(k) = rmiss
  end do

  do k = -2, mft
    ttlat(k) = tlat(k)
    ttlon(k) = tlon(k)
  end do

  do k = -2, mft
    im = 0
    i0 = 0
    ip = 0
    if (ttlat(k - 1) < rmiss) im = 1
    if (ttlat(k) < rmiss) i0 = 1
    if (ttlat(k + 1) < rmiss) ip = 1

    icode = ip + 10*i0 + 100*im

    if (icode == 101 .or. icode == 111) then
!     Use centered differences
      dlat = (ttlat(k + 1) - ttlat(k - 1))
      dlon = (ttlon(k + 1) - ttlon(k - 1))
      dt = (ttime(k + 1) - ttime(k - 1))
      alat = 0.5*(ttlat(k + 1) + ttlat(k - 1))
      cfac = cos(dtr*alat)

      cy(k) = dtnmi*dlat/dt
      cx(k) = cfac*dtnmi*dlon/dt
      cmag(k) = sqrt(cx(k)**2 + cy(k)**2)
    elseif (icode == 011) then
!     Use forward difference
      dlat = (ttlat(k + 1) - ttlat(k))
      dlon = (ttlon(k + 1) - ttlon(k))
      dt = (ttime(k + 1) - ttime(k))
      alat = 0.5*(ttlat(k + 1) + ttlat(k))
      cfac = cos(dtr*alat)

      cy(k) = dtnmi*dlat/dt
      cx(k) = cfac*dtnmi*dlon/dt
      cmag(k) = sqrt(cx(k)**2 + cy(k)**2)
    elseif (icode == 110) then
!     Use backward difference
      dlat = (ttlat(k) - ttlat(k - 1))
      dlon = (ttlon(k) - ttlon(k - 1))
      dt = (ttime(k) - ttime(k - 1))
      alat = 0.5*(ttlat(k) + ttlat(k - 1))
      cfac = cos(dtr*alat)

      cy(k) = dtnmi*dlat/dt
      cx(k) = cfac*dtnmi*dlon/dt
      cmag(k) = sqrt(cx(k)**2 + cy(k)**2)
    else
!     Two of three points are missing. Can not calculate speed components
      cx(k) = rmiss
      cy(k) = rmiss
      cmag(k) = rmiss
    end if
  end do

  return
end subroutine tspdcale