anchor.diff:        lwork = int(work(1))
         !print *,'LWORK',lwork
   C Copy observation vector, b to btmp (replaced in dgels)
  -      do 101 i = 1,nobs
  +      do 441 i = 1,nobs
             btmp(i,1) = b(i,1)
  -  101 continue
  +  441 continue
   C Compute parameter array
         call dgels(trans,m,n,nrhs,a,lda,btmp,ldb,work,lwork,info)
  -      do 102 i = 1,npar
  +      do 442 i = 1,npar
             x(i,1) = btmp(i,1)
  -  102 continue
  +  442 continue
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE anneal(soln,gf,obs,nobs,OBSMAX,nflt,FLTMAX,damp,smooth,
  +     1                  smoof,annf)
  +C----
  +C Perform a random search for the optimal solution, utilizing a
  +C simulated annealing algorithm.
  +C----
  +      IMPLICIT none
  +      REAL*8 pi,d2r
  +      PARAMETER (pi=4.0d0*atan(1.0d0),d2r=pi/1.8d2)
  +C Inversion variables
  +      INTEGER OBSMAX,FLTMAX
  +      REAL*8 obs(OBSMAX,6),gf(OBSMAX,FLTMAX,6)
  +      REAL*8 soln(FLTMAX,2),best(FLTMAX,2)
  +      INTEGER nobs,nflt
  +      CHARACTER*80 smoof
  +      REAL*8 smooth,damp
  +      INTEGER nsmoo,smarry(FLTMAX,8)
  +      REAL*8 misfit,length,rough
  +      REAL*8 best0(FLTMAX,2,50),obj0(50)
  +      INTEGER nsoln
  +      CHARACTER*2 nsolnc
  +C Annealing-specific variables
  +      REAL*8 energy(500),obj,temp,objmin
  +      REAL*8 denom,test
  +      INTEGER i,j,jj
  +C Random number variables
  +      INTEGER time(3),timeseed,idum
  +      CHARACTER*8 timec
  +      EXTERNAL ran2
  +      REAL ran2
  +C Control file parameters
  +      CHARACTER*80 annf             ! Name of annealing control file 
  +      INTEGER      inityp           ! Initial solution type: 0=zero; 1=auto; 2=user
  +      REAL*8       slip0,rake0      ! Initial slip, rake values
  +      REAL*8       slip(FLTMAX,3)   ! Slip ranges: slipmn,slipmx,nslip,dslip
  +      REAL*8       rake(FLTMAX,3)   ! Rake ranges: rakemn,rakemx,nrake,drake
  +      REAL*8       temp0            ! Initial temperature
  +      REAL*8       cool             ! Cooling rate
  +      INTEGER      nitmax           ! Maximum iterations to run algorithm
  +C Other variables
  +      REAL*8 dslip(FLTMAX),drake(FLTMAX)
  +
  +      open(unit=51,file='anneal.log',status='unknown')
  +      open(unit=52,file='solutions.log',status='unknown')
  +      write(51,'(A)') 'Starting simulated annealing search'
  +      write(51,*)
  +
  +C----
  +C Read control file
  +C----
  +      write(51,'(A)') 'Setting up control parameters'
  +      call readannf(inityp,slip0,rake0,slip,rake,temp0,cool,nitmax,
  +     1              annf,nflt,FLTMAX,nsoln)
  +      write(51,'("  Initial model (0=zero,1=auto,2=user)")')
  +      write(51,'("      inityp:   ",I10)') inityp
  +      write(51,'("      slip0:    ",F10.2)') slip0
  +      write(51,'("      rake0:    ",F10.2)') rake0
  +      write(51,'("  Parameter ranges:")')
  +      write(51,'("      Fault slip ranges:")')
  +      write(51,'("         flt    slipmn    slipmx nslip")')
  +      do 501 i=1,nflt
  +          write(51,'(8X,I4,2F10.2,F6.0)') i,(slip(i,j),j=1,3)
  +  501 continue
  +      write(51,'("      Fault rake ranges:")')
  +      write(51,'("         flt    rakemn    rakemx nrake")')
  +      do 502 i=1,nflt
  +          write(51,'(8X,I4,2F10.2,F6.0)') i,(rake(i,j),j=1,3)
  +  502 continue
  +      write(51,'("  Cooling schedule:")')
  +      write(51,'("      temp0:    ",F10.4)') temp0
  +      write(51,'("      cool:     ",F10.6)') cool
  +      write(51,'("      nitmax:   ",I10)') nitmax
  +      write(51,'("  Number of solutions stored:")')
  +      write(51,'("      nsoln:    ",I10)') nsoln
  +      write(51,*)
  +
  +C----
  +C Set up array for computing smoothing
  +C----
  +      if (smooth.gt.0.0d0) then
  +          call readsmoof(smarry,nsmoo,smoof,nflt,FLTMAX)
  +          do 503 i = 1,nsmoo
  +              print *,smarry(i,1),smarry(i,2),
  +     1                                     (smarry(i,j),j=1,smarry(i,2))
  +  503     continue
  +      endif
  +
  +C----
  +C Initial solution
  +C----
  +      write(51,'(A)') 'Loading initial solution'
  +      call model0(soln,gf,obs,nobs,OBSMAX,nflt,FLTMAX,inityp,slip0,
  +     1            rake0)
  +      write(51,*)
  +      do 505 j = 1,nsoln
  +          obj0(j) = 1d10
  +  505 continue
  +
  +C----
  +C Initial solution objective function
  +C----
  +      write(51,'(A)') 'Computing initial objective function'
  +      call calcobj(obj,misfit,length,rough,soln,gf,obs,nobs,OBSMAX,nflt,
  +     1             FLTMAX,damp,smooth,smarry,nsmoo)
  +      do 506 i = 1,nflt
  +          objmin = obj
  +          best(i,1) = soln(i,1)
  +          best(i,2) = soln(i,2)
  +          obj0(1) = obj
  +          best0(i,1,1) = soln(i,1)
  +          best0(i,2,1) = soln(i,2)
  +  506 continue
  +      write(51,1001) 
  +      write(51,1002) obj,misfit,damp,length,smooth,rough
  +      write(51,*)
  + 1001 format("   objective =      misfit + (        damp*      length",
  +     1       ") + (      smooth*       rough)")
  + 1002 format(F12.6," =",F12.6," + (",F12.6,"*",F12.6,") + (",F12.6,"*",
  +     1       F12.6,")")
  +
  +C Record solutions in solutions.log
  + 1003 format('Iteration: ',I4,4X,'Temperature: ',F9.6,4X,'Objective=',
  +     1       F12.6)
  + 1004 format("  flt      slip      rake")
  + 1005 format(I5,F10.4,F10.2)
  +      write(52,1003) 0,1e10,obj
  +      write(52,1004)
  +      do 504 i=1,nflt
  +          write(52,1005) i,soln(i,1),soln(i,2)
  +  504 continue
  +      write(52,*)
  +
  +C----
  +C Run annealing search
  +C----
  +      write(51,'(A)') 'Starting annealing'
  +C Initialize random number generator
  +      call itime(time)
  +      write(timec,'(I0.2I0.2I0.2)') time(1),time(2),time(3)
  +      read(timec,*) timeseed
  +      if (timeseed.gt.0) timeseed = -timeseed
  +      idum = timeseed
  +C Store some derived arrays
  +      do 507 i = 1,nflt
  +          dslip(i) = (slip(i,2)-slip(i,1))/slip(i,3)
  +          drake(i) = (rake(i,2)-rake(i,1))/rake(i,3)
  +  507 continue
  +C Iterate and search
  +      temp = temp0
  +      do 520 i = 1,nitmax
  +          do 514 j = 1,nflt
  +C             Gibbs sampler: create variable PDFs; directly sample
  +              denom = 0.0d0
  +              do 508 jj = 1,int(slip(j,3)+0.1)
  +                  soln(j,1) = slip(j,1) + (dble(jj)-0.5d0)*dslip(j)
  +                  call calcobj(obj,misfit,length,rough,soln,gf,obs,nobs,
  +     1                         OBSMAX,nflt,FLTMAX,damp,smooth,smarry,
  +     2                         nsmoo)
  +                  if (obj/temp.gt.50.0d0) then
  +                      energy(jj) = 0.0d0
  +                  else
  +                      energy(jj) = exp(-obj/temp)
  +                  endif
  +                  denom = denom + energy(jj)
  +  508         continue
  +              test = 0.0d0
  +              do 509 jj = 1,int(slip(j,3)+0.1)
  +                  test = test + energy(jj)/denom
  +                  if (ran2(idum).lt.test) then ! Sample PDF
  +                      goto 510
  +                  endif
  +  509         continue
  +  510         soln(j,1) = slip(j,1) + (dble(jj)-0.5d0)*dslip(j)
  +              denom = 0.0d0
  +              do 511 jj = 1,int(rake(j,3)+0.1)
  +                  soln(j,2) = rake(j,1) + (dble(jj)-0.5d0)*drake(j)
  +                  call calcobj(obj,misfit,length,rough,soln,gf,obs,nobs,
  +     1                         OBSMAX,nflt,FLTMAX,damp,smooth,smarry,
  +     2                         nsmoo)
  +                  if (obj/temp.gt.50.0d0) then
  +                      energy(jj) = 0.0d0
  +                  else
  +                      energy(jj) = exp(-obj/temp)
  +                  endif
  +                  denom = denom + energy(jj)
  +  511         continue
  +              test = 0.0d0
  +              do 512 jj = 1,int(rake(j,3)+0.1)
  +                  test = test + energy(jj)/denom
  +                  if (ran2(idum).lt.test) then ! Sample PDF
  +                      goto 513
  +                  endif
  +  512         continue
  +  513         soln(j,2) = rake(j,1) + (dble(jj)-0.5d0)*drake(j)
  +  514     continue
  +C         Compute objective function for new model
  +          call calcobj(obj,misfit,length,rough,soln,gf,obs,nobs,OBSMAX,
  +     1                 nflt,FLTMAX,damp,smooth,smarry,nsmoo)
  +C         Record results in log files
  +          write(51,1003) i,temp,obj
  +          write(52,1003) i,temp,obj
  +          write(52,1004)
  +          do 515 j = 1,nflt
  +              write(52,1005) j,soln(j,1),soln(j,2)
  +  515     continue
  +          write(52,*)
  +C         Save solution in best models if in top NSOLN obj values
  +          do 517 j = 1,nsoln
  +              if (obj.lt.obj0(j)) then
  +                  do 516 jj = 1,nflt
  +                      best0(jj,1,j) = soln(jj,1)
  +                      best0(jj,2,j) = soln(jj,2)
  +  516             continue
  +                  obj0(j) = obj
  +                  goto 518
  +              endif
  +  517     continue
  +  518     continue
  +C         If solution is best, save it
  +          if (obj.lt.objmin) then
  +              do 519 j = 1,nflt
  +                  best(j,1) = soln(j,1)
  +                  best(j,2) = soln(j,2)
  +                  objmin = obj
  +  519         continue
  +          endif
  +C         Cool off
  +          if (temp.gt.1e-4) then
  +              temp = temp*cool
  +          else
  +              temp = temp*0.995d0
  +          endif
  +  520 continue
  +      write(51,*)
  +      write(51,'("Finished annealing")')
  +      close(51)
  +      close(52)
  +
  +C----
  +C Write nsoln best solutions to file best_nsoln.dat
  +C----
  +      write(nsolnc,'(I0.2)') nsoln
  +      open(unit=53,file='best_'//nsolnc//'.dat',status='unknown')
  +      do 521 i = 1,nflt
  +          write(53,'(50F12.6)') (best0(i,1,j),j=1,nsoln)
  +  521 continue
  +      do 522 i = 1,nflt
  +          write(53,'(50F12.6)') (best0(i,2,j),j=1,nsoln)
  +  522 continue
  +      close(53)
  +
  +C----
  +C Return best overall solution
  +C----
  +      do 523 i = 1,nflt
  +          soln(i,1) = best(i,1)*dcos(best(i,2)*d2r)
  +          soln(i,2) = best(i,1)*dsin(best(i,2)*d2r)
  +  523 continue
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE readannf(inityp,slip0,rake0,slip,rake,temp0,cool,
  +     1                    nitmax,annf,nflt,FLTMAX,nsoln)
  +C----
  +C Read annealing control file if it exists
  +C----
  +      IMPLICIT none
  +      INTEGER nflt,FLTMAX
  +      CHARACTER*80 annf
  +      INTEGER inityp,nitmax
  +      INTEGER nsoln
  +      REAL*8 slip0,rake0,slip(FLTMAX,3),rake(FLTMAX,3),temp0,cool
  +      LOGICAL ex
  +      CHARACTER*120 line
  +      CHARACTER*1 dum
  +      INTEGER i,j,ns,nr
  +C If annealing file not specified, use some default values
  +  551 if (annf.eq.'none') then
  +          inityp = 0
  +          slip0 = 0.0d0
  +          rake0 = 0.0d0
  +          do 552 i = 1,nflt
  +              slip(i,1) = 0.0d0
  +              slip(i,2) = 25.0d0
  +              slip(i,3) = 10.0d0
  +              rake(i,1) = -180.0d0
  +              rake(i,2) = 180.0d0
  +              rake(i,3) = 30.0d0
  +  552     continue
  +          temp0 = 1.0d0
  +          cool = 0.900d0
  +          nitmax = 250
  +          nsoln = 1
  +          return
  +      endif
  +C Check if annealing control file exists
  +      inquire(file=annf,exist=ex)
  +      if (.not.ex) then
  +          write(*,*) '!! Warning: no annealing file '//trim(annf)
  +          write(*,*) 'Using default parameters'
  +          annf = 'none'
  +          goto 551
  +      endif
  +C Initialize slip, rake arrays
  +      do 553 i = 1,nflt
  +          slip(i,1) = 0.0d0
  +          slip(i,2) = 25.0d0
  +          slip(i,3) = 25.0d0
  +          rake(i,1) = -180.0d0
  +          rake(i,2) = 180.0d0
  +          rake(i,3) = 36.0d0
  +  553 continue
  +C Read annealing control file
  +      open(unit=59,file=annf,status='old')
  +      read(59,*) line ! auto/user/zero
  +      if (line.eq.'zero') then
  +          inityp = 0
  +      elseif (line.eq.'auto') then
  +          inityp = 1
  +      elseif (line.eq.'user') then
  +          inityp = 2
  +      else
  +          write(*,*) '!! Error: no initialize option '//trim(line)
  +          call usage('!! Use auto, user, or zero')
  +      endif
  +      read(59,*) slip0,rake0
  +      ns = 0
  +      nr = 0
  +  554 read(59,'(A)') line
  +      i = index(line,'s')
  +      if (i.eq.1) then
  +          read(line,*) dum,j,slip(j,1),slip(j,2),slip(j,3)
  +          ns = ns + 1
  +          goto 554
  +      endif
  +      i = index(line,'r')
  +      if (i.eq.1) then
  +          read(line,*) dum,j,rake(j,1),rake(j,2),rake(j,3)
  +          nr = nr + 1
  +          goto 554
  +      endif
  +      if (ns.eq.0.or.nr.eq.0) then
  +          write(*,*) '!! Error: must define slip and rake limits'
  +          write(*,*) '!! s flt slpmn slipmx nslip'
  +          call usage('!! r flt slpmn rakemx nrake')
  +      endif
  +C If only one fault parameter limits is defined, use it for all
  +      if (ns.eq.1.and.nflt.ge.2) then
  +          do 555 i = 2,nflt
  +              slip(i,1) = slip(1,1)
  +              slip(i,2) = slip(1,2)
  +              slip(i,3) = slip(1,3)
  +  555     continue
  +      endif
  +      if (nr.eq.1.and.nflt.ge.2) then
  +          do 556 i = 2,nflt
  +              rake(i,1) = rake(1,1)
  +              rake(i,2) = rake(1,2)
  +              rake(i,3) = rake(1,3)
  +  556     continue
  +      endif
  +      read(line,*) temp0,cool,nitmax
  +      nsoln = 1
  +      read(59,*,end=557) nsoln
  +  557 continue
  +      close(59)
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE readsmoof(smarry,nsmoo,smoof,nflt,FLTMAX)
  +      IMPLICIT none
  +      CHARACTER*80 smoof
  +      INTEGER nflt,FLTMAX,i,j
  +      INTEGER nsmoo,smarry(FLTMAX,8)
  +      call linect(nsmoo,smoof)
  +      if (nsmoo.gt.nflt) then
  +          call usage('!! Error: number of lines in smooth file > nflt')
  +      endif
  +C Read smoothing element linking file
  +      open(unit=63,file=smoof,status='old')
  +      do 561 i = 1,nsmoo
  +          read(63,*) smarry(i,1),smarry(i,2),
  +     1                                     (smarry(i,j),j=1,smarry(i,2))
  +  561 continue
  +      close(63)
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE model0(soln,gf,obs,nobs,OBSMAX,nflt,FLTMAX,inityp,
  +     1                  slip0,rake0)
  +C----
  +C Set up the initial solution
  +C----
  +      IMPLICIT none
  +      REAL*8 pi,d2r
  +      PARAMETER (pi=4.0d0*atan(1.0d0),d2r=pi/1.8d2)
  +      INTEGER nobs,OBSMAX,nflt,FLTMAX,inityp,i,j
  +      REAL*8 obs(OBSMAX,6),gf(OBSMAX,FLTMAX,6),soln(FLTMAX,2)
  +      REAL*8 pre(3),op,pp,slip0,rake0,ss,ds
  +      if (inityp.eq.0) then
  +          slip0 = 0.0d0
  +          rake0  = 0.0d0
  +      elseif (inityp.eq.1) then
  +          ss = dcos(rake0*d2r)
  +          ds = dsin(rake0*d2r)
  +          op = 0.0d0 ! sum of dot product obs.pre
  +          pp = 0.0d0 ! sum of dot product pre.pre
  +          do 572 i = 1,nobs
  +              do 571 j = 1,nflt
  +                  pre(1) = ss*gf(i,j,1) + ds*gf(i,j,4)
  +                  pre(2) = ss*gf(i,j,2) + ds*gf(i,j,5)
  +                  pre(3) = ss*gf(i,j,3) + ds*gf(i,j,6)
  +                  op = op + obs(i,4)*pre(1) + obs(i,5)*pre(2) +
  +     1                                                   obs(i,6)*pre(3)
  +                  pp = pp +   pre(1)*pre(1) +   pre(2)*pre(2) +
  +     1                                                     pre(3)*pre(3)
  +  571         continue
  +  572     continue
  +          slip0 = op/pp
  +      elseif (inityp.eq.2) then
  +          ! slip0 and rake0 already defined
  +      else
  +          call usage('!! Error: no annealing option '//char(inityp))
  +      endif
  +      do 573 i = 1,nflt
  +          soln(i,1) = slip0
  +          soln(i,2) = rake0
  +  573 continue
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE calcobj(obj,misfit,length,rough,soln,gf,obs,nobs,
  +     1                   OBSMAX,nflt,FLTMAX,damp,smooth,smarry,nsmoo)
  +      IMPLICIT none
  +      REAL*8 obj
  +      REAL*8 misfit,length,rough,smooth,damp
  +      INTEGER nobs,OBSMAX,nflt,FLTMAX
  +      REAL*8 obs(OBSMAX,6),gf(OBSMAX,FLTMAX,6),soln(FLTMAX,2)
  +      INTEGER nsmoo,smarry(FLTMAX,8)
  +      call calcmisfit(misfit,soln,gf,obs,nobs,OBSMAX,nflt,FLTMAX)
  +      misfit = dsqrt(misfit/nflt) ! RMS misfit
  +      length = 0.0d0
  +      rough = 0.0d0
  +      if (damp.gt.0.0d0) then
  +          call calclength(length,soln,nflt,FLTMAX)
  +          length = damp*length/nflt ! Average length
  +      endif
  +      if (smooth.gt.0.0d0) then
  +          call calcrough(rough,soln,nflt,FLTMAX,smarry,nsmoo)
  +          rough = smooth*rough/nflt ! Average roughness
  +      endif
  +      obj = misfit + damp*length + smooth*rough
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE calcmisfit(misfit,soln,gf,obs,nobs,OBSMAX,nflt,FLTMAX)
  +      IMPLICIT none
  +      REAL*8 pi,d2r
  +      PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2)
  +      INTEGER nobs,OBSMAX,nflt,FLTMAX
  +      REAL*8 obs(OBSMAX,6),gf(OBSMAX,FLTMAX,6)
  +      REAL*8 soln(FLTMAX,2)
  +      REAL*8 pre(3)
  +      REAL*8 misfit
  +      INTEGER i,j
  +      REAL*8 cr,sr,ss(FLTMAX),ds(FLTMAX)
  +C Store strike-slip and dip-slip components of faults
  +      do 581 i = 1,nflt
  +          cr = dcos(soln(i,2)*d2r)
  +          sr = dsin(soln(i,2)*d2r)
  +          ss(i) = soln(i,1)*cr
  +          ds(i) = soln(i,1)*sr
  +  581 continue
  +C Compute misfit
  +      misfit = 0.0d0
  +      do 583 i = 1,nobs
  +          pre(1) = 0.0d0
  +          pre(2) = 0.0d0
  +          pre(3) = 0.0d0
  +          do 582 j = 1,nflt
  +              pre(1) = pre(1) + gf(i,j,1)*ss(j) + gf(i,j,4)*ds(j)
  +              pre(2) = pre(2) + gf(i,j,2)*ss(j) + gf(i,j,5)*ds(j)
  +              pre(3) = pre(3) + gf(i,j,3)*ss(j) + gf(i,j,6)*ds(j)
  +  582     continue
  +          misfit = misfit + (obs(i,4)-pre(1))*(obs(i,4)-pre(1))
  +     1                    + (obs(i,5)-pre(2))*(obs(i,5)-pre(2))
  +     2                    + (obs(i,6)-pre(3))*(obs(i,6)-pre(3))
  +  583 continue
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE calclength(length,soln,nflt,FLTMAX)
  +      IMPLICIT none
  +      INTEGER nflt,FLTMAX,i
  +      REAL*8 soln(FLTMAX,2),length
  +      length = 0.0d0
  +      do 584 i = 1,nflt
  +          length = length + soln(i,1)
  +  584 continue
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE calcrough(rough,soln,nflt,FLTMAX,smarry,nsmoo)
  +      IMPLICIT none
  +      REAL*8 pi,d2r
  +      PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2)
  +      INTEGER nflt,FLTMAX,nsmoo,smarry(FLTMAX,8),i,j,jj
  +      REAL*8 soln(FLTMAX,2),rough,r(2),ss,ds
  +      r(1) = 0.0d0
  +      r(2) = 0.0d0
  +      do 586 i = 1,nsmoo
  +          ss = soln(i,1)*dcos(soln(i,2)*d2r)
  +          ds = soln(i,1)*dsin(soln(i,2)*d2r)
  +          r(1) = r(1) + smarry(i,2)*ss
  +          r(2) = r(2) + smarry(i,2)*ds
  +          do 585 j = 1,smarry(i,2)
  +              jj = smarry(i,j+2)
  +              ss = soln(jj,1)*dcos(soln(jj,2)*d2r)
  +              ds = soln(jj,1)*dsin(soln(jj,2)*d2r)
  +              r(1) = r(1) - ss
  +              r(2) = r(2) - ds
  +  585     continue
  +  586 continue
  +      rough = r(1) + r(2)
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE chksol(x,nflt)
  +      IMPLICIT none
  +      INTEGER OBSMAX
  +      PARAMETER (OBSMAX=2500)
  +      REAL*8 x(OBSMAX,1)
  +      INTEGER nflt,i
  +      do 611 i = 1,nflt
  +          if (x(i,1).gt.100.0d0) then
  +              write(*,*) '!! Error: found slip larger than 100 m'
  +              write(*,*) '!! Check inputs and units'
  +              write(*,*) '!! Try adding damping/smoothing'
  +              stop
  +          endif
  +  611 continue
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE output(soln,nflt,FLTMAX,ofile,p)
  +      IMPLICIT none
  +      INTEGER FLTMAX
  +      REAL*8 soln(FLTMAX,2)
  +      INTEGER nflt,p,i
  +      CHARACTER*80 ofile
  +      if (p.eq.0) then
  +          open(unit=23,file=ofile,status='unknown')
  +          do 601 i = 1,nflt
  +              write(23,1003) soln(i,1),soln(i,2)
  +  601     continue
  +          close(23)
  +      else
  +          do 602 i = 1,nflt
  +              write(*,1003) soln(i,1),soln(i,2)
  +  602     continue
  +      endif
  + 1003 format(2F14.6)
         RETURN
         END
   
  
positive.code:               read(tag,'(BN,F10.0)') smooth
                i = i + 1
                call getarg(i,smoof)
            elseif (tag(1:6).eq.'-debug') then
                dbg = 1
            elseif (tag(1:5).eq.'-damp') then
                i = i + 1
                call getarg(i,tag)
                read(tag,'(BN,F10.0)') damp
            elseif (tag(1:2).eq.'-2') then
                dimen = 2
            elseif (tag(1:2).eq.'-3') then
                dimen = 3
            elseif (tag(1:2).eq.'-o') then
                i = i + 1
                call getarg(i,ofile)
  
negative.code: 
  C----------------------------------------------------------------------C
  
        SUBROUTINE gcmdln(fltf,obsf,smoof,smooth,damp,ofile,p,geo,haff,
       1                  invert,annf)
        IMPLICIT NONE
        CHARACTER*80 fltf,obsf,tag,smoof,ofile,haff,annf
        INTEGER i,narg,indx
        REAL*8 smooth,damp
        INTEGER geo,p,invert,vrb
        COMMON /VERBOSE/ vrb
        fltf = 'none'
        obsf = 'none'
        smoof = 'none'
        ofile = 'none'
        haff  = 'none'
        annf  = 'none'
        smooth = -1.0d0
        damp = -1.0d0
        geo = 0
        p = 0
        invert = 0
        vrb = 0
        indx = 0
        narg = iargc()
        if (narg.eq.0) call usage('!! Error: no command line arguments')
        i = 0
    101 i = i + 1
        if (i.gt.narg) goto 102
            call getarg(i,tag)
            if (tag(1:4).eq.'-flt') then
                i = i + 1
  