c gi4 - gi3 applied to a large gravity grid through windowing
c
c written by Jeff Phillips, 7/96
c
      dimension grav(500,30),surf(500,30),cgrav(500,30)
      dimension wgt1(500,30)
      dimension g(30,30),cg(30,30),s(30,30),wgt(30,30)
c      dimension dv(500)
      character ifile*50, id*56, pgm*8, ans*1
      data ddval/1.e30/,dval/1.e38/
c
c initialize
c
      nover=10
      call pfinit('gi4')
      call askin()
c
c get input parameters
c
      ifile=' '
1     call askc('Enter gridded gravity file: ',ifile,icode)
      if(icode.lt.0) stop
      call gopen(20,ifile,'old','read',ierr)
      if(ierr.ne.0) stop 'error opening file'
      call gheader('r',20,id,nc,nr,xo,dx,yo,dy,ierr)
      if(ierr.ne.0) stop 'error reading grid header'
      if(nc.gt.500) stop 'too many columns'
      if(dx.ne.dy) stop 'dx ne dy'
      gmin=1.e37
      gmax=-1.e37
      do 20 j=1,nr
      call grow('r',20,j,grav(1,1),nc,ierr)
      do 20 i=1,nc
      if(grav(i,1).ge.ddval) go to 20
      if(grav(i,1).gt.gmax) gmax=grav(i,1)
      if(grav(i,1).lt.gmin) gmin=grav(i,1)
20    continue
2      call askf4('Enter density contrast: ',rho,icode)
      if(icode.eq.-2) go to 1
      if(rho.gt.0.) const1=-gmin
      if(rho.lt.0.) const1=-gmax
      if(rho.eq.0.) go to 2
      print*,' gravity maximum = ',gmax
      print*,' gravity minimum = ',gmin
3     call askf4('Enter constant to be added to the gravity: ',
     1const1,icode)
      if(icode.eq.-2) go to 2
4     call askf4('Enter reference surface depth: ',zs1,icode)
      if(icode.eq.-2) go to 3
5     call askc('Is this the top or bottom surface? (t or b): ',
     1ans,icode)
      if(icode.lt.0) go to 4
      iopt1=0
      if(ans.eq.'t'.or.ans.eq.'T') iopt1=1
      if(ans.eq.'b'.or.ans.eq.'B') iopt1=2
      nmax=20
6     call aski4('Enter number of iterations: ',nmax,icode)
      if(icode.eq.-2) go to 5
      eps=1.
7     call askf4('Enter desired error: ',eps,icode)
      if(icode.eq.-2) go to 6
      ifile=' '
8     call askc('Enter output surface file: ',ifile,icode)
      if(icode.lt.0) go to 7
      call gopen(21,ifile,'new','readwrite',ierr)
      if(ierr.ne.0) stop 'error opening file'
      id=' '
9     call askalt
      call askc('Title:',id,icode)
      call askalt
      if(icode.eq.-2) go to 8
      call gheader('w',21,id,nc,nr,xo,dx,yo,dy,ierr)
      if(ierr.ne.0) stop 'error writing grid header'
      ifile=' '
10    call askc('Enter output calculated gravity file: ',ifile,icode)
      if(icode.lt.0) go to 9
      call gopen(22,ifile,'new','readwrite',ierr)
      if(ierr.ne.0) stop 'error opening file'
      id=' '
      call askalt
      call askc('Title:',id,icode)
      if(icode.eq.-2) go to 10
      call gheader('w',22,id,nc,nr,xo,dx,yo,dy,ierr)
      if(ierr.ne.0) stop 'error writing grid header'
      call gopen(23,'gi4.wgt','unknown','write',ierr)
      if(ierr.ne.0) stop 'error opening file'
      id='gi4 weights '
      call gheader('w',23,id,nc,nr,xo,dx,yo,dy,ierr)
      if(ierr.ne.0) stop 'error writing grid header'
c
c assign weights
c
      do 15 j=1,nover
      do 15 i=1,30
      if (j.le.i.and.j.le.(31-i)) then
c        wgt(i,j)=float(j)/float(nover+1)
        wgt(i,j)=float(j-1)/float(nover-1)
      else if (i.le.nover) then
c        wgt(i,j)=float(i)/float(nover+1)
        wgt(i,j)=float(i-1)/float(nover-1)
      else
c        wgt(i,j)=float(31-i)/float(nover+1)
        wgt(i,j)=float(30-i)/float(nover-1)
      endif
15    continue
      do 16 j=nover+1,30-nover
      do 16 i=1,30
      if(i.le.nover) then
c        wgt(i,j)=float(i)/float(nover+1)
        wgt(i,j)=float(i-1)/float(nover-1)
      else if(i.gt.(30-nover)) then
c        wgt(i,j)=float(31-i)/float(nover+1)
        wgt(i,j)=float(30-i)/float(nover-1)
      else
        wgt(i,j)=1.0
      endif
16    continue
      do 17 j=31-nover,30
      do 17 i=1,30
      if(j.ge.i.and.j.ge.(31-i)) then
c        wgt(i,j)=float(31-j)/float(nover+1)
        wgt(i,j)=float(30-j)/float(nover-1)
      else if(i.le.nover) then
c        wgt(i,j)=float(i)/float(nover+1)
        wgt(i,j)=float(i-1)/float(nover-1)
      else
c        wgt(i,j)=float(31-i)/float(nover+1)
        wgt(i,j)=float(30-i)/float(nover-1)
      endif
17    continue
c      do 18 j=1,nover
c      do 18 i=1,nover
      do 18 j=2,nover-1
      do 18 i=2,nover-1
c      y=abs(float(j)-float(nover+1)/2.)/(float(nover+1)/2.)
c      x=abs(float(i)-float(nover+1)/2.)/(float(nover+1)/2.)
      y=abs(float(j-1)-float(nover-1)/2.)/(float(nover-1)/2.)
      x=abs(float(i-1)-float(nover-1)/2.)/(float(nover-1)/2.)
      den=2. - max(x,y)
      wgt(i,j)=wgt(i,j)/den
      wgt(31-i,j)=wgt(31-i,j)/den
      wgt(i,31-j)=wgt(i,31-j)/den
      wgt(31-i,31-j)=wgt(31-i,31-j)/den
c      wgt(30-i,j)=wgt(30-i,j)/den
c      wgt(i,30-j)=wgt(i,30-j)/den
c      wgt(30-i,30-j)=wgt(30-i,30-j)/den
18    continue
c
c initialize output arrays
c
      do 21 j=1,30
      do 21 i=1,nc
      surf(i,j)=0.
      cgrav(i,j)=0.
      wgt1(i,j)=0.
21    continue
      rewind(20)
      read(20) id
      nrowb=0
      nrow=30
      if(nr.lt.30) nrow=nr
c
c  top of loop through rows
c
22    continue
      do 30 j=1,nrow
      call grow('r',20,j,grav(1,j),nc,ierr)
30    continue
      ncoff=0
      ncol=30
      if(nc.lt.30) ncol=nc
c
c  top of loop through columns
c
32    continue
      do 35 j=1,nrow
      do 35 i=1,ncol
      if(grav(i+ncoff,j).ge.ddval) go to 41
      g(i,j)=grav(i+ncoff,j)+const1
      if(rho.lt.0.) then
        if(g(i,j).gt.0.0) g(i,j)=0.0
      else
        if(g(i,j).lt.0.0) g(i,j)=0.0
      endif
35    continue
      print*,'processing rows ',nrowb+1,' to ',nrowb+nrow
      print*,'processing cols ',ncoff+1,' to ',ncoff+ncol
      call gi3(g,cg,s,nrow,ncol,zs1,rho,iopt1,dx,nmax,eps)
      do 40 j=1,nrow
      do 40 i=1,ncol
      cgrav(i+ncoff,j)=cgrav(i+ncoff,j)+wgt(i,j)*cg(i,j)
      surf(i+ncoff,j)=surf(i+ncoff,j)+wgt(i,j)*s(i,j)
      wgt1(i+ncoff,j)=wgt1(i+ncoff,j)+wgt(i,j)
40    continue
      go to 42
41    continue
        do 33 j=1,nrow
        do 33 i=1,ncol
        cgrav(i+ncoff,j)=0.
        surf(i+ncoff,j)=0.
33      continue
42     if(ncol+ncoff.ne.nc) then
         ncoff=ncoff+30-nover
         if(ncol+ncoff.gt.nc) ncol=nc-ncoff
         go to 32
c
c  end of loop through columns
c
       endif
      print*,'writing rows ',nrowb+1,' to ',nrowb+nrow-nover
      do 45 j=1,nrow-nover
      do 455 i=1,nc
      if(wgt1(i,j).lt.0.25.or.i.eq.nc.or.j+nrowb.eq.nr) then
        cgrav(i,j)=dval
        surf(i,j)=dval
      else
        cgrav(i,j)=cgrav(i,j)/wgt1(i,j)-const1
        surf(i,j)=surf(i,j)/wgt1(i,j)
      endif
455   continue
c      do 456 i=1,2
c      cgrav(i,j)=dval
c      cgrav(nc-i+1,j)=dval
c      surf(i,j)=dval
c      surf(nc-i+1,j)=dval
c456   continue
c      if(nrowb.eq.0.and.j.le.2) then
c        do 43 i=1,nc
c        surf(i,j)=dval
c        cgrav(i,j)=dval
c43      continue
c      endif
      call grow('w',21,j,surf(1,j),nc,ierr)
      call grow('w',22,j,cgrav(1,j),nc,ierr)
      call grow('w',23,j,wgt1(1,j),nc,ierr)
45    continue
      if(nrowb+nrow.ne.nr) then
        do 46 j=1,nover
        do 46 i=1,nc
        surf(i,j)=surf(i,j+2*nover)
        cgrav(i,j)=cgrav(i,j+2*nover)
        wgt1(i,j)=wgt1(i,j+2*nover)
46      continue
        do 47 j=nover+1,30
        do 47 i=1,nc
        surf(i,j)=0.
        cgrav(i,j)=0.
        wgt1(i,j)=0.
47      continue
        nrowb=nrowb+30-nover
        if(nrowb+nrow.gt.nr) nrow=nr-nrowb
        rewind(20)
        read(20) id
        do 56 j=1,nrowb
        read(20) dlt
56      continue
        go to 22
      endif
c
c  end of loop through rows
c
      print*,'writing rows ',nrowb+nrow-nover+1,' to ',nrowb+nrow
      do 50 j=nrow-nover+1,nrow
      do 505 i=1,nc
      if(wgt1(i,j).lt.0.25.or.i.eq.nc.or.j+nrowb.eq.nr) then
        cgrav(i,j)=dval
        surf(i,j)=dval
      else
        cgrav(i,j)=cgrav(i,j)/wgt1(i,j)-const1
        surf(i,j)=surf(i,j)/wgt1(i,j)
      endif
505   continue
c      do 506 i=1,2
c      cgrav(i,j)=dval
c      cgrav(nc-i+1,j)=dval
c      surf(i,j)=dval
c      surf(nc-i+1,j)=dval
c506   continue
c      if(j.ge.nrow-1) then
c        do 48 i=1,nc
c        surf(i,j)=dval
c        cgrav(i,j)=dval
c48      continue
c      endif
      call grow('w',21,j,surf(1,j),nc,ierr)
      call grow('w',22,j,cgrav(1,j),nc,ierr)
      call grow('w',23,j,wgt1(1,j),nc,ierr)
50    continue
      stop
      end
c
c  gi3.for
c
c  Iterative 3-D inversion of gravity data based on method of
c    Cordell and Henderson, 1968, Geophysics v. 33, p. 596-601.
c
c  Driver for iterative 3-d gravity subroutine 'GRVINV'. Written 
c    by Cordell ca. 1977 and subsequently modified by Cordell,
c    R.N. Godson, et al.  This version by Cordell, Oct 1988.
c
c  To compile: use /Gt1800
c
c       arrays to separate data segments
      subroutine gi3(g,t,zs,nrow,ncol,zs1,rho,iopt1,rdx,nmax,eps)
      implicit double precision (d)
      dimension g(30,30),dgcalc(30,30),t(30,30),zs(30,30),
     &da1(32,32),da2(32,32),id(14),pgm(2),bigst(50),ib(50),
     &jb(50),rmsn(50)

c      equivalence (imax,ncol),(jmax,nrow),(iopt3,mode)
c      common /parms/ w,zs1,rho,nmax,nstar,eps,mode,iopt1,iopt5
      imax=ncol
      jmax=nrow
c      rho=0.0
      nstar=0
c      nmax=10
c      eps=0.1
      mode=0
      iopt2=1
      iopt3=mode
      iopt4=0
      iopt5=0
      iopt6=0
c      zs1=0
cc  read command and data files
c        print*
c      write(6,1)
c1     format(1x,'Enter command file name:',/,' *',$)
c      read(5,2) cfn
c2     format(a50)
c       open(10,file=cfn,status='old')
cc      read (10,parms,end=99)
c        call namemc(10)
c      close(10)
c      write(6,40)zs1,rho,mode,iopt1
c40    format(1x,'Following are the basic input parameters:',/,1x,'zs1 ='
c     & ,e12.6,10x,'rho =',e12.6,10x,/,1x,'mode =',i2,10x,
c     & 'iopt1 =',i2/)
c      if(nmax.eq.0)nmax=5
c      write(6,41) nmax,nstar,eps,iopt5
c41    format(1x,'Following are input parameters which may have been
c     & assigned by default:',/,1x,'nmax =',i2,10x,'nstar =',i2,
c     & 10x,'eps =',e12.6,10x,'iopt5 =',i2,/,1x,'Want to continue ?',
c     & 1x,$)
c      read(5,42) answer
c42    format(a1)
c      if(answer.eq.'y'.or.answer.eq.'Y') go to 44
c55    write(6,43)
c43    format(1x,'Return to command level to revise parameters.')
c      go to 999
c44    continue
c      if(mode.eq.2) go to 56
c      if(mode.eq.1.or.mode.eq.0) go to 50
c      iferr=6
c      go to 90
c50    write(6,51)
c51    format(/1x,'Enter gravity data file name:',/' *',$)
c      read(5,2) gdfn
c      open(10,file=gdfn,status='old',
c     & form='unformatted')
c      read(10)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
c        if(ncol.gt.50.or.nrow.gt.50) then
c        write(6,52)
c52    format(' Error: ncol or nrow > 50.')
c        go to 999
c        endif
c        if(rdx.ne.rdy) then
c      print *,'Error: dx must equal dy.'
c        go to 999
c        endif
c        if(nz.ne.1) then
c      print *,'Error: nz must equal 1.'
c        go to 999
c        endif
        w=rdx
c      do 61 j=1,nrow
c      read(10) rdum,(g(i,j),i=1,ncol)
c        do 662 k=1,ncol
c        if(g(k,j).ge.1.e+30) then
c        print*,'Can''t handle dvals.'
c        go to 999
c        endif
c662     continue
c61      continue
c      close(10)
c56    if(mode.ne.1.and.mode.ne.2) go to 57
c53    write(6,54)
c54    format(1x,'Enter model array file name:',/' *',$)
c      read(5,2) mafn
c       open(10,file=mafn,status='old',
c     & form='unformatted')
c      read(10)id,pgm,mcol,mrow,mz,xx0,xdx,yy0,ydy
c        if(mode.eq.2)then
c          rdx=xdx
c          rdy=ydy
c          ncol=mcol
c          nrow=mrow
c          nz=mz
c          x0=xx0
c          y0=yy0
c          if(ncol.gt.50.or.nrow.gt.50) then
c            write(6,52)
c            go to 999
c          endif
c          if(rdx.ne.rdy) then
c            print*,'Error: dx must equal dy.'
c            go to 999
c          endif
c          if(nz.ne.1) then
c            print*,'Error: nz must equal 1.'
c            go to 999
c          endif
c          w=rdx
c          go to 661
c        endif
c        if(mcol.ne.ncol.or.mrow.ne.nrow.or.mz.ne.nz) go to 660
c        if(xdx.ne.rdx.or.ydy.ne.rdy) go to 660
c        slop=.00001
c        if(xx0.gt.x0+slop.or.xx0.lt.x0-slop) go to 660
c        if(yy0.gt.y0+slop.or.yy0.lt.y0-slop) go to 660
c        go to 661
c660     print*,'Grid specifications don''t match'
c        go to 999
c661   do 62 j=1,nrow
c62    read(10) rdum,(t(i,j),i=1,ncol)
c      close(10)
c57    if(iopt5.eq.0) go to 60
c      write(6,58)
c58    format(1x,'Enter reference surface array file name:',/' *',$)
c      read(5,2) rsafn
c      open(10,file=rsafn,status='old',
c     & form='unformatted')
c      read(10)id,pgm,mcol,mrow,mz,xx0,xdx,yy0,ydy
c        if(mcol.ne.ncol.or.mrow.ne.nrow.or.mz.ne.nz) go to 664
c        if(xdx.ne.rdx.or.ydy.ne.rdy) go to 664
c        slop=.00001
c        if(xx0.gt.x0+slop.or.xx0.lt.x0-slop) go to 664
c        if(yy0.gt.y0+slop.or.yy0.lt.y0-slop) go to 664
c        go to 663
c664     print*,'Grid specifications don''t match'
c        go to 999
c663    do 63 j=1,nrow
c63    read(10) rdum,(zs(i,j),i=1,ncol)
c      close(10)
c60    continue
      gmin=1.e37
      gmax=-1.e37
      do 120 j=1,nrow
      do 120 i=1,ncol
c      if(g(i,j).ge.ddval) stop 'dvals found'
      if(g(i,j).gt.gmax) gmax=g(i,j)
      if(g(i,j).lt.gmin) gmin=g(i,j)
120   continue
c      print*,'local gmin, gmax = ',gmin,gmax
      if(rho.gt.0.) const0=-gmin
      if(rho.lt.0.) const0=-gmax
c      print*,'adding ',const0,' to local anomaly'
      do 35 j=1,nrow
      do 35 i=1,ncol
      g(i,j)=g(i,j)+const0
35    continue
      h=-const0/(41.927696*rho)
      if(iopt1.eq.1) zs1loc=zs1+h
      if(iopt1.eq.2) zs1loc=max(0.,zs1-h)
c      print*,'changing reference level to ',zs1loc,
c     1' using a slab of thickness ',h
c
      write(6,64)
64    format(1x,'iteration',3x,'largest error',3x,'at point',
     15x,'rms error')
c64    format(1x,'There may be a long pause;',
c     & ' wait for Computation completed')
c  call the GRVINV  subroutine
      call GRVINV(imax,jmax,zs1loc,rho,w,nmax,nstar,eps,iopt1,iopt2,
     & mode,iopt5,g,dgcalc,t,zs,da1,da2,bigst,ib,jb,rmsn,iferr)
90    if(iferr.eq.0) go to 20
      if(iferr.eq.1) write(6,10)
10    format(' Error: the reference plane or surface is not
     & properly defined')
      if(iferr.eq.3) write(6,12)
12    format(' Error in definition of grid width')
      if(iferr.eq.4) write(6,13)
13    format(' Error: gravity and density contrast',/
     & ' must have the same algebraic sign')
      if(iferr.eq.5) write(6,14)
14    format(' Error in definition of units')
      if(iferr.eq.6) write(6,15)
15    format(' Error: parameter mode must be 0,1,2, or undefined
     & in parms list')
	if(iferr.eq.7) write(6,116)
116	format(' Error:  base or middle of prisms cannot be
     $ referenced to 0 depth.'/'  Check zs1 and iopt1')
      go to 999
20    continue
c  output
c30    write(6,31)
c31    format(1x,'Computation completed.')
c      if(iopt3.eq.2)go to 32
c      write(6,21) nmax
c21    format(1x,'The program has completed ',i3,' iterations.')
c      do 22 i=1,nmax
c22    write(6,23)i,bigst(i),ib(i),jb(i),rmsn(i)
c23    format(1x,'At iteration ',i3, ' the largest error is ',
c     & e12.4,' at point '
c     & ,i3,',',i2,' ;'/1x,'    rms error =',e12.4,'.')
c32    write(6,100)
c100   format(1x,'Want to see some profiles? [n]')
c      read(5,42) answer
c      if(answer.ne.'y'.and.answer.ne.'Y')go to 120
c101   write(6,102)
c102   format(1x,'Row or column?')
c      read(5,42)answer
c      if(answer.eq.'r'.or.answer.eq.'R')go to 103
c      if(answer.eq.'c'.or.answer.eq.'C')go to 108
c      go to 101
c103   write(6,104)
c104   format(1x,'Enter row number.')
c      read(5,*) jrow
c      write(6,140)
c140   format(5x,'Gravity   Calc Grav   Thickness')
c      do 106 i=1,ncol
c106   write(6,107)i,g(i,jrow),dgcalc(i,jrow),t(i,jrow)
c107   format(1x,i2,3f10.3)
c      go to 111
c108   write(6,109)
c109   format(1x,'Enter column number.')
c      read(5,*)icol
c      write(6,140)
c      do 110 j=1,nrow
c110   write(6,107)j,g(icol,j),dgcalc(icol,j),t(icol,j)
c111   write(6,112)
c112   format(1x,'More profiles?')
c      read(5,42) answer
c      if(answer.eq.'y'.or.answer.eq.'Y')go to 101
c120   continue
cc  output standard files
ccaution: note that in this section array t is redefined.
c      if(iopt3.eq.2)go to 130
c121   write(6,123)
c123   format(1x,'Enter output thickness file name (or hit return)')
c      read(5,2)omfn
c        if(omfn.eq.' ') go to 126
c      write(6,127)
c127   format(1x,'Enter title (record id of standard file).')
c      read(5,128)(id(i),i=1,14)
c128   format(14a4)
c      open(11,file=omfn,status='unknown',
c     & form='unformatted')
c      write(11)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
c      do 124 j=1,nrow
c      y=j
c124   write(11)y,(t(i,j),i=1,ncol)
c      close(11)
c      write(6,125)omfn
c125   format(1x,'Segment ',a50,/' contains derived thickness',
c     & ' in standard fmt.')
c126   continue
c141   write(6,143)
c143   format(1x,'Enter derived surface file name (or hit return).')
c      read(5,2)dsfn
c        if(dsfn.eq.' ') go to 150
c      write(6,127)
c      read(5,128)(id(i),i=1,14)
c      open(11,file=dsfn,status='unknown',
c     & form='unformatted')
c      write(11)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
c      do 144 j=1,nrow
      do 146 j=1,nrow
      do 146 i=1,ncol
146   zs(i,j)=sngl(da2(i,j))-h
c146   zs(i,j)=sngl(da2(i,j))
c      y=j
c144   write(11)y,(t(k,j),k=1,ncol)
c      close(11)
c      write(6,145)dsfn
c145   format(1x,'Segment ',a50,/' contains derived surface',
c     & ' standard file.')
c150   continue
c130   write(6,132)
c132   format(1x,'Enter calculated gravity file name (or hit return).')
c      read(5,2)cgfn
c        if(cgfn.eq.' ') go to 135
c      write(6,127)
c      read(5,128)(id(i),i=1,14)
c      open(11,file=cgfn,status='unknown',
c     & form='unformatted')
c      write(11)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
c      do 136 j=1,nrow
c      print*,'subtracting ',const0,' from calculated anomaly'
      do 133 j=1,nrow
      do 133 i=1,ncol
133   t(i,j)=sngl(dgcalc(i,j))-const0
c      y =j
c136   write(11)y,(t(k,j),k=1,ncol)
c      close(11)
c      write(6,134)cgfn
c134   format(1x,'Segment ',a50,/' contains calculated gravity array in
c     & standard file format.')
c135   continue
c151   write(6,153)
c153   format(1x,'Enter gravity residual file name (or hit return).')
c      read(5,2)grfn
c        if(grfn.eq.' ') go to 160
c      write(6,127)
c      read(5,128)(id(i),i=1,14)
c      open(11,file=grfn,status='unknown',
c     & form='unformatted')
c      write(11)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
c      do 155 j=1,nrow
c      do 154 i=1,ncol
c154   t(i,j)=g(i,j)-sngl(dgcalc(i,j))
c      y=j
c155   write(11)y,(t(k,j),k=1,ncol)
c      close(11)
c      write(6,156)grfn
c156   format(1x,'Segment ',a50,/' contains residual',
c     & ' gravity standard file.')
c160   continue
c      go to 999
c99    continue
c      write(6,900)
c900   format(' Abnormal termination due to problem',
c     & ' in command file.')
999   return
      end
      subroutine GRVINV(imax,jmax,zs1,rho,w,nmax,nstar,eps,iopt1,iopt2,
     & iopt3,iopt5,g,dgcalc,t,zs,da1,da2,bigst,ib,jb,rmsn,iferr)
c  ITERATIVE THREE-DIMENSIONAL SOLUTION OF GRAVITY ANOMALY DATA.
c  LIN CORDELL, U.S. GEOLOGICAL SURVEY, WASHINGTON D.C., NOV. 1968.
c  RECAST AS A SUBROUTINE AUG 1977.
          implicit double precision (d)
        dimension g(30,1),dgcalc(30,1),t(30,1),zs(30,1),
     & da1(32,1),da2(32,1)
          dimension bigst(1),ib(1),jb(1),rmsn(1)
	dpi=3.1415926536d0
	da=w/2.0d0
	total=imax*jmax
	drsq=4.0d0/dpi
c  DEFAULT OPTIONS AND BOUNDS OF INPUT DATA ARRAY(S)
          iferr=0
        gmin=500000.0
	gmax=-gmin
	do 7 j=1,jmax
	do 6 i=1,imax
	if(iopt3.eq.2) g(i,j)=1.0
	if(iopt5.eq.0) zs(i,j)=zs1
	if(g(i,j).gt.gmax) gmax=g(i,j)
	if(g(i,j).lt.gmin) gmin=g(i,j)
	z=zs(i,j)
	if(z.lt.0.0) iferr=1
	if(iopt3.eq.0) go to 6
	t1=t(i,j)
	if(t1.lt.0.0) iferr=1
	if(iopt1.eq.2.and.(z-t1).lt.0.0) iferr=1
	if(iopt1.eq.3.and.(z-t1/2.0).lt.0.0) iferr = 1
6	continue
7	continue
        if(nmax.le.0.or.nmax.gt.30) nmax=5
      if(iopt3.eq.2) nmax=1
8	if(nstar.le.0.or.nstar.gt.nmax) nstar=nmax+1
c   TEST FOR PARAMETER ERRORS
      if(imax.lt.1.or.imax.gt.30.or.jmax.lt.1.or.jmax.gt.30) iferr=2
	if(w.le.0.0) iferr=3
	if(gmax.gt.0.0.and.gmin.lt.0.0) iferr=4
	if(gmin.lt.0.0.and.rho.ge.0.0) iferr=4
	if(gmin.ge.0.0.and.rho.le.0.0.and.iopt3.ne.2) iferr=4
	eps=abs(eps)
          if(iopt1.ne.1.and.iopt1.ne.2.and.iopt1.ne.3)iferr=1
17	if(iopt2.eq.1) go to 21
	if(iopt2.eq.2) go to 23
	if(iopt2.eq.3) go to 25
          iferr=5
21	c=1.0
	go to 27
23        c=1.60935
	go to 27
25	c=0.304801
27	continue
	dcnst=6.673d0*c*rho*da
	dcnsta=4.0d0*dcnst
          if(iopt3.ne.0.and.iopt3.ne.1.and.iopt3.ne.2) iferr=6
          if(iopt5.ne.0.and.iopt5.ne.1) iferr=1
          if(zs1.lt.0.0) iferr=1
          if(iopt3.eq.2) go to 95
100       if(nstar.eq.1) go to 95
          namax=nstar-1
          if(iferr.ne.0) go to 999
95	continue
c  BEGIN COMPUTATIONS
c  FIRST APPROXIMATION T(I,J),BY BOUGUER SLAB T=G/(2*PI*GAMMA*RHO)
	n=1
	if(iopt3.ne.0) go to 400
	bougr=c*rho*41.9277
	do 35 j=1,jmax
	do 34 i=1,imax
34	t(i,j)=g(i,j)/bougr
35	continue
c  COMPUTATIONS.  NOTE -- (I1,J1)=BODY POINT,(I3,J3)=FIELD POINT,
c   (I2,J2)=FIELD POINT  IN TRANSFORMED FRAME. DA1 = GRAVITY ARRAY,
c   DA2 = STORAGE ARRAY (FOR EXACT FORMULA).
      bigold=1.e38
400	continue
c  CLEAR DGCALC ARRAY
	do 403 j3=1,jmax
	do 402 i3=1,imax
402	dgcalc(i3,j3)=0.0
403	continue
	do 415 j1=1,jmax
	l1=j1+1
	j2max1=jmax-j1+1
	if(j2max1.lt.j1) j2max1=j1
	do 414 i1=1,imax
	l2=i1+1
	i2max=imax-i1+1
	if(i2max.lt.i1) i2max=i1
	j2max=j2max1
	if(i2max.ge.j2max) go to 416
	j2max=i2max
	i2max=j2max1
416	continue
c  COMPUTE DG OF T(I1,J1), IN (I2,J2) FRAME
c  DEFINE PRISM TOP AND BASE
	tt=t(i1,j1)
	zz=zs(i1,j1)
	if(iopt1-2) 418,419,420
418	dz1=zz
	dz2=dz1+tt
	go to 421
419	if(zz-tt) 430,431,431
430	t(i1,j1)=zz
	tt=zz
431	dz1=zz-tt
	dz2=zz
	go to 421
420	th=tt/2.0
	if(zz-th) 432,433,433
432	t(i1,j1)=2.0*zz
	th=zz
433	dz1=zz-th
	dz2=zz+th
421	dz1=dz1/da
	dz2=dz2/da
	if(dz2-dz1) 414,414,434
434	dz1sq=dz1*dz1
	dz2sq=dz2*dz2
c  SELECT EXACT OR APPROXIMATE PRISM FORMULAE
	if(nstar-n) 435,435,450
c  GRAVITY OF N-TH MODEL BY EXACT FORMULA
c  AXIAL POINT (1,1)
435	dr1=dsqrt(1.0+dz1sq)
	dr2=dsqrt(1.0+dz2sq)
	db1=1.0+dsqrt(2.0+dz1sq)
	db2=1.0+dsqrt(2.0+dz2sq)
	da2(1,1)=dcnst*(2.0*dlog((db1*dr2)/(db2*dr1))+
     & dz2*(dasin((dz2sq+db2)/(db2*dr2))-dasin(dz2/dr2)))
	if(dz1)438,438,437
437	da2(1,1)=da2(1,1)+dcnst*dz1*(dasin(dz1/dr1)-
     & dasin((dz1sq+db1)/(db1*dr1)))
438	da1(1,1)=da2(1,1)*4.0
c  LINE (I2,1), I2=2,3...I2MAX
	do 401 i2=2,i2max
	dx=2*i2-1
	dxsq=dx*dx
	dp=1.0+dxsq
	dq1=dsqrt(dp+dz1sq)
	dq2=dsqrt(dp+dz2sq)
	dr1=dsqrt(dxsq+dz1sq)
	dr2=dsqrt(dxsq+dz2sq)
	ds1=dsqrt(1.0+dz1sq)
	ds2=dsqrt(1.0+dz2sq)
	db1=dx+dq1
	db2=dx+dq2
	dc1=1.0+dq1
	dc2=1.0+dq2
	da2(i2,1)=dcnst*(dx*dlog((dc1*dr2)/(dc2*dr1))+dlog((db1*ds2)/
     & (db2*ds1))+dz2*(dasin((dc2+dz2sq)/(dc2*ds2))-dasin(dz2/dr2)))
	if(dz1) 440,440,439
439	da2(i2,1)=da2(i2,1)+dcnst*dz1*(dasin(dz1/dr1)-
     & dasin((dc1+dz1sq)/(dc1*ds1)))
440	da1(i2,1)=2.0*(da2(i2,1)-da2(i2-1,1))
	da1(1,i2)=da1(i2,1)
401	continue
c  DIAGONAL (K,K), K=2,3...J2MAX
	j2=1
422	j2=j2+1
          dy=2*j2-1
          dysq=dy*dy
          dp=2.0d0*dysq
          dq1=dsqrt(dp+dz1sq)
          dq2=dsqrt(dp+dz2sq)
          ds1=dsqrt(dysq+dz1sq)
          ds2=dsqrt(dysq+dz2sq)
          dc1=dy+dq1
          dc2=dy+dq2
          da2(j2,j2)=dcnst*(2.0d0*dy*dlog((dc1*ds2)/(dc2*ds1))+dz2
     &    *(dasin((dz2sq+dy*dc2)/(dc2*ds2))-dasin(dz2/ds2)))
          if(dz1) 442,442,441
441       da2(j2,j2)=da2(j2,j2)+dcnst*dz1*(dasin(dz1/ds1)-dasin((
     &    dz1sq+dy*dc1)/(dc1*ds1)))
442	da1(j2,j2)=da2(j2,j2)-2.0*da2(j2,j2-1)+da2(j2-1,j2-1)
c  WEDGE AND RECTANGLE, (I2,J2), J2=2,3...J2MAX; I2=3,4...I2MAX
	i2=j2
417	i2=i2+1
	dx=2*i2-1
	dxsq=dx*dx
	dp=dxsq+dysq
	dq1=dsqrt(dp+dz1sq)
	dq2=dsqrt(dp+dz2sq)
	dr1=dsqrt(dxsq+dz1sq)
	dr2=dsqrt(dxsq+dz2sq)
	db1=dx+dq1
	db2=dx+dq2
	dc1=dy+dq1
	dc2=dy+dq2
      da2(i2,j2)=dcnst*(dx*dlog((dc1*dr2)/(dc2*dr1))+dy*dlog((db1*ds2)/
     &(db2*ds1))+dz2*(dasin((dz2sq+dy*dc2)/
     & (dc2*ds2))-dasin(dz2/dr2)))
	if(dz1) 444,444,443
443	da2(i2,j2)=da2(i2,j2)+dcnst*dz1*(dasin(dz1/dr1)-dasin((dz1sq+
     & dy*dc1)/(dc1*ds1)))
444	da1(i2,j2)=da2(i2,j2)-da2(i2,j2-1)-da2(i2-1,j2)+da2(i2-1,j2-1)
	da1(j2,i2)=da1(i2,j2)
	if(i2-i2max) 417,436,436
436	if(j2-j2max) 422,407,407
c  GRAVITY OF N-TH MODEL BY APPROXIMATE FORMULAE
450	j2=1
	da1(1,1)=(2.0*dpi*dcnst)*(dsqrt(drsq+dz1sq)-dsqrt(drsq+dz2sq)
     & +dz2-dz1)
	do 451 i2=2,i2max
	dp=4*(i2-1)*(i2-1)
	da1(i2, 1)=dcnsta*((1.0/dsqrt(dp+dz1sq))-(1.0/dsqrt(dp+dz2sq)))
451	da1(1,i2)=da1(i2,1)
452	j2=j2+1
	dy=2*(j2-1)
	dp=2.0*dy*dy
	i2=j2
	da1(i2,j2)=dcnsta*((1.0/dsqrt(dp+dz1sq))-(1.0/dsqrt(dp+dz2sq)))
453	i2=i2+1
	dx=2*(i2-1)
	dp=dx*dx+dy*dy
	da1(i2,j2)=dcnsta*((1.0/dsqrt(dp+dz1sq))-(1.0/dsqrt(dp+dz2sq)))
	da1(j2,i2)=da1(i2,j2)
	if(i2-i2max) 453,454,454
454	if(j2-j2max) 452,407,407
c  TRANSFORM G(I2,J2) INTO G(I3,J3), IN FOUR SECTORS
407	do 410 j3=1,j1
	j2=j1-j3+1
	do 408 i3=1,i1
408	dgcalc(i3,j3)=dgcalc(i3,j3)+da1(i1-i3+1,j2)
	do 409 i3=l2,imax
409	dgcalc(i3,j3)=dgcalc(i3,j3)+da1(i3-i1+1,j2)
410	continue
	do 413 j3=l1,jmax
	j2=j3-j1+1
	do 411 i3=1,i1
411	dgcalc(i3,j3)=dgcalc(i3,j3)+da1(i1-i3+1,j2)
	do 412 i3=l2,imax
412	dgcalc(i3,j3)=dgcalc(i3,j3)+da1(i3-i1+1,j2)
413	continue
414	continue
415	continue
c  COMPLETE DGCALC, TEST N FOR OUTPUT
c  COMPUTE LARGEST AND RMS ERROR
	rms=0.0
	big=0.0
	do 42 j=1,jmax
	do 41 i=1,imax
	da1(i,j)=g(i,j)-dgcalc(i,j)
	adiff=dabs(da1(i,j))
	if(adiff.lt.big) go to 40
	big=adiff
	ibig=i
	jbig=j
40	rms=rms+adiff*adiff
41	continue
42	continue
          rmsn(n)=sqrt(rms/total)
          bigst(n)=big
          ib(n)=ibig
          jb(n)=jbig
          write(6,222) n,bigst(n),ib(n),jb(n),rmsn(n)
222   format(4x,i3,7x,g12.4,4x,i3,',',i2,4x,g12.4)
          if(iopt3.eq.2) go to 50
          if(big.lt.eps) go to 50
          if(n.eq.nmax) go to 50
          if(big.gt.bigold) go to 50
          bigold=big
c  REVISE MODEL AND RECYCLE
	n=n+1
	do 45 j=1,jmax
	do 44 i=1,imax
44	t(i,j)=(t(i,j)*g(i,j))/dgcalc(i,j)
45	continue
	go to 400
c  OUTPUT
50        continue
	do 58 j=1,jmax
	do 57 i=1,imax
	da2(i,j)=zs(i,j)+t(i,j)
	if(iopt1.eq.2) da2(i,j)=zs(i,j)-t(i,j)
	if(iopt1.eq.3) da2(i,j)=zs(i,j)-t(i,j)/2.0d0
57	continue
58	continue
999	return
	end

