When we moved from ark to fits files and extensions, the 3d capability was
lost.  This file is to help you restore it, by preserving the old routines


        integer function inpark_3d(naxis, data, axdata)
                                                                       
!       Outputs.
!       The lengths of each axis.
        integer, dimension(3) :: naxis
!       The data.
        real, dimension(:,:,:) :: data
!       The axis scales.
        real, dimension(:,:) :: axdata
                                
!       Locals           
        logical :: there
        character(len=8) :: arktyp
        integer :: i, ifail, iunit, iorder
        real, external :: ripoly
        character(len=80) :: input
        double precision, dimension(0:9) :: coefs
        real :: junk1, junk2

!       The header information.
        character*80 dross(180)
        common /pasdrs/ dross

        logical :: debug
        common /bug/ debug

        if (debug) then
          print*, '@ Entering inpark_3d, with;'
          print*, '@    naxis being ', naxis(1), naxis(2), naxis(3)
          print*, '@    the shape of data being ', &
          size(data,1), size(data,2), size(data,3)
          print*, '@    and axdata being ', &
          size(axdata,1), size(axdata,2)
        end if
                   
100     call inpprm('INPARK Version 2.0.', file_in, inter_in, filex, lstfil)
        if (lstfil.eq.'end' .or. lstfil.eq.'END') then 
          inpark_3d = -1
          goto 900
        end if

!       Inquire to see if the file is there, and at the same time find
!       out some other useful facts.
        call arkinq(lstfil, filex, there, arktyp)
        if (.not. there) then
          if (inter_in) then
            print*, 'File not found, please try again.'
            goto 100
          else                            
            inpark_3d=-2
            goto 900
          end if       
        end if
        print*, 'Reading file '//lstfil(1:50)
        call clear_header()
              
!       Now go to the appropriate read routine.
        if (arktyp .eq. 'ARK_FITS') then
          call arkfit_3d(data, naxis, axdata, ifail, inter_in, lstfil)
          if (ifail .ge. 0) then
            itypef=3
            inpark_3d=3
          else                    
            inpark_3d=ifail
          end if
        else if (ArkTyp .eq. 'ARK') Then  ! ARK file
          call arkopn(iunit, lstfil, ' ', 'old', 'readonly', &
          'unformatted', 'sequential', 0)
          read(iunit) naxis(1), (dross(i), i = 1, 36)
          if (naxis(1) .gt. size(data,1)) then
            if (inter_in) then
              print*, 'Array DATA too small in inpark.  The first'
              print*, 'dimension of the array is ', size(data,1)
              print*, 'but the data has ', naxis(1), ' points.'
            end if
            inpark_3d = -99
            goto 900
          end if
          if (naxis(1) .gt. size(axdata,1)) then
            if (inter_in) then      
              print*, 'Array AXDATA too small in inpark.  The first'
              print*, 'dimension of the array is ', size(axdata,1)
              print*, 'but the data has ', naxis(1), ' points.'
            end if
            inpark_3d = -99
            goto 900
          end if
          read(iunit) (data(i,1,1), i = 1, naxis(1))
          read(iunit) coefs
          close(iunit)
          do i = 0, 9
            if (coefs(i) .ne. 0.0) iorder = i
          end do      
          do i = 1, naxis(1)
            axdata(i,1) = RIPOLY(i, coefs, iorder)
          end do      
          inpark_3d = 1
          itypef=1
        else
          ! IUE File     
          call arkopn(iunit, lstfil, ' ', 'old', 'readonly', &
          'formatted', ' ', 0)
          ! Read the header.
          read(iunit,'(a80)') dross(1)
          read(iunit,'(a80)') dross(2)
          read(iunit,'(a80)') dross(3)
          dross(4)='END'
          i=0
          do                                                             
220         Read(iunit, '(A80)', End = 230 ) Input
            read(input, *, end = 220, err = 220 ) junk1, junk2
            i = i + 1
            if (i .gt. size(data,1)) then
              if (inter_in) then
                print*, 'Array DATA too small in inpark.  The first'
                print*, 'dimension of the array is ', size(data,1)
                print*, 'but the input file is bigger.'
              end if     
              inpark_3d = -99
              goto 900
            end if
            if (i .gt. size(axdata,1)) then                
              if (inter_in) then
                print*, 'Array AXDATA too small in inpark.  The first'
                print*, 'dimension of the array is ', size(axdata,1)
                print*, 'but the input file is bigger.'
              end if
              inpark_3d = -99
              goto 900
            end if
            Read(input, *, end = 220, err = 220 ) axdata(i,1), data(i,1,1)
          End Do              
230       close(iunit)
          naxis(1) = i
          naxis(2)=1
          naxis(3)=1
         inpark_3d = 2
          itypef=2
        end if
                    
900       inter_in=.true.
        end function inpark_3d



        subroutine arkfit_3d(data, naxis, axdata, ifail, inter, lstfil)
                                                 
!       ifail is returned as 3 for success.
                                           
!       Passed variables.
        integer, dimension(3) :: naxis
        real, dimension(:,:) :: axdata
        real, dimension(:,:,:) :: data
        integer ifail
        logical inter
        character(len=*) :: lstfil
                     
!       Locals.
        integer*2 inblk(1440)
        real, dimension(2880/4) :: relblk
        double precision bscale, bzero, crval(3), cdelt(3)
        logical scale, reals
        integer i, j, k, iblk, icount, max1, max2, max3, nbit, iunit
                                 
!       Functions called.
        logical, external :: arkgsm
                                                                
        common /bug/ debug
        logical debug

        iunit=0
        call fithed(lstfil, '   ', inter, iunit, ifail)
        if (ifail .ne. 0) goto 900
        ifail=3
        inquire(iunit, nextrec=iblk)
             
        i=get_header_i('BITPIX', nbit)
        if (nbit .eq. -32) then
          reals=.true.
        else if (nbit .eq. 16) then
          reals=.false.
        else
          if (inter) then
            print*, '* Warning BITPIX =', nbit
          else
            ifail=-30
            goto 900
          end if
        end if
        if (get_header_i('NAXIS1', naxis(1)).lt.1 .or. &
          get_header_i('NAXIS2', naxis(2)).lt.1) then
          if (inter) then
            print*,'* Failed to read image size from FITS file header.'
            print*,'* Currently NAXIS1 =', naxis(1), 'NAXIS2 =', &
            naxis(2)
410         print*,'> Please give NAXIS1, NAXIS2 for this file.'
            read(*,*,err=410) naxis(1), naxis(2)
          else
            ifail=-3
            goto 900
          end if
        end if
        if (get_header_i('NAXIS3', naxis(3)) .lt. 1) naxis(3)=1
        if (debug) print*, '@ Array is ', naxis(1), 'x', naxis(2),&
        'x', naxis(3)
        scale=.false.
        if (get_header_d('BSCALE', bscale) .gt. 0) then
          if (debug) print*, '@ bscale is ', bscale
          scale=.true.
        else
          bscale=1.0d0
        end if
        if (get_header_d('BZERO', bzero) .gt. 0) then
          if (debug) print*, '@ bzero is ', bzero
          scale=.true.
        else
          bzero=1.0d0
        end if
              
!       Check that reasonable naxis(1) and naxis(2) values have been found.
        if (naxis(1).gt.size(data,1) .or. naxis(2).gt.size(data,2)) then
          if (inter) then
            print*,'* Possible error in data array dimensions.'
            print*,'* Program data array size is ', size(data,1), &
            'x', size(data,2), 'x', size(data,3)
            print*,'* Data file header gives size as ', &
            naxis(1), 'x', naxis(2), 'x', naxis(3)
          else
            ifail=-99
            goto 900
          end if
        end if
        if (size(axdata,1) .lt. max(naxis(1), naxis(2), naxis(3))) then
          if (inter) then
            print*, '* Possible error in axdata array dimensions.'
            print*, '* Program size is ', size(axdata,1)
            print*, '* Data file header gives longest axis as ', &
            max(naxis(1), naxis(2), naxis(3))
          else
            ifail=-97
            goto 900                        
          end if
        end if  

!       Create the axis information.
        if (get_header_d('CRVAL1', crval(1)) .lt. 1) crval(1)=1.0d+00
        if (get_header_d('CDELT1', cdelt(1)) .lt. 1) cdelt(1)=1.0d+00
        if (get_header_d('CRVAL2', crval(2)) .lt. 1) crval(2)=1.0d+00
        if (get_header_d('CDELT2', cdelt(2)) .lt. 1) cdelt(2)=1.0d+00
        if (get_header_d('CRVAL3', crval(3)) .lt. 1) crval(3)=1.0d+00
        if (get_header_d('CDELT3', cdelt(3)) .lt. 1) cdelt(3)=1.0d+00
        do j=1, 3
          do i=1, naxis(j)
            axdata(i,j)=real(crval(j) + dble(i-1)*cdelt(j))
          end do
        end do
                
!       Read the data.
        max1=naxis(1)
        max2=naxis(2)
        max3=naxis(3)
        if (reals) then                 
          icount=2880/4
          do k=1, max3
            do j=1, max2
              do i=1, max1
                if (icount .eq. 2880/4) then
                  if (debug) print*, '@ About to read block ', iblk
                  read(unit=iunit, rec=iblk) relblk
                  iblk=iblk+1
                  icount=0
                end if
                icount=icount+1
                data(i,j,k)=relblk(icount)
              end do
            end do
          end do
        else
          icount=2880/2
          do k=1, max3
            do j=1, max2
              do i=1, max1
                if (icount .eq. 2880/2) then
                  if (debug) print*, '@ About to read block ', iblk
                  read(unit=iunit, rec=iblk) inblk
                  iblk=iblk+1
                  icount=0
                end if
                icount=icount+1
                data(i,j,k)=inblk(icount)
              end do                         
            end do
          end do
        end if                

        if (scale) then
          if (debug) print*, '@ Applying scale factors...'
          do k=1, naxis(3)
            do j=1, naxis(2)
              do i=1, naxis(1)
                data(i,j,k)=sngl(dble(data(i,j,k))*bscale + bzero)
              end do
            end do
          end do
        end if
              
        if (debug) then
          print*, '@ Test data point (1,1,1) = ', data(1,1,1)
        end if
              
900     close(iunit)
        end subroutine arkfit_3d

