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 = 17280, njmx = 21600, nsub = 800000)
      character*80 filein1,filein2,fileout,cneg
      integer*2 idata(2*nimx*njmx)
      dimension gravs(nsub),topos(nsub),weight(nsub)
      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=-80.738
      rltf=80.738
      rln0=0.
      rlnf=360.
      rdln=0.0166666666
      rdlt=0.
      iproj=2
      iland=9999
      idum=9998
c
c  get values from command line
c
      narg = iargc()
      if(narg.ne.4) then
        write(*,'(a)')'  '
        write(*,'(a)')
     &  'Usage: nettleton4 neg grav.gips topo.gips slope.out '
        write(*,'(a)')
     &  '       neg (-1 slopes g<0, +1 slope g>0, 0 all g)'
        write(*,'(a)')'  '
        stop
      else
        call getarg(1,cneg)
        call getarg(2,filein1)
        call getarg(3,filein2)
        call getarg(4,fileout)
        read(cneg,*)ineg
      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
      ni00=50
      do 1000 irow=ni00,nlat,ni00
      i00=max(1,irow-ndi)
      iff=min(nlat,irow+ndi)
      jcol=1
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   calculate the latitude-dependent size of the box
c   this is set for a 0.5 diameter of 160 km and a 2-minute grid.
c
c     ni0=160./(3.66*cost)
      ni0=2*160./(3.66*cost)
      nj0=ni0/2.0
      ndi=2*ni0
      ndj=2*nj0
c
      sig0=ni0*ni0/cost2
      sig0=ni0*ni0
      siglarge=sig0*2.
      sigsmall=sig0/2.
c
      do 2000 jcol=nj0,nlon,nj0
      j00=jcol-ndj
      jff=jcol+ndj
      sig=sig0
 550  continue
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
      itestg=idata(indx)
      itestt=idata(indx+nj)
      iuncr=mod(itestt,2)
c
c  skip unconstrained or land topo data
c
      if(iuncr.eq.0.or.itestt.gt.9998) go to 200
      if(ineg.eq.0.or.itestg.le.ineg.or.itestg.ge.ineg) then
      if(ineg.lt.0.and.itestg.gt.0) go to 200
      if(ineg.gt.0.and.itestg.lt.0.or.itestg.gt.1500.) go to 200
        ns=ns+1
        if(ns.gt.nsub) then
          write(*,'(a)')' increase dimension nsub '
          go to 200
        endif
        jd=jsave-jcol
        id=i-irow
        dist=sqrt(id*id/sig+jd*jd/sig)
        if(dist.le.1.0) then
          wht=0.5*(1.+cos(pi*dist))
        else
          wht=0.
        endif
        weight(ns)=wht
        wsum=wsum+weight(ns)
        gravs(ns)=idata(indx)/10.
        topos(ns)=idata(indx+nj)
      endif
  200 continue
c
c  see if there are a good number of points in the window
c
c     write(*,*)irow,jcol,wsum,sqrt(sig),cost,sqrt(sig)*cost*3.66
      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,sqrt(sig)*cost*3.66
  908 format(2i6,8e13.5,i7,f8.0)
c
c  write out the appropriate area
c
c     rlatt=18.5
c     rlont=190.5
c     if(abs(rlon-rlont).lt.1.0.and.abs(rlat-rlatt).lt.1.0) then
c       do 1500 k=1,ns
c       tpred=gravs(k)*slope 
c       write(lus,*)gravs(k),topos(k),tpred
c1500   continue
c     endif
      endif
 2000 continue
 1000 continue
      close(lus)
      call ihbcl(iin1)
      call ihbcl(iin2)
      stop
      end