To: wen@climate
From: Kent Hills <hills@mail630.gsfc.nasa.gov>
Subject: Fortran program HFE.FOR to process Apollo HFE data
Cc: 
Bcc: 

      program hfe

c H. Kent Hills, QSS, (NSSDC)
c June 7, 2005
c this reads and converts Apollo 15 and 17 HFE data sets:
c   NSSDC ID: PSPG-00093 (old ID = 71-063C-06A)
c   NSSDC ID: PSPG-00022 (old ID = 72-096C-01A)

      character*4 buff(700),dummy_char,blanks(8)
      integer*4 number(700), cuff(700), buff4(700), skip
      integer*4 sign, mantissa,ch,power, whole, fraction
      real*4 value(700)
      integer*4 stringend
      character*1 string(2800), reverse(2800)
      equivalence(buff4,buff)
      equivalence(cuff,string)
      equivalence(reverse,number)
      data blanks/8*'    '/
      nrec = 0  ! Counts the input records read
      inrec = 1200 ! Limits the number of input blocks read (1200 will get all)
c     iproc = records to process; program stops after selected record is done.

c  A15P1#3:
      infile1_1 = 290  ! The number of records in Tape 1, File 1
      infile1_2 = 291  ! Tape1, File 2
      infile1_4 = 7
      infile1_5 = 7
      infile1_3 = 299

c A15P2#3:
c     infile2_1 = 280  ! The number of records in Tape 1, File 1
c     infile2_2 = 281  ! Tape1, File 2
c     infile2_4 = 6
c     infile2_5 = 7
c     infile2_3 = 289

      nzero = 0  ! number of cases where number(m) = 0 input.

c      The above numbers of records are from the documentation, but corrected
c  to match the tapes in the cases of 280 and 290. The Tape 1 here
c  is DS005595, which contains only 2 files, apparently one for probe 1 and one
c  for probe 2 of the HFE on Apollo 15. Within this one file, the first N
c  records are for documented file 1, the next N2 records are for documented
c  file 2, etc.

c      To process a given file,  skip n1 records, then process n2 records.
c  There is a record-counting discrepancy: NSSDC's RCOUNT program finds
c  one less record than shown in the documentation, for each file.

    3 continue
      write(6,51)
   51 format(1x,'Enter the file number desired (from documentation.',/
     1 ' Program will skip to that location by counting records, and ',/
     1 ' process those records. 5595.file1 = 11;  5595.file2 = 21')
      read(5,52) kk
   52 format(I2)
      write(13,28) kk,blanks
      write(14,28) kk,blanks
   28 format(1x,I3,8A4)

C  Apollo 15 File 1 ("A15P2#3" in documentation)

      if(kk.eq.11) then
         iwds = 3
         nend = 100 * iwds
         stringend = 1200
         skip = 0
         iproc = 280
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
  873 format(1x,'iwds,nend,stringend,skip,iproc = ',5I6)
         go to 87
      endif
      if(kk.eq.12) then
         iwds = 3
         nend = 100 * iwds
         stringend = 1200
         skip = 280
         iproc = 281
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 85
      endif
      if(kk.eq.13) then
         iwds = 7
         nend = 100 * iwds
         stringend = 2800
         skip = 574
         iproc = 289
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 85
      endif
      if(kk.eq.14) then
         iwds = 3
         nend = 100 * iwds
         stringend = 1200
         skip = 561
         iproc = 6
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 85
      endif
      if(kk.eq.15) then
         iwds = 3
         nend = 100 * iwds
         stringend = 1200
         skip = 567
         iproc = 7
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 85
      endif


C  Apollo 15 File 2 ("A15P1#3" in documentation)

      if(kk.eq.21) then
         iwds = 3
         nend = 100 * iwds
         stringend = 1200
         skip = 0
         iproc = 290
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 87
      endif
      if(kk.eq.22) then
         iwds = 3
         nend = 100 * iwds
         stringend = 1200
         skip = 290
         iproc = 291
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 85
      endif
      if(kk.eq.23) then
         iwds = 7
         nend = 100 * iwds
         stringend = 2800
         skip = 595
         iproc = 299
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 85
      endif
      if(kk.eq.24) then
         iwds = 3
         nend = 100 * iwds
         stringend = 1200
         skip = 581
         iproc = 7
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 85
      endif
      if(kk.eq.25) then
         iwds = 3
         nend = 100 * iwds
         stringend = 1200
         skip = 588
         iproc = 7
         nrec = 0
        write(6,873) iwds,nend,stringend,skip,iproc
        write(14,873) iwds,nend,stringend,skip,iproc
         go to 85
      endif

      go to 3  ! If didn't find a match

   85 continue  ! Skip

      do 86, ii=1,skip
      read(11,899,end=30) dummy_char
      iskip  = skip
  899 format(A)
   86 continue
      write(6,*) 'Finished skipping.'
      write(14,*) 'Finished skipping.'

   87 continue  ! Read and Process

    1 read(11,901,end=330) (buff(j),j=1,nend)
  901 format(700A4)

      nrec = nrec + 1
      if(nrec.gt.inrec) then
         stop 'Record Count'
         write(6,1001) nrec
         write(14,1001) nrec
 1001    format(1x,I5,' records processed. Stopped on record count.')
      endif
      if(nrec.gt.iproc) stop 'Counted End of Selected File'

  903 format(1x,3(Z8, 2x),I5)


C Reverse the 32-bit word order, to get chronologic order in this block.
      do 102 k=1,nend
      cuff(k) = buff4(nend+1-k)
  102 continue

C Within each word, reverse the byte order
      do 129 k=1,stringend,4
      reverse(k) = string((k-1)+4)
      reverse(k+1) = string((k-1) +3)
      reverse(k+2) = string((k-1) +2)
      reverse(k+3) = string((k-1) +1)
  129 continue



C***
c  call mvbits(isource,ifirst,len,idest,ifirst)
c  moves bit string of length len from isource, starting with bit ifirst,
c  and putting it into idest, starting with ifirst.
c***

      do 120 m=1,nend

      if(number(m).eq.0) then
         nzero = nzero + 1
         write(6,875) skip, nrec,m,number(m)
         write(14,875) skip, nrec,m,number(m)
  875    format(1x,'ZERO: skip,nrec,m,number = ',3I5,Z10)
         value(m) = -9999.
         go to 120
      endif

      call mvbits(number(m),8,23,mantissa,0)
      call mvbits(number(m),0,8,ch,0)
      rexp = float(ch - 128 -23)
      rnum = 2.**rexp
      value(m) = float(mantissa) * rnum
      if(bjtest(number(m),31)) value(m) = -value(m)

  120 continue

      do 121 n=1,nend,iwds
      write(13,122) (value(m),m=n,n+iwds-1)
  121 continue
  122 format(7(1x,1PE14.7))

      if(nrec.ge.iproc) go to 130   ! End of requested number processed.

      go to 1

   30 continue   ! EOF while skipping records (blocks)
      write(6,1301) iskip
      write(14,1301) iskip
 1301 format(1x,'EOF after',I5,' blocks, while skipping blocks.')
      stop

  130 continue  ! Reached requested number of processed physical records
      write(6,1303) iproc,nzero
      write(14,1303) iproc,nzero
 1303 format(1x,'Processed specified limit of ',I5,' physical records.',/
     1 1x, I6,' data values were zero, changed to -9999. on output.')
      stop

  330 continue  ! EOF while reading records (blocks)
      write(6,1302) nrec
      write(14,1302) nrec
 1302 format(1x,'EOF after',I5,' blocks, while reading blocks.'
     1 1x, I6,' logical records had zero for a value.')


      stop
      end