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
      implicit double precision (d)
      dimension g(50,50),dgcalc(50,50),t(50,50),zs(50,50),
     &da1(52,52),da2(52,52),id(14),pgm(2),bigst(50),ib(50),
     &jb(50),rmsn(50)
c      character*50 cmd_file_name
        character*50 cfn
c      character*50 ref_surface_array_file_name
        character*50 rsafn
c      character*50 gravity_data_file_name
        character*50 gdfn
c      character*50 model_array_file_name
        character*50 mafn
c      character*50 output_model_file_name
        character*50 omfn
c      character*50 derived_surface_file_name
        character*50 dsfn
c      character*50 gravity_residual_file_name
        character*50 grfn
c      character*50 calculated_gravity_file_name
        character*50 cgfn
      character*1 answer
      equivalence (imax,ncol),(jmax,nrow),(iopt3,mode)
      common /parms/ w,zs1,rho,nmax,nstar,eps,mode,iopt1,iopt5       
      imax=0
      jmax=0
      rho=0.0
      nstar=0
      nmax=0
      eps=0.0
      mode=0
      iopt2=1
      iopt4=0
      iopt5=0
      iopt6=0
      zs1=0
c  read command and data files
	print*
      write(6,1)
1     format(1x,'Enter command file name:',/,' *',$)
      read(5,2) cfn
2     format(a50)
       open(10,file=cfn,status='old')
c      read (10,parms,end=99)
        call namemc(10)
      close(10)
      write(6,40)zs1,rho,mode,iopt1
40    format(1x,'Following are the basic input parameters:',/,1x,'zs1 ='
     & ,e12.6,10x,'rho =',e12.6,10x,/,1x,'mode =',i2,10x,
     & 'iopt1 =',i2/)
      if(nmax.eq.0)nmax=5
      write(6,41) nmax,nstar,eps,iopt5
41    format(1x,'Following are input parameters which may have been
     & assigned by default:',/,1x,'nmax =',i2,10x,'nstar =',i2,
     & 10x,'eps =',e12.6,10x,'iopt5 =',i2,/,1x,'Want to continue ?',
     & 1x,$)
      read(5,42) answer
42    format(a1)
      if(answer.eq.'y'.or.answer.eq.'Y') go to 44
55    write(6,43)
43    format(1x,'Return to command level to revise parameters.')
      go to 999
44    continue
      if(mode.eq.2) go to 56
      if(mode.eq.1.or.mode.eq.0) go to 50
      iferr=6
      go to 90
50    write(6,51)
51    format(/1x,'Enter gravity data file name:',/' *',$)
      read(5,2) gdfn
      open(10,file=gdfn,status='old',
     & form='unformatted')
      read(10)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
	if(ncol.gt.50.or.nrow.gt.50) then
	write(6,52)
52    format(' Error: ncol or nrow > 50.')
	go to 999
	endif
	if(rdx.ne.rdy) then
      print *,'Error: dx must equal dy.'
	go to 999
	endif
	if(nz.ne.1) then
      print *,'Error: nz must equal 1.'
	go to 999
	endif
	w=rdx
      do 61 j=1,nrow
      read(10) rdum,(g(i,j),i=1,ncol)
	do 662 k=1,ncol
	if(g(k,j).ge.1.e+30) then
	print*,'Can''t handle dvals.'
	go to 999
	endif
662	continue
61	continue
      close(10)
56    if(mode.ne.1.and.mode.ne.2) go to 57
53    write(6,54)
54    format(1x,'Enter model array file name:',/' *',$)
      read(5,2) mafn
       open(10,file=mafn,status='old',
     & form='unformatted')
      read(10)id,pgm,mcol,mrow,mz,xx0,xdx,yy0,ydy
	if(mode.eq.2)then
	  rdx=xdx
	  rdy=ydy
	  ncol=mcol
	  nrow=mrow
	  nz=mz
	  x0=xx0
	  y0=yy0
	  if(ncol.gt.50.or.nrow.gt.50) then
	    write(6,52)
	    go to 999
	  endif
	  if(rdx.ne.rdy) then
	    print*,'Error: dx must equal dy.'
	    go to 999
	  endif
	  if(nz.ne.1) then
	    print*,'Error: nz must equal 1.'
	    go to 999
	  endif
	  w=rdx
	  go to 661
	endif
	if(mcol.ne.ncol.or.mrow.ne.nrow.or.mz.ne.nz) go to 660
	if(xdx.ne.rdx.or.ydy.ne.rdy) go to 660
	slop=.00001
	if(xx0.gt.x0+slop.or.xx0.lt.x0-slop) go to 660
	if(yy0.gt.y0+slop.or.yy0.lt.y0-slop) go to 660
	go to 661
660	print*,'Grid specifications don''t match'
	go to 999
661   do 62 j=1,nrow
62    read(10) rdum,(t(i,j),i=1,ncol)
      close(10)
57    if(iopt5.eq.0) go to 60
      write(6,58)
58    format(1x,'Enter reference surface array file name:',/' *',$)
      read(5,2) rsafn
      open(10,file=rsafn,status='old',
     & form='unformatted')
      read(10)id,pgm,mcol,mrow,mz,xx0,xdx,yy0,ydy
	if(mcol.ne.ncol.or.mrow.ne.nrow.or.mz.ne.nz) go to 664
	if(xdx.ne.rdx.or.ydy.ne.rdy) go to 664
	slop=.00001
	if(xx0.gt.x0+slop.or.xx0.lt.x0-slop) go to 664
	if(yy0.gt.y0+slop.or.yy0.lt.y0-slop) go to 664
	go to 663
664	print*,'Grid specifications don''t match'
	go to 999
663    do 63 j=1,nrow
63    read(10) rdum,(zs(i,j),i=1,ncol)
      close(10)
60    continue
      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,zs1,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
30    write(6,31)
31    format(1x,'Computation completed.')
      if(iopt3.eq.2)go to 32
c      write(6,21) nmax
21    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)
23    format(1x,'At iteration ',i3, ' the largest error is ',
     & e12.4,' at point '
     & ,i3,',',i2,' ;'/1x,'    rms error =',e12.4,'.')
32    write(6,100)
100   format(1x,'Want to see some profiles? [n]')
      read(5,42) answer
      if(answer.ne.'y'.and.answer.ne.'Y')go to 120
101   write(6,102)
102   format(1x,'Row or column?')
      read(5,42)answer
      if(answer.eq.'r'.or.answer.eq.'R')go to 103
      if(answer.eq.'c'.or.answer.eq.'C')go to 108
      go to 101
103   write(6,104)
104   format(1x,'Enter row number.')
      read(5,*) jrow
      write(6,140)
140   format(5x,'Gravity   Calc Grav   Thickness')
      do 106 i=1,ncol
106   write(6,107)i,g(i,jrow),dgcalc(i,jrow),t(i,jrow)
107   format(1x,i2,3f10.3)
      go to 111
108   write(6,109)
109   format(1x,'Enter column number.')
      read(5,*)icol
      write(6,140)
      do 110 j=1,nrow
110   write(6,107)j,g(icol,j),dgcalc(icol,j),t(icol,j)
111   write(6,112)
112   format(1x,'More profiles?')
      read(5,42) answer
      if(answer.eq.'y'.or.answer.eq.'Y')go to 101
120   continue
c  output standard files
caution: note that in this section array t is redefined.
      if(iopt3.eq.2)go to 130
121   write(6,123)
123   format(1x,'Enter output thickness file name (or hit return)')
      read(5,2)omfn
	if(omfn.eq.' ') go to 126
      write(6,127)
127   format(1x,'Enter title (record id of standard file).')
      read(5,128)(id(i),i=1,14)
128   format(14a4)
      open(11,file=omfn,status='unknown',
     & form='unformatted')
      write(11)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
      do 124 j=1,nrow
      y=j
124   write(11)y,(t(i,j),i=1,ncol)
      close(11)
      write(6,125)omfn
125   format(1x,'Segment ',a50,/' contains derived thickness',
     & ' in standard fmt.')
126   continue
141   write(6,143)
143   format(1x,'Enter derived surface file name (or hit return).')
      read(5,2)dsfn
	if(dsfn.eq.' ') go to 150
      write(6,127)
      read(5,128)(id(i),i=1,14)
      open(11,file=dsfn,status='unknown',
     & form='unformatted')
      write(11)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
      do 144 j=1,nrow
      do 146 i=1,ncol
146   t(i,j)=sngl(da2(i,j))
      y=j
144   write(11)y,(t(k,j),k=1,ncol)
      close(11)
      write(6,145)dsfn
145   format(1x,'Segment ',a50,/' contains derived surface',
     & ' standard file.')
150   continue
130   write(6,132)
132   format(1x,'Enter calculated gravity file name (or hit return).')
      read(5,2)cgfn
	if(cgfn.eq.' ') go to 135
      write(6,127)
      read(5,128)(id(i),i=1,14)
      open(11,file=cgfn,status='unknown',
     & form='unformatted')
      write(11)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
      do 136 j=1,nrow
      do 133 i=1,ncol
133   t(i,j)=sngl(dgcalc(i,j))
      y=j
136   write(11)y,(t(k,j),k=1,ncol)
      close(11)
      write(6,134)cgfn
134   format(1x,'Segment ',a50,/' contains calculated gravity array in
     & standard file format.')
135   continue
151   write(6,153)
153   format(1x,'Enter gravity residual file name (or hit return).')
      read(5,2)grfn
	if(grfn.eq.' ') go to 160
      write(6,127)
      read(5,128)(id(i),i=1,14)
      open(11,file=grfn,status='unknown',
     & form='unformatted')
      write(11)id,pgm,ncol,nrow,nz,x0,rdx,y0,rdy
      do 155 j=1,nrow
      do 154 i=1,ncol
154   t(i,j)=g(i,j)-sngl(dgcalc(i,j))
      y=j
155   write(11)y,(t(k,j),k=1,ncol)
      close(11)
      write(6,156)grfn
156   format(1x,'Segment ',a50,/' contains residual',
     & ' gravity standard file.')
160   continue
      go to 999
99    continue
      write(6,900)
900   format(' Abnormal termination due to problem',
     & ' in command file.')
999   stop
      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(50,1),dgcalc(50,1),t(50,1),zs(50,1),
     & da1(52,1),da2(52,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=1000000.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.50) 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.50.or.jmax.lt.1.or.jmax.gt.50) 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).
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
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
      subroutine namemc(icmd)
c
c     namelist simulator for pc's with no namelist statement
c     this subroutine is program independent but it
c     calls subroutine check which is program dependent
c     The program does not handle the repeat indicator *
c     nor does it handle subscripted array variables or
c     substrings of character variables.
c     It also does not handle complex or logical variables.
c     nn = number of characters in value name(tvar)
c     chv = character variable indicator(logical)
c     inum = an array index number
c     nvar = a number returned by subroutine check to
c     determine whether a variable is an array
c     numa = a number returned by subroutine check
c     which is compared against nvar
c     r.godson,usgs,denver,co., 11/87
c
      parameter(ivar=80)
      character var(ivar),comma,apos,blank,equal,amp,dollar
      character*6 pvar
      character*56 tvar
      logical chv,vset
      data comma/','/,apos/#27/,blank/' '/,equal/'='/
      data amp/'&'/,dollar/'$'/,pvar/' '/,tvar/' '/
c
c     get namelist start name(e.g.&parms)
c
    5 read(icmd,1000,end=910) var
      do 7 i=1,ivar
      if(var(i).ne.amp.and.var(i).ne.dollar) go to 7
      do 6 m=i+1,ivar
      if(var(m).eq.blank) go to 15
    6 continue
    7 continue
      go to 5
c
c     start processing variables
c
   10 read(icmd,1000,end=900) var
 1000 format(80a1)
      m=1
   15 do 20 i=m,ivar
      if(var(i).eq.amp.or.var(i).eq.dollar) go to 900
      if(var(i).ne.blank.and.var(i).ne.comma) go to 30
   20 continue
      go to 10
   30 mm=0
c
c     check for continuation of array values on succeding lines
c
      item=ichar(var(i))
      if(item.eq.43.or.item.eq.45.or.item
     & .eq.46.or.(item.gt.47.and.item.lt.58).or.var(i).eq.apos) then
      k=i
      go to 95
      endif
c
c     get program variable name
c
   35 inum=1
      pvar=blank
      vset=.false.
      do 40 j=i,ivar
      if(var(j).eq.blank.or.var(j).eq.equal) go to 50
      mm=mm+1
      pvar(mm:mm)=var(j)
   40 continue
c
c     get variable value
c
     

   50 do 80 k=j+1,ivar
      if(var(k).eq.blank.or.var(k).eq.equal) go to 80
      if(var(k).eq.apos) then
c
c     character variable
      chv=.true.
      vset=.true.
      nn=0
      do 60 l=k+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
   60 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      else
c
c     not a character variable
      chv=.false.
      nn=0
      do 70 l=k,ivar
      if(var(l).eq.blank.or.var(l).eq.comma
     & .or.var(l).eq.amp.or.var(l).eq.dollar) go to 90
c
c     check for non numeric character
c
      item=ichar(var(l))
      if((item.lt.48.and.(item.ne.46
     & .and.item.ne.45.and.item.ne.43)) .or. (item .gt. 57 .and.
     & (item .ne. 69 .and. item .ne. 101 .and. item.ne.68
     &  .and. item.ne.100))) then
      write(*,*) 'missing first apostropy or a non numeric',
     & ' character in namelist variable ',pvar
      stop
      endif
      nn=nn+1
      tvar(nn:nn)=var(l)
   70 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no delimiter'
      stop
      endif
   80 continue
      write(*,*) 'error in namelist variable ',pvar,' no  value'
      stop
c
c     call program dependent subroutine to assign values
c
   90 call check(pvar,tvar,nn,chv,nvar,numa,inum)
      tvar=blank
      k=l+1
c
c     check to see if variable is an array
c     the variable numa used is program dependent
c
      if(nvar.lt.numa) go to 110
   95 if(.not.chv) then
c
c     array variable
c
      nn=0
      do 100 l=k,ivar
      if(var(l).eq.blank.and.nn.eq.0) go to 100
      if(var(l-1).eq.blank.and.var(l).eq.comma.and.nn.eq.0) go to 100
      if((var(l).eq.amp.or.var(l).eq.dollar).and.nn.eq.0) go to 900
c
c     check for consecutive commas
c
      if(var(l).eq.comma.and.vset) then
      vset=.false.
      inum=inum+1
      go to 90
      else if(var(l).eq.comma) then
      inum=inum + 1
      go to 100
      endif
c
      if(var(l).eq.blank
     & .or.var(l).eq.amp.or.var(l).eq.dollar) then
      inum=inum+1
      go to 90
      endif
      if(ichar(var(l)).gt.57.and.nn.eq.0) go to 120
      nn=nn+1
      tvar(nn:nn)=var(l)
      vset=.true.
  100 continue
      else
c
c     character array variable
c
      do 105 j=k,ivar
      if(var(j).eq.blank) go to 105
      if(var(j).eq.amp.or.var(j).eq.dollar) go to 900
      if(ichar(var(j)).gt.57) then
      l=j
      go to 120
      endif
c
c     check for consecutive commas
c
      if(var(j).eq.comma.and.vset) then
      vset=.false.
      go to 105
      else if(var(j).eq.comma) then
      inum=inum + 1
      go to 105
      endif
      if(var(j).eq.apos) then
      nn=0
      inum=inum+1
      vset=.true.
      do 103 l=j+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
  103 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      endif
  105 continue
c
c     end of array values
c
      endif
  110 m=l+1
      if(var(l).eq.amp.or.var(l).eq.dollar) go to 900
      go to 130
  120 m=l
  130 if(m.lt.81) then
      go to 15
      else
      go to 10
      endif
  900 return
  910 write(*,*) ' error in namelist;no beginning & or $'
      stop
      end
        subroutine check(pvar,tvar,nn,chv,nvar,numa,inum)
c
c       assigns values to proper variables
c       variables are passed to program gi3 through common blocks
c       numr=position in the array var where real varialbes start
c       numa=position in the array var where arrays start
c       nnvar=number of variables in program gi3
c
        parameter (nnvar=9,numr=6)
        character*6 pvar,var(nnvar)
        character*56 tvar, kvar, cfmt
        logical chv  
        common/parms/w,zs1,rho,nmax,nstar,eps,mode,iopt1,iopt5
        data var/'nmax','nstar','mode','iopt1','iopt5',
     1  'w','zs1','rho','eps'/
        inum=1
        numa= 10
        do 190 i=1,nnvar
        if (pvar .ne. var(i)) goto 190
        if (.not. chv) then
c
c       non character value
c       right justify the number in variable kvar
c
        m = 57 - nn
        im = m - 1
        kvar(m:56) = tvar(1:nn)
        if (i .lt. numr) then
c
c       integer value
c
        write(cfmt, 50) im, nn
50      format ('(',i2,'x,i',i2,')')
        read(kvar, cfmt) jvar
        else
c
c       real value
c
        write (cfmt, 60) im, nn
60      format ('(',i2,'x,g',i2, '.0)')
        read (kvar, cfmt) xvar
        endif
        endif
        goto (101,102,103,104,106,107,108,109,110),i
101     nmax = jvar
        goto 200
102     nstar = jvar
        goto 200
103     mode = jvar
        goto 200
104     iopt1 = jvar
        goto 200
106     iopt5 = jvar
        goto 200
107     w = xvar
        goto 200
108     zs1 = xvar
        goto 200
109     rho = xvar
        goto 200
110     eps = xvar
        goto 200
190     continue
        write(*,*)' error in namelist - ',pvar,' variable not included'
        stop
200     nvar = i
        return
        end


