program extrapolation
!------------------------------------------------------------------------------
!   NAME: extrap.f

!   LANGUAGE:  Fortran 77/90/95

!   PURPOSE: perform track extrapolation

!   DESCRIPTION: creates track extrapolation using the current position and
!            12 hr hold position to determine distance and direction.
!            Position information is then extrapolated from 00 HRS out
!            to 120 HRS.

!   INFORMATION:
!    National Hurricane Center / Tropical Prediction Center
!    Department of Commerce - NOAA / NWS

!   CATEGORY: crude track guidance

!   CALLING SEQUENCE: extrap.x al882007.com

!   INPUTS: CARQ or Compute File

!   OUTPUTS: extrapolated track data called "extrap.dat"

!   COMMON BLOCKS: See code

!   EXTERNAL or INTERNAL CALLS:

!   "External" ATCF data I/O routines
!   getARecord      - get Adeck record for specific Techname
!   getSingleTAU    - get single Tau value
!   newWriteAidRcd  - write out techname and tech data

!   "Internal"
!XTRAP - general extrapolation subroutine
!rhdst - rhumb line dist and direction between 2 points
!rltlg - extrapolate from set of points using dir and distance

!   SIDE EFFECTS: Unknown
!
!   RESTRICTIONS: Unexpected results across prime meridian

!   EXAMPLE or USAGE: extrap.x al882007.com
!                 outputs -> extrap.dat "to current directory"

!   OTHER NOTES: modified to produce results for all basins except as
!            listed in RESTRICTIONS section.

!   MODIFICATION HISTORY:
!    Written by:     UNKNOWN
!    Modified by:    Chris Sisko,  31-OCT-2006 (NOAA-NWS)

!   REVISION INFORMATION (RCS Keyword):
! $Id$

! Modifications:
! 6/2/21 (J. Dostalek) Started this Fortran 90 version
! 1/31/23 (J. Dostalek) This Fortran 90 version was created somewhat
! quickly to include in the February 2023 delivery to NHC. It gives the same 
! results as the FORTRAN 77 version using al072022.com, ep122021.com, and 
! cp012018.com as input, but will need some additional work to complete the 
! goals of the f77 to f90 conversion.
! ------------------------------------------------------------------------------
  implicit none

  include 'dataformats.inc'

  common/extrap/ifpxtrp(15, 3)

  real :: latcur, loncur, latm12, lonm12

  character(len=8) :: strmid
  character(len=1) ::  NS, EW
  character(len=4) ::  techname
  character(len=10) :: aymdh
  character(len=50) :: input_file

  integer :: ifpxtrp
  integer :: i, j
  integer :: ios
  integer :: istat
  integer :: iymdh

  type(AID_DATA) comRcd, tauData

! **   Zero the final array
  do i = 1, 15
    do j = 1, 3
      ifpxtrp(i, j) = 0
    end do
  end do

! **   Get the command line parameter and the storms directory path
  call getarg(1, input_file)
  strmid = input_file(1:8)

! **   Open the input file
  open (21, file=input_file, status='old', iostat=ios, err=1010)

! **   Read in the compute data and close the file
  call getARecord(21, "CARQ", comRcd, istat)
  close (21)

  if (istat == 0) then
    print *, ' ERROR - reading file = ', input_file, ' istat = ', istat
    stop
  end if

! **   Read in the current data
  call getSingleTAU(comRcd, 0, tauData, istat)
  if (istat /= 1) then
    print *, ' ERROR - reading data = ', input_file, ' istat = ', istat
    stop
  end if

  aymdh = tauData%aRecord(1)%DTG
  read (tauData%aRecord(1)%DTG, '(i10)') iymdh

  latcur = tauData%aRecord(1)%lat
  loncur = tauData%aRecord(1)%lon

  NS = trim(tauData%aRecord(1)%NS)
  EW = trim(tauData%aRecord(1)%EW)
  if (NS == 'S') latcur = -(latcur)
  if (EW == 'W') loncur = -(loncur)

! **   Read in the minus 12-hour old data
  call getSingleTAU(comRcd, -12, tauData, istat)
  if (istat /= 1) then
    print *, ' ERROR - reading data = ', input_file, ' istat = ', istat
    stop
  end if

  latm12 = tauData%aRecord(1)%lat
  lonm12 = tauData%aRecord(1)%lon

  NS = trim(tauData%aRecord(1)%NS)
  EW = trim(tauData%aRecord(1)%EW)
  if (NS == 'S') latm12 = -(latm12)
  if (EW == 'W') lonm12 = -(lonm12)

  write (6, 1) strmid
1 format(' eXTRaPolation forecast for ', a8, //)

  call XTRAP(iymdh, latcur, loncur, latm12, lonm12, STRMID)

! **   Open the output file, write out the data and close the file
  open (31, file='extrap.dat', status='unknown')

!**   Write the forecasts to a temporary adeck for merging with the
! **     actual adeck and close the file
  techname = 'XTRP'
  call newWriteAidRcd(31, strmid, aymdh, techname, ifpxtrp)

  close (31)

  stop '***** EXTRAPOLATION FORECASTS ARE FINISHED *****'

!**   Error messages
1010 print *, ' ERROR - opening file = ', input_file, ' istat = ', istat
  stop

end
!***********************************************************************
subroutine XTRAP(iymdh, latcur, loncur, latm12, lonm12, STRMID)

!**   EXTRAPOLATION FORECAST USING THE PAST 12 HOUR MOTION

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

  implicit none

  integer, intent(in) :: iymdh
  real, intent(in) :: latcur
  real, intent(in) :: loncur
  real, intent(in) :: latm12
  real, intent(in) :: lonm12
  character(len=8), intent(in) :: strmid

  integer :: ifpxtrp
  real, dimension(15, 2) :: fpxtrp
  integer, dimension(15) :: itime
  character(len=1) ::  EW, NS
  character(len=3), dimension(15) :: ftime
  integer :: i
  real :: fprite
  real :: dir
  real :: dst

  common/extrap/ifpxtrp(15, 3)

  data FTIME/'  0', ' 12', ' 24', ' 36', ' 48', ' 60', ' 72', ' 84', ' 96', &
  '108', '120', '132', '144', '156', '168'/
  data itime/0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168/

! **   WRITE THE INPUT PARAMETER TO THE PRINTER
  write (*, '('' EXTRAPOLATION FORECAST'',/)')
  write (*, '('' INPUT PARAMETERS FOR TROPICAL CYCLONE '', A,'' ON YMDH = '', &
  I10.10, / )') STRMID, iymdh

  call rhdst(latm12, lonm12, latcur, loncur, dir, dst)

  fpxtrp(1, 1) = latcur
  fpxtrp(1, 2) = loncur

  call rltlg(latcur, loncur, fpxtrp(2, 1), fpxtrp(2, 2), dir, dst*1.0)
  call rltlg(latcur, loncur, fpxtrp(3, 1), fpxtrp(3, 2), dir, dst*2.0)
  call rltlg(latcur, loncur, fpxtrp(4, 1), fpxtrp(4, 2), dir, dst*3.0)
  call rltlg(latcur, loncur, fpxtrp(5, 1), fpxtrp(5, 2), dir, dst*4.0)
  call rltlg(latcur, loncur, fpxtrp(6, 1), fpxtrp(6, 2), dir, dst*5.0)
  call rltlg(latcur, loncur, fpxtrp(7, 1), fpxtrp(7, 2), dir, dst*6.0)
  call rltlg(latcur, loncur, fpxtrp(8, 1), fpxtrp(8, 2), dir, dst*7.0)
  call rltlg(latcur, loncur, fpxtrp(9, 1), fpxtrp(9, 2), dir, dst*8.0)
  call rltlg(latcur, loncur, fpxtrp(10, 1), fpxtrp(10, 2), dir, dst*9.0)
  call rltlg(latcur, loncur, fpxtrp(11, 1), fpxtrp(11, 2), dir, dst*10.0)
  call rltlg(latcur, loncur, fpxtrp(12, 1), fpxtrp(12, 2), dir, dst*11.0)
  call rltlg(latcur, loncur, fpxtrp(13, 1), fpxtrp(13, 2), dir, dst*12.0)
  call rltlg(latcur, loncur, fpxtrp(14, 1), fpxtrp(14, 2), dir, dst*13.0)
  call rltlg(latcur, loncur, fpxtrp(15, 1), fpxtrp(15, 2), dir, dst*14.0)

  NS = 'N'
  if (latcur < 0.0) NS = 'S'
  EW = 'E'
  if (loncur < 0.0) EW = 'W'

  write (*, '( '' FORECAST POSITIONS'', / )')

  write (*, '( '' INITIAL = '', F5.1, A2 ,F7.1, A2 )') latcur, NS, loncur, EW
  write (*, '(''----------------------------'')')

  do i = 1, 15

!   **   Form the integer values for the adeck
    ifpxtrp(i, 1) = nint(fpxtrp(i, 1)*10.0)
    ifpxtrp(i, 2) = -(nint(fpxtrp(i, 2)*10.0))

!   ** high latitude check
    if (abs(fpxtrp(i, 1)) > 88.0) then

      ifpxtrp(i, 1) = 0
      ifpxtrp(i, 2) = 0
      fpxtrp(i, 1) = 0.0
      fpxtrp(i, 2) = 0.0

      write (*, '(''Error: High Latitiude error!!!'')')

    end if

!   **   Write out the forecast for the correct hemisphere
    EW = 'E'
    if (FPXTRP(i, 2) < 0.0) EW = 'W'

    NS = 'N'
    if (FPXTRP(i, 1) < 0.0) NS = 'S'

    FPRITE = FPXTRP(i, 2)

    if (FPXTRP(i, 2) > 180.0) then

      FPRITE = FPXTRP(i, 2) - 360.0
      EW = 'W'

    else if (FPXTRP(i, 2) < -180.0) then

      FPRITE = 360.0 + FPXTRP(i, 2)
      EW = 'E'

    end if

    write (*, '( 1X, A, '' HR = '',F5.1,A2,F7.1,A2,I5,2(f7.1))') FTIME(i), &
    FPXTRP(i, 1), NS, FPRITE, EW, itime(i), fpxtrp(i, 1), fpxtrp(i, 2)

  end do

  return
end subroutine XTRAP
!***********************************************************************
subroutine rhdst(flat, flng, tlat, tlng, dir, dst)

!**   THIS IS A "GHDST" SUBROUTINE; FROM POINT (FLAT,FLNG) TO POINT
!**     (TLAT,TLNG) FINDS THE DIRECTION (DIR) AND DISTANCE (DST) IN NM
!**     DIRECTION IS CALCULATED ACCORDING TO RHUMB-LINE DIRECTION.

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

  implicit none

  real, intent(in) :: flat
  real, intent(in) :: flng
  real, intent(in) :: tlat
  real, intent(in) :: tlng

  real, intent(out) :: dir
  real, intent(out) :: dst

  real :: rnum
  real :: td1, td2, rltd1, rltd2, denom, rmag
  integer :: icrs

  dir = 0.0
  dst = 0.0

  rnum = (flng - tlng)*3.1415926535898/180.0
  td1 = tan((45.0 + (0.5*tlat))*3.1415926535898/180.0)
  td2 = tan((45.0 + (0.5*flat))*3.1415926535898/180.0)
  rltd1 = alog(td1)
  rltd2 = alog(td2)
  denom = rltd1 - rltd2
  rmag = (rnum*rnum) + (denom*denom)

  if (rmag /= 0.0) dir = atan2(rnum, denom)*180.0/3.1415926535898
  if (dir <= 0.0) dir = 360.0 + dir

  dst = 60.0*abs(flat - tlat)/abs(cos(dir*3.1415926535898/180.0))
  icrs = dir + 0.5
  if (icrs == 90 .or. icrs == 270) dst = 60.0*abs(flng - tlng)* &
  cos(flat*3.1415926535898/180.0)

  dir = float(icrs)
  icrs = dst*10.0 + 0.5
  dst = float(icrs)/10.0

  return
end subroutine rhdst
!*******************************************************************
subroutine rltlg(flat, flng, tlat, tlng, dir, dst)

!**   THIS IS A "GLTLG" SUBROUTINE; FROM POINT (FLAT,FLNG) AND DIR,DST
!**       FINDS THE END POINT (TLAT,TLNG)
!**       DIRECTION IS THE RHUMB-LINE DIRECTION.

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

  implicit none

  real, intent(in) :: flat
  real, intent(in) :: flng
  real, intent(in) :: dir
  real, intent(in) :: dst

  real, intent(out) :: tlat
  real, intent(out) :: tlng

  integer :: icrs
  real :: crpd
  real :: td1, td2, rltd1, rltd2, denom
  real :: dlon

  tlat = 0.0
  tlng = 0.0

  icrs = dir
  if (icrs == 90 .or. icrs == 270) then

    dlon = dst/(60.0*cos(flat*3.1415926535898/180.0))
    if (icrs == 90) tlng = flng - dlon
    if (icrs == 270) tlng = flng + dlon
    tlat = flat

  else

    crpd = dir*3.1415926535898/180.0
    tlat = flat + (dst*cos(crpd)/60.0)
    if (tlat > 89.0) tlat = 89.0
    if (tlat < -89.0) tlat = -89.0

    td1 = tan((45.0 + (0.5*tlat))*3.1415926535898/180.0)
    td2 = tan((45.0 + (0.5*flat))*3.1415926535898/180.0)
    rltd1 = alog(td1)
    rltd2 = alog(td2)
    denom = rltd1 - rltd2
    tlng = flng - ((tan(crpd)*denom)*180.0/3.1415926535898)

  end if

  icrs = tlat*10.0 + 0.5
  tlat = float(icrs)/10.0
  icrs = tlng*10.0 + 0.5
  tlng = float(icrs)/10.0

  return
end subroutine rltlg
!------------------------------------------------------------------------------

! EXPANDED REVISION INFORMATION (RCS Keywords):

!   $Id$
!   $Author$
!   $Log$
!   $Locker$

!------------------------------------------------------------------------------