program retglobal
c
c  Program to read altimeter data from alt.reg.
c
c  This version is tuned for ERS-1 data obtained on Fast Delivery form
c  Each ERM is stored in a separate file.
c  HERE FOR THE GLOBAL STORED DATA
c
c  Ole B. Andersen. 20 / 9 - 92 , modiefied after P.Knudsen 1990
c  KMS - Denmark.
c
c
      implicit real*8(a-h)
      character*60 rfile
      character*60 name
      character*3 sta
c
c  Altimeter data are stored 'sequentially' on a direct access
c  file.
c  A catalogue and a key in each record describe how data inside
c  subareas are found. The catalogue are stored in the beginning
c  of the file together with:
c
c   nlngth    length of reg.file in records
c   ntotal    total number of data
c   nhead     length of headder in records
c   time1     start time of data
c   time2     end time of data
c   irev1     first revolution number
c   irev2     last revolution number
c   ilat0     lower latitude bound of area
c   ilon0     lower longitude bound of area
c   idlat     size in latitude of subareas
c   idlon     size in longitude of subareas
c   idimi     number of subareas in latitude
c   idimj     number of subareas in longitude
c
c  The catalogue contains for subarea (i,j)
c
c   istore(1,i,j)  first record
c   istore(2,i,j)  last record
c   istore(3,i,j)  number of records in subarea
c
      real*8 time1,time2
      integer*4 irev1,irev2
      integer*4 istore(3,34,36)
c
c  Description of data record:
c
c   time       time in seconds from the beginning of 1985
c   irev       revolution number of track
c   ilat       latitude
c   ilon       longitude
c   issh       sea surface height after orbit and geophys. corr.
c   igeoid     OSU89B geoid heights
c   istd       Standard deviation of issh
c   itide      Ocean tide, =32767 when no information
c   next       number of next record inside sub-area
c              (=0: end)
c
      real*8 time,timelast
      integer*4 ilat,ilon,next,irev0
      integer*2 istd,itide,iswh,idepth,imssh,iwind,issh,igeoid
      integer*2 i2rec(19),icorr,irev
        equivalence (i2rec(1),time),
     .              (i2rec(5),ilat),
     .              (i2rec(7),ilon),
     .              (i2rec(9),issh),
     .              (i2rec(10),imssh),
     .              (i2rec(11),igeoid),
     .              (i2rec(12),idepth),
     .              (i2rec(13),next),
     .              (i2rec(15),istd),
     .              (i2rec(16),irev),
     .              (i2rec(17),itide),
     .              (i2rec(18),iwind),
     .              (i2rec(19),iswh)
 
c
      write(*,*)
      write(*,*) '****************************************************'
      write(*,*)
      write(*,*) 'Retrieve ERS-1 geodetic data from data base.'
      write(*,*)
      write(*,*) 'Notice : Data are available with geoid height as '
      write(*,*) '         provided by ESA and JGM-2 orbits '
      write(*,*)
      write(*,*) ' reg.ID        type '
      write(*,*) ' 1-60          3 day portions '
      write(*,*) ' '
      write(*,*) ' Present output format : '
      write(*,*) ' itrack,latitude,longitude,geoid,ssh,dssh,otide'
      write(*,*) ' --- Format can be changed by user ----'
      write(*,*)
      write(*,*) ' Notice :  geoid model is unknown '
      write(*,*) '           longitude -180<lon<180 '
      write(*,*)
c     write(*,*) 'Enter first and last register.ID '
c     read(*,*) ireg1,ireg2
      ireg1 = 139
      ireg2=143
      ireg3=1
      write(*,*) 'Enter min and max lat and'
      write(*,*) '      min and max lon (deg):'
      read(*,*) alat1,alat2,alon1,alon2
      write(*,*) 'Processing Area: ',alat1,alat2,alon1,alon2
      if (alon1.gt.179.999) then
	 alon1 = alon1-360.0
	 alon2 = alon2-360.0
      end if
      write(*,*) 'Enter NAME of output file'
      read(*,'(A)') name
c     write(*,*) 'Enter STATUS of output file (old or new)'
c     read(*,'(A)') sta
      open(20,file=name,form='formatted')
      write(*,*)
     .'Enter max search criterias (dep from geoid,std (m), SWH(m)'
      write(*,*) 'Usual values (10 m, 0.3 m, 10 m) ' 
      read(*,*) amaxgeo,amaxstd,amaxswh
      imaxgeo = idint(amaxgeo*100)
      imaxstd = idint(amaxstd*100)
      imaxswh = idint(amaxswh*1000)
      write(*,*) 'Enter output longitude range 1=-180,180;2=0,360'
      read(*,*) iout
c
c
c  find area extentions in micro degrees
c
      islat1=alat1*1000000.0
      islat2=alat2*1000000.0
      islon1=alon1*1000000.0
      islon2=alon2*1000000.0
c
      ndat=0
      timelast = 0.0d0
      ireg3 = 1
c
      irevold=0
      do 999 ireg=IREG1,IREG2,IREG3
c
      rfile = 'ers1rads000.reg'
      write(rfile(9:11),'(I3.3)') ireg
      write(*,495) rfile
  495 format(/,' Reg.file: ',a55)
      open(12,file=rfile,access='direct',
     .        recl=38,status='old',err=998)
c
c  read headder and catalogue
c
      jrec=1
      read(12,rec=jrec) nlngth,ntotal,nhead,time1,time2,irev1,irev2
      write(*,*) ' rev.s and ntotal  ',irev1,irev2,ntotal
      jrec=2
      read(12,rec=jrec) ilat0,ilon0,idlat,idlon,idimi,idimj
      do 5 i=1,idimi
      do 5 j=1,idimj
      jrec=jrec+1
      read(12,rec=jrec)(istore(k,i,j),k=1,3)
    5 continue
c
c  Find subareas of interest
c
      ilat1=ilat0+idimi*idlat
      ilon1=ilon0+idimj*idlon

      iarea1=(islat1-ilat0)/idlat+1
      iarea1=max(iarea1,1)
      iarea2=(islat2-ilat0)/idlat+1
      iarea2=min(iarea2,idimi)
      jarea1=(islon1-ilon0)/idlon+1
      jarea1=max(jarea1,1)
      jarea2=(islon2-ilon0)/idlon+1
      jarea2=min(jarea2,idimj)
c
      irec=0
       idat=0
c
c  Start reading in current register file
c
   10 ifirst=ntotal+nhead+3
c
c  Find the subarea containing the first set of observations
c  from the current possition
c
      do 20 i=iarea1,iarea2
      do 20 j=jarea1,jarea2
      if((istore(1,i,j).gt.0).and.(istore(1,i,j).lt.ifirst)) then
        ifirst=istore(1,i,j)
        iarea=i
        jarea=j
        endif
   20 continue
      if(ifirst.eq.(ntotal+nhead+3)) go to 998
      if(ifirst.lt.irec) write(*,*) 'ifirst>irec',ifirst,irec
c
      irec=ifirst
      icount=0
c
c  read track of observations in subarea (iarea,jarea)
c  which ends when next.ne.(irec+1) (goto 41)
c
      do 40 k=1,istore(3,iarea,jarea)
c
      read(12,rec=irec) i2rec
c     if((time-timelast).lt.0.5) go to 30
      icount=icount+1
      if((time-timelast).gt.2.5) then
	timelast=time
	go to 30
	endif
c
      if(ilat.lt.islat1.or.ilat.gt.islat2.or.
     .   ilon.lt.islon1.or.ilon.gt.islon2) go to 30
         if(irev.ne.irevold) then
           irevold=irev
           go to 30
           endif
c        if(istd.gt.imaxstd) go to 30
c        if(iswh.gt.imaxswh) go to 30
         idat=idat+1
         ndat=ndat+1
  
c        jrev = mod(irev,501)
c        if (irev.lt.1) goto 30
c        if (irev.gt.4823) goto 30
c        if (irev.lt.2411) goto 30
	 if (iout.gt.1.and.ilon.lt.0) ilon = ilon+360000000

c	 if (itide.lt.32000) issh = issh + itide 
c	 if (itide2.lt.32000) issh = issh - itide2

         ialt = issh + igeoid*10
c	 goto 922
	 dh = 0.0
          if(abs(issh).gt.imaxgeo) go to 30
          if(abs(issh-issh0).gt.450) go to 30
c         if(abs(issh-issh0).gt.150) go to 30
c         if(istd-abs(issh-issh0)*0.25.gt.60) go to 30
c         if(istd-abs(issh-issh0)*0.25.gt.20) go to 30

	 if(irev.ne.irevold) icut=0
	 if(irev.eq.irevold.and.(time-timelast).gt.5.0) icut=icut+1
	 if(istd.lt.10) istd = 10
c 922     continue
         write(20,444) (irev+icut*10000),(ilat*1.0d-6),
     .    (ilon*1.0d-6),dh,(ialt*1.0d-3),(istd*1.0d-3),time
  444    format(i6,2f11.5,f4.1,2f9.3,f12.1)
c
   30 continue
	 irevold=irev
	 timelast=time
	 issh0=issh
c
      if(next.ne.(irec+1)) go to 41
      irec=next
   40 continue
   41 continue
c
c  Update catalogue
c
      istore(1,iarea,jarea)=next
      istore(3,iarea,jarea)=istore(3,iarea,jarea)-icount
      go to 10
  998 write(*,510) idat
  510 format(50x,'number of data:',i6)
      close(12)
  999 continue
      write(*,520) ndat
  520 format(/,44x,' Total number of data:',i6)
      close(20)
      end