program et2_subset
c
c     read formatted version of etopo2, select subset, and
c     reduce resolution

c     Rows of data (latitude lines) are read one by
c     one starting at 89 58' North and stepping southward
c     by 2'. There are 5400 latitude lines.  The 89 58'
c     South line is not given, elevation 2810 meters can
c     be assumed.
c     There are 10800 longitude points along each latitude
c     line. The first is 180 West, second 179 58' West, etc.
c
c     Depth and elevations are in meters
 
      
      parameter (ilon=10800,ilat=5400)
      dimension itopo(ilon),alon(ilon),alat(ilat)
      character fileout*80,filein*80
c      dimension idepth(ilon,ilat)

c================ USER SPECIFICATION GO HERE ====================
      filein = '/huron/ftp/datasets/ds759.3/data/etopo2.dat' !input data file
      open(20,file=filein)
      alonl = -80.0  !specify regional subset limits (longitude left)
      alonr = -70.0   ! longitude right
      alatt = -30.0   ! latitude top
      alatb = -40.0   ! latitude bottom
      minute = 60     !decimation in minutes, e.g. 30=.5 degrees
      fileout = 'et2_1x1_chile' !output data file name
      open(21,file=fileout)
C===============================================================      
      
      idec = minute/2
      if(idec*2.ne.minute)then
        print*,'The number of minutes chosen as the decimation factor'
	print*,'is not an even multiple of 2 minutes.  The next lower'
	print*,'decimation interval was chosen.'
      endif
      	
      if(alonl.lt.-180..or.alonr.gt.180.)then
        print*,' longitude specification error +- 180'
	print*,' alonl, alonr : ',alonl,alonr
	stop
      endif	
      if(alatt.gt.90..or.alatb.lt.-90.)then
        print*,' latitude specification error +- 90'
	print*,' alatt, alatb : ',alatt,alatb
	stop
      endif
      
      write(21,30)fileout,alatb,alonl,alatt,alonr,minute,idec*2
30    format('Output from program et2_subset.f, file name : ',a/,
     *       'regional limits :'/,
     *       'SW corner (lat/lon)     :',2f6.0/,
     *       'NE corner (lat/lon)     :',2f6.0/,
     *       'Grid selected (minutes) :',i6/,
     *       'Grid computed (minutes) :',i6)
      write(*,30)fileout,alatb,alonl,alatt,alonr,minute,idec*2 

c     compute the latitude and longitude points
      do i = 1,ilon
        alon(i) = -180.0 + (float(i)-1.0)*2.0/60.0 
      enddo
      do j = 1,ilat
        alat(j) = 90.0 - (float(j) -1.0)*2.0/60.0
      enddo
      		               
      ilonl = nint(alonl + 180.)*30 + 1
      ilonr = nint(alonr + 180.)*30 + 1
      nlon = ilonr - ilonl + 1
      ilatb = nint(90.- alatb)*30 + 1
      ilatt = nint(90.- alatt)*30 + 1
      if(ilonr.eq.ilon+1)ilonr=ilon  !boundary check
      if(ilatb.eq.ilat+1)ilatb=ilat  !boundary check
      nlon = ilonr - ilonl + 1
      nlat = ilatb - ilatt + 1
      ndlon = nlon/idec + 1
      ndlat = nlat/idec + 1
      if(ndlon*idec.gt.ilon)ndlon = ndlon - 1
      if(ndlat*idec.gt.ilat)ndlat = ndlat - 1
      
      
      write(21,35)ndlon, ndlat, (int(alon(i)),
     *   abs(nint(mod(alon(i),1.0)*60.0)),i=ilonl,ilonr,idec)
35    format(2i6,'  nlon and nlat '/,
     *      (15(i4,":",i2,"' ")))
c      write(*,35)ndlon,ndlat, (int(alon(i)),
c     *   abs(nint(mod(alon(i),1.0)*60.0)),i=ilonl,ilonr,idec)      

      jnext = 0
      do 45 j=1,ilat
         if(j.gt.ilatb)go to 60  !stop read data once min. latitude found 
         read(20,'(15i6)')(itopo(i),i=1,ilon)
	 if(j.lt.ilatt.or.j.lt.jnext)go to 45	 
	 write(21,40)int(alat(j)),abs(nint(mod(alat(j),1.0)*60.0)),
     *        (itopo(i),i=ilonl,ilonr,idec)
40       format(i4,":",i2,"'",14i8/,(15i8))	 
c	 write(*,40)int(alat(j)),abs(nint(mod(alat(j),1.0)*60.0)),
c     *        (itopo(i),i=ilonl,ilonr,idec)
	 jnext = j + idec
	  	 
45    enddo

60    print*,'Processing completed: '      
      close(20)
      close(21)     
      stop
      end