program nettleton4
c
c  program to determine the mapping from gravity to topography
c  in sub regions.
c
c  The standard good window contains 60 points.
c
c  If sumw is less than 42 the window radius is increased
c  by a factor of square root of 2.
c
c  If sumw is greater than 170 the window radius is decreased
c  by a factor of square root of 2.
c 
c
      parameter (nimx = 6336, njmx = 10800, nsub = 40000)
      parameter (nrow = 64, ncol = 108)
      character*80 filein1,filein2,fileout
      integer*2 idata(2*nimx*njmx)
      dimension gravs(nsub),topos(nsub),weight(nsub)
      dimension w(-nimx/nrow:nimx/nrow,-njmx/ncol:njmx/ncol)
      common/grid/nlat,nlon,rdlt,rdln,iproj
      common/bounds/rlt0,rltf,rln0,rlnf,rlnm
      data pi,rad,gamma/3.1415927,.01745329,6.67e-11/
      data lus /20/
c
c  setup the grid parameters
c
      nlon=njmx
      nj=njmx
      nlat=nimx
      ni=nimx
      rlt0=-72.006
      rltf=72.006
      rln0=0.
      rlnf=360.
      rdln=0.033333333
      rdlt=0.
      iproj=2
      iland=9999
      idum=9998
c
c  get values from command line
c
      narg = iargc()
      if(narg.ne.3) then
        write(*,'(a)')'  '
        write(*,'(a)')
     &  'Usage: nettleton4 grav.gips topo.gips slope.out '
        write(*,'(a)')'  '
        stop
      else
        call getarg(1,filein1)
        call getarg(2,filein2)
        call getarg(3,fileout)
      endif
c
c  add string terminator to filenames and decode sigma
c
      nc=index(filein1,' ')
      filein1(nc:nc)='\0'
c
c  open the gravity file and read header info
c
      iin1=ihread(filein1,"/mgg/gips/tables/cprofile")
      if (iin1.eq.0) then
        write(*,*) '****header incorrectly formatted '
        stop
      endif
c
c  check the file size
c
      ig=ihget("xnum",njg)
      ig=ihget("ynum",nig)
      if(njg.ne.nj.or.nig.ne.ni) then
       write(*,'(a)') 'incorrect dimension for grav file'
       stop
      endif
c
      nc=index(filein2,' ')
      filein2(nc:nc)='\0'
c 
c  open the topography file and read header info 
c 
      iin2=ihread(filein2,"/mgg/gips/tables/cprofile") 
      if (iin2.eq.0) then 
        write(*,*) '****header incorrectly formatted ' 
        stop 
      endif 
c 
c  check the file size 
c 
      ig=ihget("xnum",njt) 
      ig=ihget("ynum",nit) 
      if(njt.ne.nj.or.nit.ne.ni) then
       write(*,'(a)') 'incorrect dimension for topo file'
       stop
      endif
c
c  open the output file
c
      open(unit=lus,file=fileout,err=9000)
      go to 9001
 9000 write(*,*)' could not open ',name
 9001 continue
c
c  read the gravity file into row i and the
c  topography file into row i+1 this will place them
c  close together in memory.
c
      do 100 i=1,ni
      istart=2*nj*(i-1)+1
      in = ihbin(idata(istart),iin1,2*nj)
      if(in.lt.0) then
       write(*,*)' could not read, rec# ',filein1,i
       stop
      endif
      istart=istart+nj
      in = ihbin(idata(istart),iin2,2*nj)
      if(in.lt.0) then
       write(*,*)' could not read, rec# ',filein2,i
       stop
      endif
  100 continue
c
c  calculate these coefficients for overlapping areas
c
      ndi=ni/nrow
      ni0=ndi/2
      ndj=nj/ncol
      nj0=ndj/2
c
      do 1000 irow=ni0,nlat,ni0
      i00=max(1,irow-ndi)
      iff=min(nlat,irow+ndi)
      jcol=ndj
c     call proj(rlat,rlon,irow,jcol,-1)
      call mercator(rlat,rlon,irow,jcol,-1)
      cost=cos(abs(rlat*rad))
      cost2=cost*cost
c
c   fill the weight array
c
      sig0=ni0*ni0/cost2
      siglarge=sig0*2.
      sigsmall=sig0/2.
      sig=sig0
 550  continue
      do 700 i=-ndi,ndi
      do 600 j=-ndj,ndj
      dist=sqrt(i*i/sig+j*j/sig)
      if(dist.le.1.0) then
        w(i,j)=0.5*(1.+cos(pi*dist))
      else
        w(i,j)=0.
      endif
 600  continue
 700  continue
c
      do 2000 jcol=nj0,nlon,nj0
      j00=jcol-ndj
      jff=jcol+ndj
c
c  loop over the area and get the data
c
      ns=0
      wsum=0.
      do 200 i=i00,iff
      ioff=2*nj*(i-1)
      do 200 j=j00,jff
c
c wrap the data arrays around Grenwich
c
      jsave=j
      jdata=1+mod(j-1+nj,nj)
      indx=jdata+ioff
c
c  odd point is constrained by ship data
c  even point has no constraint
c
      itest=idata(indx+nj)
      iuncr=mod(itest,2)
      if(iuncr.eq.0.or.itest.ge.0) go to 200
      ns=ns+1
      if(ns.gt.nsub) then
        write(*,'(a)')' increase dimension nsub '
        go to 200
      endif
      weight(ns)=w(i-irow,jsave-jcol)
      wsum=wsum+weight(ns)
      gravs(ns)=idata(indx)/10.
      topos(ns)=idata(indx+nj)
  200 continue
c
c  see if there are a good number of points in the window
c
      if(wsum.lt.42.and.sig.eq.sig0) then
       sig=siglarge
       go to 550
      endif
      if(wsum.gt.170.and.sig.eq.sig0) then
       sig=sigsmall
       go to 550
      endif
c
c   compute the nettleton slope
c
      sumw=0.
      gbar=0.
      tbar=0.
      r=0.
      prob=0.
      slope=-100.
      if(ns.gt.6) then
        int1=nslope3(gravs,topos,weight,ns,sumw,
     +               gbar,sigg,tbar,sigt,r,prob)
      if(sigg.le.0..or.r.eq.0..or.sumw.lt.2.) then
        slope=-100.
      else
        slope=sign(sigt/sigg,r)
      endif
c
c  compute the lat and lon
c
c     call proj(rlat,rlon,irow,jcol,-1)
      call mercator(rlat,rlon,irow,jcol,-1)
      irowbot=ni-irow+1
      write(lus,908)jcol,irowbot,rlon,rlat,sumw,sigg,
     +              sigt,slope,r,prob,ns,sig
  908 format(2i6,8e13.5,i7,f8.0)
c
c  write out the appropriate area
c
c     if(abs(rlon-rlont).lt.0.5.and.abs(rlat-rlatt).lt.0.5) then
c       do 1500 k=1,ns
c       tpred=gravs(k)*slope 
c       write(*,*)gravs(k),topos(k),tpred
c1500   continue
c     endif
      endif
 2000 continue
 1000 continue
      close(lus)
      call ihbcl(iin1)
      call ihbcl(iin2)
      stop
      end