if [ $# -lt 3 ]
then
  echo " "
  echo "            *****     EXT2GAD     *****      "
  echo " "
  echo " The program transforms an EXTRA file in a GRADS file"
  echo " "
  echo "   Syntax: ext2gad <file1> <file2> <x_dimension>"
  echo " "
  exit 0
else
 echo " "
 echo "            *****     EXT2GAD     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "Columns number in the input file="$3
 cp $1 fort.1
fi
cat > ext2gad.f << END
      PROGRAM EXT2GADX
      parameter(mxdim=22000)
      real data(mxdim)
      character*50 infile,ofile,gadfile
      character*10 type
      data idim,jdim/ 1,1/
      data infile,ofile,gadfile,type
     $     /'                                                  ',
     $     '                                                  ',
     $     '                                                  ',
     $     
     $     '          '/
      data amiss/ 0.9e+10/
      data minlev,maxlev/999999,-999999/

c --Read command inputs:
      print*,'Input file=?'
      read(5,*) infile
      print*,'Output file basename=?'
      read(5,*) ofile
      print*,'Data type=? or idim=?'
      read(5,*) type
      lin=index(infile,' ')-1
      lout=index(ofile,' ')-1
      ltype=index(type,' ')-1
      gadfile(1:lout+4)=ofile(1:lout)//'.gad'
      print*,'Input file:           ',infile(1:lin)
      print*,'Output file basename: ',ofile(1:lout)
      print*,'Datatype:             ',type(1:ltype)
      print '(30a1)', ('-',i=1,30)
      print*
      OPEN(1,file=infile(1:lin),form='unformatted')
      irec=0
      
      READ(1,END=99) IDAT, IVAR, ILEV, LEN
      print*, 'First field: idat=',IDAT,'  ivar=',
     $        IVAR,'  ilev=',ILEV, '  ilen=',LEN
      idat1=idat
      if (type(1:4).eq.'ukmo') then
         idim=96
      elseif (type(1:6).eq.'nphope') then
         idim=61
      elseif (type(1:4).eq.'hope') then
         len=2*len
      else
         print*,'Reading x-dim from 3rd input variable (=type)'
         read(type,'(i7)') idim
         print*,'x-dim= ',idim
      endif
      if (irec.eq.1) then
         idat2=idat
      endif
      rewind 1
      jdim=len/idim
      nrecl=4*len
      OPEN(2,file=gadfile(1:lout+4),form='unformatted',
     $     access='direct',recl=nrecl)
      if (type(1:4).ne.'hope')
     $      print*,'idim=',idim,'  jdim=',jdim,' len=',len

c -------Start while-not-end-of-file unit 9:
 10   continue
         read(1,end=99) idat,ivar,ilev,len
         read(1) (data(i),i=1,len)
         minlev=min(ilev,minlev)
         maxlev=max(ilev,maxlev)
         irec=irec+1
         if (type(1:4).eq.'hope')
     $        call hp2gad(data,ivar,len,idim,jdim,imin,jmax)

         if (irec.eq.1)
     $        print*, '     Field No: ', irec,
     $        '   idat=',idat,'  ivar=', ivar,'  ilev=',ilev,
     $        '  ilen=',len,'  idim=',idim,'  jdim=',jdim,
     $        '  imin=',imin, '  jmax=',jmax

         write(2,rec=irec) 
     $        (((min(data(I+jm*idim),amiss)),I=1,idim),jm=jdim-1,0,-1)
      goto 10
 99   continue
c ------End while loop

      print*, 'Last field No: ', irec,'   idat=',idat,
     $        '  ivar=', ivar,'  ilev=',ilev, '  ilen=',len,
     $        '  idim=',idim,'  jdim=',jdim,'  imin=',imin,
     $        '  jmax=',jmax
      call wrtdes(ofile,lout,idat1,idat2,ivar,ilev,len,type,
     $     idim,jdim,imin,jmax,irec,maxlev-minlev+1)

      print*, '----- normal end of ext2gadx -----'
      stop
      end

c =====================================================================
      subroutine wrtdes(ofile,lout,idat1,idat2,ivar,ilev,len,type,
     $     idim,jdim,imin,jmax,irec,nlevz)
c*** writing description file on unit 11.

      character*50 ofile,desfile
      character*10 type
      character*3 cmonth(12)

      data (cmonth(i),i=1,12) / 'JAN','FEB','MAR','APR','MAY','JUN',
     $     'JUL','AUG','SEP','OCT','NOV','DEC'/
      data imnth,iyr1,ndmnth / 1, 1000, 1/
 
      desfile(1:lout+4)=ofile(1:lout)//'.des'
      open(3,file=desfile(1:lout+4))
      write(3,'(3a)') 'dset  ^', ofile(1:lout),'.gad'
      write(3,'(a,i8,3(a,i6),a)') 'title Conv. Extra data (idat=',
     $    idat1,' ivar= ',ivar,'  ilev= ', ilev,' len= ', len,')'
      write(3,'(a)') 'undef 0.9E+10'
C      write(3,'(a)') 'options cray_32bit_ieee big_endian'
      if (type(1:4).eq.'ukmo') then
         dlon=3.75
         dlat=2.5
         if(jdim.eq.72)then
            xlons=1.875
         elseif (jdim.eq.73) then
            xlons=0.
         endif
         xlats=-90.
      elseif(type(1:4).eq.'hope') then
         dlon=2.5
         dlat=2.5
         xlons=1.25+(imin-1)*5.
         xlats=98.75-(jmax-1)*5.-2.5
      else
         
         dlon=360./idim
         dlat=180./jdim
         xlons=0.
         xlats=-90.
      endif
      write(3,'(a,i5,a,2(f8.3,1x))')
     $     'xdef ',idim,' linear ',xlons,dlon
      write(3,'(a,i5,a,2(f8.3,1x))')
     $     'ydef ',jdim,' linear ',xlats,dlat
      write(3,'(a,i3,a)') 'zdef ',nlevz,' linear  0 1'
      if(idat1 .gt. 100000) then
         iyr1=idat1/10000
         iyr2=idat2/10000
         imnth=(idat1-iyr1*10000)/100
         imnth2=(idat2-iyr2*10000)/100
         if(imnth.gt.12)imnth=1
         if(imnth2.gt.12)imnth2=1
         iday=idat1-iyr1*10000-imnth*100
         ndmnth=(iyr2-iyr1)*12+imnth2-imnth
      endif
      ivar=mod(ivar,1000)
      if (iyr1.eq.0) iyr1=1000
      write(3,'(a,i4,2a,i4.4,1x,i2,a)') 'tdef ', irec/nlevz,
     $    ' linear ',cmonth(imnth),iyr1,1,'mo'
C      write(3,'(a,i4,2a,i4.4,1x,i2,a)') 'tdef ', irec/nlevz,
C     $    ' linear ',cmonth(imnth),iyr1, max(ndmnth,1),'mo'
      write(3,'(a,i2)') 'vars ',1
         if(ivar.eq. 101) then
            write(3,'(a,1x,a)') 'zeta',
     $        '0 0 Sea level'
         elseif(ivar.eq. 103) then
            write(3,'(a,1x,i4,a)') 'u',
     $           ilev,' 0 zonal surface velocity'
         elseif(ivar.eq. 104) then
           write(3,'(a,1x,i4,a)') 'v',
     $           ilev,' 0 meridional surface velocity'
         elseif(ivar.eq. 105) then
            write(3,'(a,1x,i4,a)') 's',
     $           ilev,' 0 surface salinity'
         elseif(ivar.eq. 113) then
            write(3,'(a,1x,a)') 'sicth',
     $           '0 0 ice thickness'
         elseif(ivar.eq. 115) then
            write(3,'(a,1x,a)') 'sicom',
     $           '0 0 ice compactness'
         elseif(ivar.eq. 52) then
            write(3,'(a,1x,a)') 'taux',
     $           '0 0 zonal wind stress'
         elseif(ivar.eq. 53) then
            write(3,'(a,1x,a)') 'tauy',
     $           '0 0 meridional wind stress'
         elseif(ivar.eq. 162) then
            write(3,'(a,1x,i4,a)') 'T',
     $        ilev,' 0 surface temperature'
         else
           write(3,'(a,i1,i3.3,1x,i5,a,3i6)') 'cd',0,abs(ivar),
     $              min(nlevz,ilev),' 0 ',idat1,ivar,ilev
        endif
      write(3,'(a)') 'endvars'
      return
      end

c =====================================================================
      subroutine  hp2gad(data,ivar,len,nx2,ny2,imin,jmax)

C$$$ Convert hope odd/even fields in one(!) EXTRA file to a grads
C$$$ file for plotting

c for HOPE T21 data: 
      parameter(nxm=74, nym=37)
      parameter(NX2m=2*nxm,NY2m=2*nym,nxnym=nxm*nym,nx2ny2=nx2m*ny2m)

      INTEGER IWETP(nxnym),iweto(nxnym),iwete(nxnym), iwet(nx2ny2)
      REAL    DATAO(NXNYM),DATAE(NXNYM),wet(nx2ny2)
      real data(*)
      logical lfirst
      save nx,ny,imino,imax,jmin,jmaxo,nxny,lfirst
     $     ,iweto,iwete

      data lfirst /.true./
      
      if(lfirst) then
         if(ivar.ge.1000)then
            nx= ivar/1000
            ny=len/nx/2
            ivar=mod(ivar,1000)
         else
            print*,'x dimension not stored in ivar'
            print*,' Enter x dimension'
            read(5,*) nx
            ny=len/nx/2
         endif
         print*,'Using: nx=',nx,' ny=',ny
      
         if(nx.eq.72 .and.ny.eq.35) then
            if(lfirst) print*, 'Assuming World region'
            imin=1
            imax=72
            jmin=3
            jmax=37
         elseif(nx.eq.74 .and.ny.eq.37) then
            if(lfirst) print*, 'Assuming World region'
            imin=1
            imax=74
            jmin=1
            jmax=37
         elseif(nx.eq.36 .and.ny.eq.15) then
            if(lfirst) print*, 'Assuming North Pacific'
            imin=22
            imax=57
            jmin=8
            jmax=22
         else
	    print*,'[imin, imax, jmin, jmax]=?'
            read(5,*) imin,imax,jmin,jmax
	    print*,'Region [imin,imax,jmin,jmax] = [',
     $           imin, ',', imax, ',', jmin, ',', jmax, ']'
         endif

         nxny=nx*ny
         imino=imin
         jmaxo=jmax
CCCC  WETMASK        
      OPEN(8,FILE='MASK',FORM='FORMATTED',STATUS='OLD')
ccc  skip scalar masks in masks file if variable is on vector grid:      
      if ( ivar.eq.103 .or.
     $     ivar.eq.104 .or.
     $     ivar.eq.052 .or.
     $     ivar.eq.053 ) read(8,'(37(74i1/)/)') (idummy,i=1,2*nxm*nym)
      read(8,'(37(74i1/)/)') ((iwetp(i+nxm*(j-1)),i=1,nxm),j=1,nym)
      do i=imin,imax
         do j=jmin,jmax
            iweto(i-imin+1+nx*(j-1-jmin+1))=iwetp(i+nxm*(j-1))
         enddo
      enddo
      read(8,'(37(74i1/)/)') ((iwetp(i+nxm*(j-1)),i=1,nxm),j=1,nym)
      do i=imin,imax
         do j=jmin,jmax
            iwete(i-imin+1+nx*(j-1-jmin+1))=iwetp(i+nxm*(j-1))
         enddo
      enddo
      CLOSE(8)

      endif

      imin=imino
      jmax=jmaxo
      nx2=2*nx
      ny2=2*ny



      do j=1,ny
	 jj=nx*(j-1)
         do i=1,nx
            datao(i+jj)=data(i+jj)
            datae(i+jj)=data(i+jj+nxny)
	 enddo
      enddo
      
      call e2agrd(datao,datae,iweto,iwete,nx,ny,nx,data,iwet)

c$$$      if (lfirst) call wrmask(iwet,nx,ny,nx2,ny2)

      lfirst=.false.
      return
      END

c =====================================================================
      
      subroutine e2agrd(ddo,dde,iwo,iwe,i1,j1,idim,
     $     odat,iowet)
c *********************************************************************
c Interpolates the staggered HOPE odd and even 2d fields 
c ddo dde [dim (i1,j1)] to a regular A-type grid field odat. 
c idim gives the i-dim of the input fields.
c iwo/e are the wet masks. 
c The output field has dimension 4*i1*j1.
c 
c Written by Matthias Muennich
c **********************************************************************
      
      dimension ddo(idim,*),dde(idim,*),iwo(idim,*),iwe(idim,*),
     $     odat(2*i1,*),iowet(2*i1,*)
c$$$      data amiss/ 9.e27/
c$$$      data cmiss/ 9.e99/
      data amiss/ 0./
      data cmiss/ 0./
      
      do j=1,2*j1
         do i=1,2*i1
            odat(i,j)=amiss
         enddo
      enddo
      
      do i=1,i1
         im=i-1
         ip=i+1
         i2=2*i
         i2m=i2-1
         do j=1,j1
            jm=j-1
            jp=j+1
            j2=2*j
            j2m=j2-1
c upper right:
            iowet(i2,j2m)= iwo(i,j)
            odat(i2,j2m)=  ddo(i,j)*iwo(i,j)+(1-iwo(i,j))*amiss
c lower left:
            iowet(i2m,j2)= iwe(i,j)
            odat(i2m,j2)=  dde(i,j)*iwe(i,j)+(1-iwe(i,j))*amiss
c upper left:
            if (im.eq.0) then
               iwoimj=0
            else
               if (ddo(im,j).eq.cmiss) then
                  iwoimj=0
               else
                  iwoimj=iwo(im,j)
               endif
            endif
            if(jm.eq.0) then
               iweijm=0
            else
               if (dde(i,jm).eq.cmiss) then
                  iwoijm=0
               else
                  iweijm=iwe(i,jm)
               endif
            endif
            
c nw:  # wet neighbors:
            nw= iwo(i,j)+iwe(i,j)+iwoimj+iweijm
            niw= iwo(i,j)+iwoimj
c Set wet mask to wet if one or more neighbors are wet:
            iowet(i2m,j2m)= min(nw,1)
c Interpolate if 3 or more neighbors are wet:
            if ((nw+1)/4.eq.1) then
               if(im.eq.0) then
                  ddoimj=0.
               else
                  if(ddo(im,j).eq.cmiss) then
                     ddoimj=0.
                  else
                     ddoimj=ddo(im,j)*iwo(im,j)
                  endif
               endif
               
               if(jm.eq.0) then
                  ddeijm=0.
               else
                  if(dde(i,jm).eq.cmiss) then
                     ddeijm=0.
                  else
                     ddeijm=dde(i,jm)*iwe(i,jm)
                  endif
               endif
               
c$$$               odat(i2m,j2m)=
c$$$     $              ( ddo(i,j)*iwo(i,j)+dde(i,j)*iwe(i,j)
c$$$     $              +ddoimj+ddeijm ) / float(nw)
               odat(i2m,j2m)=
     $              ( ddo(i,j)*iwo(i,j)
     $                  +ddoimj ) / float(niw)
            else
               odat(i2m,j2m)=amiss
            endif
            
c lower right:
            if(jp.eq.j1+1) then
               iwoijp=0
            else
               iwoijp=iwo(i,jp)
            endif
            if(ip.eq.i1+1) then
               iweipj=0
            else
               iweipj=iwe(ip,j)
            endif

            nw= iwo(i,j)+iwe(i,j)+iwoijp+iweipj
            niw= iwe(i,j)+iweipj
c Set wet mask to wet if one or more neighbors are wet:
            iowet(i2,j2)= min(nw,1)
c Interpolate if 3 or more neighbors are wet:
            if ((nw+1)/4.eq.1) then
               if(jp.eq.j1+1) then
                  ddoijp=0.
               else
                  ddoijp=ddo(i,jp)*iwo(i,jp)
               endif
               if(ip.eq.i1+1) then
                  ddeipj=0.
               else
                  ddeipj=dde(ip,j)*iwe(ip,j)
               endif
c$$$               odat(i2,j2)=
c$$$     $              ( ddo(i,j)*iwo(i,j)+dde(i,j)*iwe(i,j)
c$$$     $              +ddoijp+ddeipj ) / float(nw)
               odat(i2,j2)=
     $              ( dde(i,j)*iwe(i,j)
     $                 +ddeipj ) / float(niw)
            else
               odat(i2,j2)=amiss
            endif

         enddo
      enddo
      return
      end

c =====================================================================
c$$$      subroutine wrmask(iwet,nx,ny,nx2,ny2)
c$$$      parameter(nxm=74, nym=37)
c$$$      parameter(NX2m=2*nxm,NY2m=2*nym,nxnym=nxm*nym,nx2ny2=nx2m*ny2m)
c$$$
c$$$      INTEGER iwet(nx2ny2)
c$$$      REAL    wet(nx2ny2)
c$$$
c$$$c real data point: wet=+1
c$$$c interpolated data points: wet=0
c$$$c land points points: wet=-1
c$$$c interpolated land points points: wet=-2
c$$$c                           Wenn nx ungerade:
c$$$      if((nx/2)*2.ne.nx) then
c$$$         do i=1,nx2*ny
c$$$            wet(2*i-1)=float(2*iwet(2*i-1)-1)
c$$$            wet(2*i)  =float(2*iwet(2*i)-2)
c$$$         enddo
c$$$      else
c$$$         do j=1,ny 
c$$$            j2m=2*j-1
c$$$            j2m2=2*j-2
c$$$            do i=1,nx
c$$$c even row:
c$$$c odd column:
c$$$               wet(2*i-1+nx2*j2m)=float(2*iwet(2*i-1+nx2*j2m)-1)
c$$$c even column:
c$$$               wet(2*i+nx2*j2m)  =float(2*iwet(2*i+nx2*j2m)-2)
c$$$c odd row:
c$$$c odd column:
c$$$               wet(2*i-1+nx2*j2m2)=float(2*iwet(2*i-1+nx2*j2m2)-2)
c$$$c even column:
c$$$               wet(2*i+nx2*j2m2)  =float(2*iwet(2*i+nx2*j2m2)-1)
c$$$            enddo
c$$$         enddo
c$$$      endif
c$$$         
c$$$      
c$$$    
c$$$      open(12,file='hopewet.gad',form='unformatted',
c$$$     $        access='direct', recl=4*nx2*ny2)
c$$$      write(12,rec=1) ((wet(i+j*nx2),i=1,nx2),j=ny2-1,0,-1)
c$$$      return
c$$$      end

END
f77 ext2gad.f -o ext2gad.x
ext2gad.x $1 $2 $3 << M
$1
$2
$3
M
#cp fort.2 $2
rm fort.1
rm ext2gad.x ext2gad.f
exit