Fortran 95 source code to convert JPL source cord from ascii to binary

asc2eph.f95

!***********************************************************************
! ASC2EPH creates a binary format JPL Planetary Ephemeris file from
! one or more ascii text files.
!
! date          name            version
! 2018.08.22    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!***********************************************************************
!
program main
  implicit none
  ! Unix プラットフォームでなければ nrecl = 1
  integer, parameter :: nrecl  = 4
  integer, parameter :: oldmax = 400
  integer, parameter :: nmax   = 1000
  character(len= 6) :: cnam(nmax)
  character(len= 6) :: ttl(14, 3)
  character(len=12) :: header
  real(8) :: au, cval(nmax), db(3000), db2z
  real(8) :: emrat, ss(3), t1, t2
  integer :: i, j, irecsz, in, out
  integer :: k, kk, kp2, ksize
  integer :: ipt(3, 12)  ! Pointers to number of coefficients for bodies
  integer :: lpt(3)      ! Pointer to number of coefficients for lunar librations
  integer :: rpt(3)      ! Pointer to number of coefficients for lunar euler angel rates
  integer :: tpt(3)      ! Pointer to number of coefficients for TT-TDB
  integer :: n, ncon, nrout, ncoeff, nrw, numde
  logical :: first = .true.


  ! By default, the output ephemeris will span the same interval as the
  ! input ascii data file(s).  The user may reset these to other JED's.
  !
  db2z = -99999999_8
  t1   = -99999999_8
  t2   =  99999999_8
  if (nrecl /= 1 .and. nrecl /= 4) then
    write (*,*) '*** ERROR: User did not set NRECL ***'
    stop
  end if

  ! Write a fingerprint to the screen.
  !
  write (*,*) ' JPL ASCII-TO-DIRECT-I/O program. ' // &
            & ' Last modified 15-Aug-2013.'

  ! Read the size and number of main ephemeris records.
  ! (from header.xxx)
  !
  read (*,'(6X,I6)') ksize
  write (*,'(/(A,I6))') 'KSIZE = ', ksize
  irecsz = nrecl * ksize

  ! Now for the alphameric heading records (GROUP 1010)
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1010') then
    call errprt(1010, 'NOT HEADER')
  end if
  read (*,'(14A6)') ttl
  write (*,'(/(14A6))') ttl

  ! Read start, end and record span  (GROUP 1030)
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1030') then
    call errprt(1030, 'NOT HEADER')
  end if
  read (*,'(3D12.0)') ss

  ! Read number of constants and names of constants (GROUP 1040/4).
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1040 ') then
    call errprt(1040, 'NOT HEADER')
  end if
  read (*,'(I6)') n
  read (*,'(10A8)') (cnam(i), i=1,n)
  ncon = n

  ! Read number of values and values (GROUP 1041/4)
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1041') then
    call errprt(1041, 'NOT HEADER')
  end if
  read (*,'(I6)') n
  read (*,*) (cval(i), i=1,n)
  do i = 1, n
    if (cnam(i) == 'AU    ') au    = cval(i)
    if (cnam(i) == 'EMRAT ') emrat = cval(i)
    if (cnam(i) == 'DENUM ') numde = cval(i)
  end do
  write (*,'(500(/2(A8,D24.16)))') (cnam(i), cval(i), i=1,n)

  ! Zero out pointer arrays
  !
  do i = 1, 3
    do j = 1, 12
      ipt(i, j) = 0
    end do
    lpt(i) = 0
    rpt(i) = 0
    tpt(i) = 0
  end do

  ! Read pointers needed by INTERP (GROUP 1050)
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1050') then
    call errprt(1050, 'NOT HEADER')
  end if
  write (*,'(/)')
  do i = 1, 3
    read  (*,'(15I6)') (ipt(i, j), j=1,12), lpt(i), rpt(i), tpt(i)
    write (*,'(15I5)') (ipt(i, j), j=1,12), lpt(i), rpt(i), tpt(i)
  end do

  ! Open direct-access output file ('JPLEPH')

  open (unit   = 12,             &
      & file   = 'JPLEPH',       &
      & access = 'DIRECT',       &
      & form   = 'UNFORMATTED',  &
      & recl   = irecsz,         &
      & status = 'NEW')

  ! Read and write the ephemeris data records (GROUP 1070).
  !
  call nxtgrp(header)
  if (header /= 'GROUP   1070') then
    call errprt(1070, 'NOT HEADER')
  end if
  nrout = 0
  in    = 0
  out   = 0
  do
    read (*,'(2I6)') nrw, ncoeff
    if (nrw /= 0) exit
  end do
  do k = 1, ncoeff, 3
    kp2 = min(k + 2, ncoeff)
    read (*,*, iostat = in) (db(kk), kk=k,kp2)
    if (IN /= 0) stop ' Error reading 1st set of coeffs'
  end do

  do
    if ((in /= 0) .or. (db(2) >= t2)) exit
    if (2 * ncoeff /= ksize) then
      call errprt(ncoeff, ' 2*NCOEFF not equal to KSIZE')
    end if

    ! Skip this data block if the end of the interval is less
    ! than the specified start time or if the it does not begin
    ! where the previous block ended.
    !
    if  ((db(2) >= t1) .and. (db(1) >= db2z)) then
      if (first) then
        ! Don't worry about the intervals overlapping
        ! or abutting if this is the first applicable
        ! interval.
        !
        db2z  = db(1)
        first = .false.
      end if
      if (db(1) /= db2z) then
        ! Beginning of current interval is past the end
        ! of the previous one.
        !
        call errprt(nrw, 'Records do not overlap or abut')
      end if
      db2z  = db(2)
      nrout = nrout + 1
      write (12, rec=nrout+2, iostat=out) (db(k), k=1,ncoeff)
      if (out /= 0) then
        call errprt(nrout, 'th record not written because of error')
      end if

      ! Save this block's starting date, its interval span, and its end 
      ! date.
      !
      if (nrout == 1) then
        ss(1) = db(1)
        ss(3) = db(2) - db(1)
      end if
      ss(2) = db(2)

      ! Update the user as to our progress every 100th block.
      !
      if (mod(nrout, 100) == 1) then
        if (db(1) >= t1) then
          write (*,'(I6,A,F12.2)') &
            & nrout, ' EPHEMERIS RECORDS WRITTEN.  LAST JED = ', db(2)
        else
          write (*,*) ' Searching for first requested record...'
        end if
      end if
    end if
    read (*,'(2I6)', iostat=in) nrw, ncoeff
    if (in == 0)then
      do K = 1, ncoeff, 3
        kp2 = min(k + 2, ncoeff)
        read (*,*, iostat=in) (db(kk), kk=k,kp2)
        if (in /= 0) stop ' Error reading nth set of coeffs'
      end do
    end if
  end do
  write (*,'(I6,A,F12.2)') &
    & nrout, ' EPHEMERIS RECORDS WRITTEN.  LAST JED = ', db(2)

  ! Write header records onto output file.
  nrout = 1
  if (ncon <= oldmax) then
    write (12, rec=1, iostat=out) &
      & ttl, (cnam(i), i=1,oldmax), &
      & ss, ncon, au, emrat, ipt, numde, lpt, rpt, tpt
  else
    k = oldmax + 1
    write (12, rec=1, iostat=out) &
      & ttl, (cnam(i), i=1,oldmax), &
      & ss, ncon, au, emrat, ipt, numde, lpt, (cnam(j), j=k,ncon), rpt, tpt
  end if
  if (out /= 0 ) then
    call errprt(nrout, 'st record not written because of error')
  end if
  nrout = 2
  if (ncon <= oldmax) then
    write (12, rec=2, iostat=out) (cval(i), i=1,oldmax)
  else
    write (12, rec=2, iostat=out) (cval(i), i=1,ncon)
  end if
  if (out /= 0) then
    call errprt(nrout, 'nd record not written because of error')
  end if

  ! We're through.  Wrap it up.
  close (12)
  stop ' OK'

contains

  subroutine errprt(i, msg)
    implicit none
    integer, intent(IN) :: i
    character(len=*), intent(IN) :: msg

    write (*,'(A,I8,2X,A50)') 'ERROR #', i, msg

    stop ' ERROR '
  end subroutine errprt

  subroutine nxtgrp(header)
    implicit none
    character(len=*), intent(OUT) :: header
    character(len=12) :: blank

    ! Start with nothing.
    !
    header = ' '

    ! The next non-blank line we encounter is a header record.
    ! The group header and data are seperated by a blank line.
    !
    do
      if (header /= ' ') exit
      read (*,'(A)') header
    end do

    ! Found the header.  Read the blank line so we can get at the data.
    !
    if (header /= 'GROUP   1070') read (*,'(A)') blank

    return
  end subroutine nxtgrp
end program main