# QSUB -eo
# QSUB -s /bin/csh
# QSUB -r read -q reg -lT 4000 -lM 2Mw

ja
echo start "$user"$$
set x = $TMPDIR
cd $x

# move data file from NCAR Mass Store to local disk file
#msread indat /JOSEPH/REANAL/DS463/SF710430
msread indat /JOEY/AFC.DEC68.NEWFLG2

cat << "EOF" > tsrc.f

      program chkhds
c
c  this program read the 1962-1972 NMC B-3 ADP Upper Air data (ds353.0)
c  with the NCAR reanalysis header on the front end of each data record
c
c
c  nbf    buffer to hold the ncar header and the data record
c  rec    buffer that holds the data record in the original format
c nbytes  length in bytes of the data buffer (rec)
c  iobuf  work buffer used by rptin
c  ista   6 byte ascii string used to hold station id
c hdrcb   this common block contains the information from the ncar reanalysis
c         header

      character nbf*96
      character rec*65
      dimension iobuf(1006)
      character fname*64,ista*6

      common /hdrcb/lthr(1),nub,nhtyp,nbhsup,nfn,ndatyp,ndasrc,nstat,
     2 ndty,nyr,nmo,ndy,nhr,nlty,nlon,nlat,nele,npty,npid,ntysup

      DATA IOBUF /1006*0/

       fname='indat'
       open(10,file=fname,form='unformatted',status='old')

      iun=10
      kmax=12
      nrec=0
      npr = 0

  10  continue
      do k =1,96
	 nbf(k:k) =' '
      enddo
      do k =1,65
	 rec(k:k) =' '
      enddo
c
c  read in one data record
c
      call rptin(iun,iobuf,nbf,nwds,jj,kmax,jeof)
      if (jeof .gt. 0 ) go to 99
      nrec=nrec+1
       if (nrec .gt. 1000) goto 99
c
c   decode the ncar reanalysis header
c
      call hdrrd(nbf,ista,xlat,xlon,xele)
c
c  put the original data record into its own buffer
c

      nbytes = (nwds-3)*8 - (nub*4)/8
c
c     rec(1:65) = nbf(25:89)
      rec(1:nbytes) = nbf(25:25+nbytes)


c
c   decode the original data record
c   neither the sea level pressure nor the surface pressure are decoded
c   the original 3 or 4 byte  integer value is returned in the buffer,or
c   -9999 if the field is missing
c
c    dd =  wind direction
c    ff =  wind speed
c  islp =  undecoded sea level pressure
c     t =  temperature
c    td =  dew point
c isfcp =  undecoded surface pressure
ciunits =  units flag that is stored in byte 14 of the data record
c the missing value for all variables is -9999.
c

      call decode_afc(rec, dd,ff,islp,t,td,isfcp,iunits)

      goto 10

  99  continue


      end

      subroutine decode_afc(line,dd,ff,islp,t,td,isfcp,iunits)
      character*65 line

      idd = ival(line(17:18),2)
      iff = ival(line(19:20),2)
      islp = ival(line(26:28),3)
      itt = ival(line(29:30),2)
      itd = ival(line(36:37),2)

      iunits = ival(line(15:15),1)
      if (iunits .eq. -9999) iunits = 9

      if (itt .gt. -9999) then
	 t = float(itt)
	 if (itt .ge. 50) t = 50. - float(itt)
      else
	 t = -9999.
      endif
      if (itd .gt. -9999) then
	 td = float(itd)
	 if (itd .ge. 50) td = 50. - float(itd)
c        if (itd .eq. 50 .and. itt .lt. 50) td = 00.
      else
	 td = -9999.
      endif

c     if (ipp .gt. -9999) then
c         slp = float(ipp)/10. + 900
c        if (ipp .lt. 600 ) slp = float(ipp)/10. + 1000
c     else
c        slp = -9999.
c     endif

      if (idd .ne. -9999 .and. idd .le.  36 ) then
	 dd = float(idd) * 10.
      else
	 dd = -9999.
      endif
      if (iff .ne. -9999 ) then
	 ff = float(iff)
      else
	 ff = -9999.
      endif
      ii=0
      isup = 0
      kk = 41
      isfcp = -9999
      sfcp = -9999.0


  200 ii = ii + 1
      jj = kk + (ii-1)*5
      if (line(jj:jj) .ne. ' ' ) then
	 isup = isup + 1
	 if (line(jj:jj) .eq. '6' ) then
c         isfcp = -9999
c         sfcp = -9999.0
	  isfcp = ival(line(jj+1:jj+4),4)
c         if (isfcp .eq. -9999 ) then
c            sfcp = -9999.
c         else if (line(jj+1:jj+1) .eq. '3' ) then
c            sfcp = float(mod(isfcp,3000))/10. + 1000
c         else
c            sfcp = float(isfcp)/10.
c         endif
	 endif
      else
	  return
      endif
      if (ii .eq. 5 ) return
      goto 200
 299  continue

      return
      end

      function ival(rec,nc)
      character rec*(*),fmt*4


      do k = 1,nc
	ic = ichar(rec(k:k))
	if (ic .lt. 48 .or. ic .gt. 57 ) then
	   ival = -9999
	   return
	endif
      enddo
      write(fmt,6) nc
   6  format('(i',i1,')')
      read(rec(1:nc),fmt) iv
      ival = iv

      return
      end

      subroutine hdrrd(nbf,ista,xlat,xlon,xele)
      save
      character npid*8,ista*6
      parameter (numl=20)
c     dimension nbf(*)
      character nbf*96
      common /hdrcb/lthr(1),nub,nhtyp,nbhsup,nfn,ndatyp,ndasrc,nstat,
     2 ndty,nyr,nmo,ndy,nhr,nlty,nlon,nlat,nele,npty,npid,ntysup
      dimension nsz(numl),nhs(8)
      data nsz/12,4,4,4,10,8,9,6,4,9,4,5,12,3,16,15,15,4,48,8/
      data nhs/4,8,13,17,19,0,0,0/
      iof=0
      call gbyte(nbf,nhtyp,16,4)
      nump=nhs(mod(nhtyp,8)+1)
      do 12 i=1,nump
      call gbyte(nbf,lthr(i),iof,nsz(i))
      iof=iof+nsz(i)
   12 continue
      nyr=nyr+1700
      xhr=float(nhr)*.01
      xlon=(nlon-18001)*.01
      xlat=(nlat- 9001)*.01
      xele=nele-1001.
      ista='      '
      call sbyte(ista,npid,0,48)
      return
      end






"EOF"


cft77 -e h tsrc.f

segldr tsrc.o -L /lib,/usr/lib,/usr/local/lib -l ncarm,ncaro,ncaru,net

a.out
# copy output back
#rcp fort.20 xxxx.yyy:/tmp/output$$

ja -s