#
# lcase
# '#' @
# QSUB -eo -s /bin/csh
# QSUB -r rdNVAP -q econ -lT 600 -lM 6Mw

cd $TMPDIR
ja
cat << "EOF" > tsrc.f
      program rdNVAP_cray
c
c This is to read the ds722.0 data files.
c To use this program, you also need gbytes routine from our anonymous ftp
c area under `ftp/libraries/gbytes/f77.f'   --- most systems
c            `ftp/libraries/gbytes/f77hp.f' --- HP systems
c
c On VAX systems, uncomment the lines that have 'VAX'.
c On SGI systems, uncomment the lines that have 'SGI'.
c
      parameter (nbit=64)                      ! number of bits in a word,
      parameter (nnx=360,nny=180,nmax=nnx*nny) ! global, 1 degree
      parameter (mxrcl=720)                    ! record length in bytes
      parameter (nword=nbit/8)                 ! number of bytes in a word,
      parameter (mxone=mxrcl/nword)            ! max number of words of ione
      dimension ione(mxone),itmp(nnx)
      dimension grd(nnx,nny)
      character one*(mxrcl),inf*72
c
      integer*4 headi(11), inu
      real*4    offset,scale,headr(4),zindef
      character tapec*4,label*40,blank6*6,blank5*5
C
      equivalence (ione,one)
c
c if VAX/VMS, SGI, change recl into words
c
c     krecl=mxrcl/4                            ! VAX,SGI
      krecl=mxrcl
c
c set iprnt=0 if you DO NOT want to see each header
c
      iprnt=1
c
c get the input filename
c
      inu=11
c     write(*,*) 'enter input NVAP data file?'
c     read(*,'(a)') inf
      inf='nvap.dat'
      open(unit=inu,file=inf,access='direct',
     &     form='unformatted',recl=krecl)
c
c-----------------------------------------------
c
      jj=2**16
      jrow=0
      irec=0
   10 continue
      irec=irec+1
      read(inu,rec=irec,err=9998) (one(j:j),j=1,mxrcl)
      if((mod(irec-1,90).eq.0) .and. (irec.gt.1))
     &  write(*,*) '..record ',irec,' done'
      if(mod(irec,nny+1).eq.1) then        ! header
        read(one(1:144),11) tapec,(headi(j),j=1,3),blank5,
     &      (headi(j),j=4,11),(headr(j),j=1,4),offset,
     &      scale, zindef, blank6, label
        if(iprnt.ne.0) then
          write(*,*)' TAPE CODE   = ',tapec
          write(*,*)' SOURCE CODE = ',headi(1)
          write(*,*)' PARAMETER   = ',headi(2)
          write(*,*)' DATA FORMAT = ',headi(3)
          write(*,*)' START DATE  = ',(headi(i),i=4,6)
          write(*,*)' END DATE    = ',(headi(i),i=7,9)
          write(*,*)' X,Y DIMENEN = ',(headi(i),i=10,11)
          write(*,*)' X,Y RESOL   = ',(headr(i),i=1,2)
          write(*,*)' START LATLON= ',(headr(i),i=3,4)
          write(*,*)' OFFSET,SCALE= ',offset,scale
          write(*,*)' INDEF VALUE = ',zindef
          write(*,*)' CHAR LABEL  = ',label
        endif
c
 11     format(a4,2i3,i1,a5,i2,i3,2i2,i3,i2,2i4,2f6.2,
     &         f7.3,f8.3,3e11.5,a6,a40)
        do j=1,nny
          do i=1,nnx
            grd(i,j)=zindef
          enddo
        enddo
      else                                 ! data
        jrow=jrow+1
        call swap2(ione,ione,mxrcl)
        call gbytes(ione,itmp,0,16,0,nnx)
        kzindef=int(zindef)
        if(kzindef.ne.-32768) write(*,*) ' BAD zindef = ',zindef,kzindef
        do i=1,nnx
          if(itmp(i).ne.kzindef) then
            if(itmp(i).lt.32768) then
              grd(i,jrow)=float(itmp(i))*scale+offset
            else
              grd(i,jrow)=float(itmp(i)-jj)*scale+offset
            endif
          endif
        enddo
        if(jrow.eq.nny) then               ! finish one set
          rmax=-99999.
          rmin=99999.
          kmxi=-99
          kmxj=-99
          kmni=-99
          kmnj=-99
          do j=1,nny
            do i=1,nnx
              if(grd(i,j).gt.rmax) then
                rmax=grd(i,j)
                kmxi=i
                kmxj=j
              endif
              if(grd(i,j).lt.rmin) then
                rmin=grd(i,j)
                kmni=i
                kmnj=j
              endif
            enddo
          enddo
          write(*,111) (headi(i),i=4,9)
  111     format('..DATE yy ddd hh ',i2,i4,i3,' to ',i2,i4,i3,' DONE..')
          write(*,*) '  max grid: ',rmax,' at i,j: ',kmxi,kmxj
          write(*,*) '  min grid: ',rmin,' at i,j: ',kmni,kmnj
          write(*,*) '  grid 1-5 at NP : ',(grd(i,1),i=1,5)
          write(*,*) '  grid 1-5 at SP : ',(grd(i,180),i=1,5)
          write(*,*) 
          write(*,*) 
          jrow=0
        endif
      endif
      go to 10
 9998 continue
      write(*,*) '..ALL DONE, last record ',irec-1
      close(unit=inu)
      stop
      end
c
c        
      subroutine swap2(in,io,nn)
c swaps bytes in groups of 2 to compensate for byte swapping within
c    words which occurs on DEC (VAX) and PC machines.
c
c in - input array to be swapped
c io - ouput array with bytes swapped
c nn - number of bytes to be swapped
c     logical*1 in(1),io(1),ih
      character*1 in(1),io(1),ih                ! DEC alpha, VAX
      do 10 i=1,nn,2
      ih=in(i)
      io(i)=in(i+1)
      io(i+1)=ih
   10 continue
      return
      end
c
c
"EOF"
#
f90 -c tsrc.f
segldr tsrc.o -L /usr/local/lib -l ncarm,ncaro,mss
#
if ( ! -x a.out ) then
  exit
endif
#
set intape = 'Y24958'
set myfile = 'tpw88.monavg_nvap.nat'

msread -fBI $intape /DSS/$intape
tar xvf $intape $myfile
mv $myfile nvap.dat

a.out
#
ls -l
ja -cst
exit