program dplus
c  ================  D+ models with 2-norm fit =========================
c  Solves E-M induction inverse problem for the class of conductivity
c  models D+ based on theory of R L Parker JGR 1980, v85, p4241.
c  See also Parker & Whaler JGR 1981, v86, p9574.
      implicit double precision (a-h,o-z)
c$$$$ calls ascan getint getvec getdat setef norm2 misfit lsoadal model
c$$$$ sample 
c  The following limitations are implicit in the array sizes reserved in
c  the following common statements:
c  Up to  maxf frequencies;  up to maxunk samples in lambda.
      parameter (maxf=50, maxunk=500, maxe=22000)
      common //e(maxe),f(2*maxf+2),x(maxunk),iwrk(maxunk)
      common /soln/ nonzer(maxunk),al(2,maxunk),nal
      common /datbas/ omega(maxf),cc(2,maxf),err(2,maxf),nfreq
      common /tactic/ omo,els(maxunk),kind
      common /cases/ plist(25), chisq(25), kase
      common /option/ ifreq,itype,kount,iverb
      common /slabs/ dsq,bigh0,thick
      common /prolix/ iprint
      common /mks/ emuo
      common /bottom/ czero
c
c
      iverb=0
      dsq=0.0
      emuo=4.0d-7 * 3.1415926535d0
c
c  Get command set and set to be worked on
 1000 call ascan
      call getint('print', iprint, ignor)
      call getdat
c
c  In normal operation (iverb=0)  skip multi-case system
      ncase=1
      plist(1)=0.0
      call getvec('base', plist(2), 25, nfound)
      ncase=ncase + max(0,nfound)
      iverb=ncase - 1
      if (ncase .eq. 1) write(*,'(//a/)')
     $' Construction of the best fitting model'
      if (ncase .gt. 1) write(*,'(//a)')
     $' Analysis of penetration depth'
c
c  Run through the list of parameters (czero)
      do 1900 kase=1, ncase
        kount=1 + kount
        czero=plist(kase)
        if (iverb .gt. 0) write(*,'(a,1p,g15.5)')
     $   ' Termination depth (m)',czero
         m=2*nfreq
        if (czero .ne. 0.0) m=m + 2
        me=m
        oldo=1.0d10
c  Finds solution, refining the lambda sampling with subroutine sample
c  the 2nd cycle of the inner loop only executes once, giving final
c  refinement, free of doublets in lambda
        do 1600 lcycle=1,2
          do 1500 kcycle=1,3
            call sample(n, els, kcycle-lcycle+1)
            if (me*n+n+me .gt. maxe) then
              write(*,*)
     $ '>>>> Too many lambdas for array space - refinement skipped'
              oldo=obj
              goto 1500
            endif
            call setef(me, m, n, leq, e, f)
            call norm2(me, m-leq, n, e, f, x, e(me*n+1),iwrk, obj)
            if (obj.gt. oldo) then
              write(*,'(/a,1p,g15.6/a)') '>>>> Objective has increased:'
     $        ,obj,'     Solution before this one accepted.'
              obj=oldo
              goto 1700
            endif
            call loadal(n, x, els, npos)
            write(*,'(/a,i5/a,i5/a,1p,g15.6)')
     $      ' Number of samples in lambda',n,
     $      ' Number of positive elements',npos,
     $      ' Quadratic objective  ',obj
            oldo=obj
            if (obj.le. 0.0001) goto 1700
            if (lcycle.eq.2) goto 1600
*           if (kcycle.gt.1 .and. obj.gt. 0.99*oldo) goto 1600
 1500     continue
 1600   continue
c
c  Display performance of the model
 1700   chisq(kase)=obj
        call misfit
c  Fit a delta-function model
        call model
 1900 continue
c
c  Summarize the results
      if (ncase.gt.1) write(*,200) (plist(k),chisq(k), k=1,ncase)
 200  format(/// 23x, 'Base depth',6x,'misfit'/(19x,1p,2g15.5))
c
      write(*,'(//a)') 'Enter new commands, or enter quit'
      goto 1000
      end
c_______________________________________________________________________
      subroutine norm2(me, m, n, e, f, x, dual, iw, obj)
      implicit double precision (a-h,o-z)
c$$$$ calls nnls
c  The array e is m x n, the 1st dimension in calling routine is me.
c  See comments in nnls
c
      dimension e(me,*),f(*),x(*),dual(*),iw(*)
c
      call nnls(e, me, m, n, f, x, obj, dual, dual(n+1), iw, mode)
c
      obj=obj**2
      if (mode.ne.1) write(*,*)'>>>> Error',mode,' in nnls'
      return
      end 
c_______________________________________________________________________
      subroutine setef(me, m, n, leq, e, f)
      implicit double precision (a-h,o-z)
c$$$$ calls rotate
c  Sets up constraint (sub) matrix for Cauer representation of the e-m
c  induction inverse problem.  D+ conductivities  are constructed.
c  Uses admittance data in /datbas/ and
c  determines sampling tactics from /tactic/
c
      parameter (maxf=50, maxunk=500)
      dimension e(me,n),f(me)
      common /prolix/ iprint
      common /option/ ifreq,itype,kount,iverb
      common /datbas/ omega(maxf),cc(2,maxf),err(2,maxf),nfreq
      common /tactic/ omo,els(maxunk),kind
      common /bottom/ czero
c
      lf=0
      leq=0
c  Fill  e  array and  f  vector
      do 1500 ii=1,m,2
        if(ii.gt. 1 .or. czero.eq. 0.0) then
          lf=1 + lf
          wt1=1.0/err(1,lf)
          wt2=1.0/err(2,lf)
          f(ii)  =wt1*cc(1,lf)
          f(ii+1)=wt2*cc(2,lf)
        else
          wt1=10000.0/czero
          wt2=wt1
          f(ii)  =10000.0
          f(ii+1)=0.0
        endif
c
        e(ii,1)=wt1
        e(ii+1,1)=0.0
c
        if (lf.eq.0) then
          omlf=0.0
        else
          omlf=omega(lf)
        endif
        do 1400 jj=2,n
          alam=els(jj)
c
          e(ii  ,jj)= wt1*alam/(alam**2 + omlf**2)
          e(ii+1,jj)=-wt2*omlf/(alam**2 + omlf**2)
 1400   continue
 1500 continue
c  Rotate matrices if these are apparent resistivity/phase data
      if (itype .eq. 2) goto 2000
      i1=0
      if (czero.ne. 0.0) i1=2
      call rotate(me, m-i1, n, e(i1+1, 1))
      call rotate(me, m-i1, 1, f(i1+1))
c  List the arrays if print priority is high enough
 2000 if (iprint.le.2) return
      write(*,'(/a/(1p,8g13.4))')
     $'Scaled data vector', (f(j),j=1,m)
      write(*,'(/a)') ' Scaled constraint array'
      do 2100 i=1,m
        write(*,'(a,i4/(1p,8g13.4))') 'row',i,(e(i,j),j=1,n)
 2100 continue
      return
      end
c_______________________________________________________________________
      subroutine forwrd(omega, c)
      implicit double precision (a-h,o-z)
c$$$$ calls nothing
c  Finds real-imag parts of admittance  c  at freq  omega  given
c  array of residues-poles  al , for solutions in D+
c
      parameter (maxunk=500)
      common /soln/ nonzer(maxunk),al(2,maxunk),nal
      dimension c(2)
c
      c(1)=al(1,1)
      c(2)=0.0
      do 1100 j=2,nal
        d=al(1,j)/(al(2,j)**2 + omega**2)
        c(1)=c(1) + d*al(2,j)
        c(2)=c(2) - d*omega
 1100 continue
      return
      end
c_______________________________________________________________________
      subroutine model
      implicit double precision (a-h,o-z)
c$$$$ calls qd plot
c  Derives and lists D+ models given partial-fraction
c  representation of admittance in /soln/
      parameter (maxunk=500)
      common /soln/ nonzer(maxunk),al(2,maxunk),nal
      common /prolix/ iprint
      common /option/ ifreq,itype,kount,iverb
      common /mks/ emuo
      dimension h(2*maxunk),cf(2*maxunk)
      dimension z(maxunk),tau(maxunk),dz(maxunk)
c
      nal2=2*nal - 2
c  Usings Rutishauser's qd algorithm, find equivalent continued fraction
      do 1100 j=1, nal2, 2
        jj=(j+3)/2
        h(j)  =al(1,jj)
        h(j+1)=al(2,jj)
 1100 continue
      call qd(nal-1, h, cf(2))
c  Note indexing runs from surface to bottom, contrary to convention
c  of Parker (JGR 1980, vol 85, p4421).  Also applies to model 
c  description.
      cf(1)=al(1,1)
      if (al(2,2).eq. 0.0) cf(nal2+1)=0.0
      if (iprint .ge. 1) write(*,120) cf(1),(cf(j+1),j=1,nal2)
 120  format(/' Continued fraction coeffs'/(1p,8g10.3))
c  Find and print the actual model from continued fraction coeffs
      dz(1)=cf(1)
      z(1)=dz(1)
c  Sort out length units.  If mu0 is standard scale by 1000 for km
      if ( abs(emuo*795775.0 - 1.0) .gt. .0001) then
        skilo=1.0
        write(*,*) 'Nonstandard units are in use'
      else
        skilo=0.001
        write(*,'(/a)') ' Units are kilometers and siemens'
      endif
      write(*,146)
 146  format( ' level',2x,'depth',6x,'conductance',3x,'separation')
      do 1500 j=2,nal2,2
        l=j/2
        if (j.gt.2) tau(l)=1.0/(emuo*dz(l)*cf(j))
        if (j.eq.2) tau(l)=1.0/(emuo*cf(j))
        write(*,150) l,skilo*z(l),tau(l),skilo*dz(l)
 150    format(i4,1p,3g14.5)
        if (cf(j+1).eq. 0.0) goto 1500
        dz(l+1)=1.0/(emuo*tau(l)*cf(j+1))
        z(l+1)=z(l) + dz(l+1)
 1500 continue
c  Print out the depth to bottom if finite
c  If perfect conductor, z(nal)=bot and tau(nal)=-1
c  If insulating base,   z(nal)=0   and tau(nal)=0
      if (al(2,2).gt.0.0)  then
        bot=z(l+1)
        tau(l+1)=-1.0d0
        write(*,155) skilo*z(l+1),skilo*dz(l+1)
 155    format(4x,1p,g14.5,5x,'----',5x,g14.5)
        write(*,*) 'Model terminates with a perfect conductor'
      else
        bot=z(l)
        z(l+1)=0.0d0
        tau(l+1)=0.0d0
        write(*,*) 'Model terminates with an insulator'
      endif
c
      call plot(nal, z,tau,bot, skilo)
      return
      end
c_______________________________________________________________________
      subroutine loadal(n, x, els, npos)
      implicit double precision (a-h,o-z)
c$$$$ calls nothing
c  Given QP solution  x  and associated  lambdas  els, forms compact
c  array  al  containing non-zero  x's  and their  lambdas.
c  Also lists the solution in compact form.
      parameter (maxunk=500)
      common /soln/ nonzer(maxunk),al(2,maxunk),nal
      common /prolix/ iprint
      dimension x(n),els(n)
c
      if (iprint .ge. 2) write(*,130) x(1)
 130  format(/'a(0) = ',1p,g13.4/2x,'i',3x,'k',4x,'lambda',8x,'a(i)')
c
      nal=1
      al(1,1)=x(1)
      al(2,1)=0.0
      do 1500 j=2,n
        if (x(j).le.0.0) goto 1500
        nal=1+nal
        al(1,nal)=x(j)
        al(2,nal)=els(j)
        nonzer(nal)=j
        if (iprint .ge. 2) write(*,150) nal-1,j,els(j),x(j)
 1500 continue
 150  format(i3,i4,1p,2g13.4)
      npos = nal
      if (x(1).le. 0.0) npos=nal-1
      return
      end
c_______________________________________________________________________
      subroutine getdat
      implicit double precision (a-h,o-z)
c$$$$ calls getchr getone getvec
c  Loads /datbas/ with response data from disk file
c
      parameter (pi=3.14159265359d00,  maxf=50)
      common /prolix/ iprint
      common /mks/ emuo
      common /option/ ifreq,itype,kount,iverb
      common /datbas/ omega(maxf),cc(2,maxf),err(2,maxf),nfreq
      common /fildat/ zz(maxf,5)
      character*80 name, row, kind*60
      dimension iover(3),scales(5)
      character*6 fmode(3)
      data fmode /' freq ','omega','period'/
      data name/' '/, kind/'ad'/
      data iover /1,1,-1/, scales/5 * 1.0 /
c
      radian=pi/180.0
c  Get file name
      call getchr('file', name, ignor)
      if (name(1:1) .eq. ' ') then
        write(*,*) '>>>> Data file name must be supplied'
        stop
      elseif (name(1:1).eq. '*') then
        in=5
      else
        in=7
        open (unit=7, file=name, status='OLD', err=5000)
      endif
c
c  Change mu0
      call getone('mu_0', emuo, muo)
      if (muo.eq. 1) write(*,*) 'mu_0 changed to: ', emuo
c
c  See if we have rho or admit, period or freq or omega
      call getchr('options', kind, ikind)
      ifreq=1
      if (index(kind, 'om') .ne. 0) ifreq=2
      if (index(kind, 'pe') .ne. 0) ifreq=3
c
c  What kind of input row expected: itype 1 => rho, 2 => admit=default
      if (index(kind, 'rh') .ne. 0) itype=1
      if (index(kind, 'rh') .eq. 0) itype=2
c
      nii=5
      call getone('uncertainty', uncert, imun)
      if (imun.eq. 1) nii=3
      read (in,'(a)') row
      if (itype.eq. 1) goto 1350
c  For admittance data, look at 1st data row to determine number of cols
      do 1300 ii = nii, 1, -1
        read(row,*,end=1300, err=5100) (zz(1,i), i=1,ii)
        goto 1310
 1300 continue
 1310 nii=ii
      if (nii.eq. 3 .and. imun.ne. 1) then
        write(*,'(2a)') '>>>> No uncertainties supplied.',
     $             ' Dplus assigns 10% of |c| in Re & Im'
        uncert=0.1
      endif
c
c  Read first data row here from character string
 1350 read(row,*,err=5100) (zz(1,i), i=1, nii)
c
c  Scaling factors for input data
      call getvec('scale', scales, 5, nsc)
      nsc = min(nii, nsc)
      if (nsc.gt. 0) write(*,'(/a/6x, 1p, 5g14.4)')
     $'Each data row has been scaled by:',(scales(j),j=1,nsc)
      if (nsc.gt.1 .and. nsc.ne.nii) write(*,*)
     $'>>>> Something fishy about the scaling: not enough factors!'
c
c  Read response data into zz; first row already read in
      do 1500 j=1, maxf
        if (j.ge. 2)
     $  read (in,*,end=1600,err=5100)  (zz(j,i),i=1,nii)
        if (zz(j,1).lt. 0.0) goto 1600
        do 1450 i=1, nsc
          zz(j,i)=scales(i)*zz(j,i)
 1450   continue
        if (nii.eq. 3) zz(j,4)=sqrt(zz(j,2)**2+zz(j,3)**2) * uncert
 1500 continue
      j=1 + maxf
      write(*,*)'>>>> Too many frequencies.  Maximum =', maxf
 1600 nfreq=j-1
      nii = max(nii,4)
c
c  Closes the file and prints the data as read
      close (unit=7)
c  Negative sign in wrong column => swapped columns
      if (zz(1,4).lt. 0.0 .and. itype.eq. 2) then
        write(*,*)'>>>> Negative value in col 4: swapping cols 3 & 4'
        do 1625 i=1, nfreq
          tmp=zz(i,3)
          zz(i,3)=zz(i,4)
          zz(i,4)=tmp
 1625   continue
      endif
c
c  Rho-app data
      if (itype .eq. 1) write(*,
     $ '(/t11,a,t24,a,t36,a,t53,a,t60,a /(i5,3g14.4,2f10.3))')
     $ fmode(ifreq),'rho-app','uncertainty','phase','uncertainty',
     $ (j,(zz(j,i), i=1, nii), j=1, nfreq)
c
c  Admittances
      if (itype.eq. 2) then
        write(*,  '(/t11,a,t24,a,t36,a,t50,a,x,a)'),
     $  fmode(ifreq),'c real',' c imaginary',(' uncertainty',i=4,nii)
        do 1650 j=1, nfreq
          write(*,'(i5,1p,3g14.4,2g13.3)') j,(zz(j,i), i=1, nii)
 1650   continue
      endif
c
c
c  Construct radian frequency  omega  from given data, which may be
c  frequency omega or period
      do 2100 i=1, nfreq
        omega(i)=zz(i,1)**iover(ifreq)
        if (ifreq .ne. 2) omega(i)=2.0*pi*omega(i)
 2100 continue
c
c  Construct data set for admittances.
      if (itype .eq. 2) then
        do 3100 l=1,nfreq
         cc(1,l)=zz(l,2)
         cc(2,l)=zz(l,3)
         err(1,l)=zz(l,4)
         err(2,l)=zz(l,nii)
 3100 continue
      return
c
c  Construct real, imag parts of  c  from  rho-phase data.  Errors 
c  assumed uncorrelated in the magnitude and phase of  c . These errors
c  are saved in err, erp , both with dimensions length.
      else
        do 3500 l=1, nfreq
          cmag=sqrt(zz(l,2)/(emuo*omega(l)))
          cpha=(min(180.0-abs(zz(l,4)), abs(zz(l,4))) - 90.0)*pi/180.0
          cc(1,l)=cmag*cos(cpha)
          cc(2,l)=cmag*sin(cpha)
          err(1,l) =zz(l,3)*cmag/(2.0*zz(l,2))
          err(2,l) =zz(l,5)*cmag*radian
 3500   continue
c  List equivalent admittances if print priority high enough
        if (iprint.ge.1) write(*,
     $ '(/a/10x,a,9x,a,12x,a/(i5,1p,5g14.4))')
     $ 'internally generated equivalent complex admittance',
     $ 'omega','complex admittance','uncertainties',
     $ (i,omega(i),cc(1,i),cc(2,i),err(1,i),err(2,i), i=1, nfreq)
      endif
      return
c  Errors
 5000 write(*,*)'>>>> Nonexistent data file or you lack read permission'
      stop
c
 5100 write(*,*)'>>>> Error reading data'
      if (in.eq. 5) stop
      backspace (unit=7)
      read (7,'(a)') row
      write(*,*) '>>>> Unreadable numbers in file: '//name,
     $'Bad line: ',row
      stop
c
      end
c_______________________________________________________________________
      subroutine misfit
      implicit double precision (a-h,o-z)
c$$$$ calls forwrd
c  Computes and displays misfits and original data.  Uses /datbas/
c  to un-scale variables
      parameter (pi=3.14159265359d00,  maxf=50)
      complex ctilde
      real d1,d2,p
      common /datbas/ omega(maxf),cc(2,maxf),err(2,maxf),nfreq
      common /cases/ plist(25), chsq(25), kase
      common /option/ ifreq,itype,kount,iverb
      common /fildat/ zz(maxf,5)
      common /prolix/ iprint
      common /mks/ emuo
      dimension c(2)
      character*6 fmode(3)
      data fmode /' freq ','omega','period'/
c
      radian=pi/180.0
      obj=chsq(kase)
      chisq=obj
c
      if (itype .eq. 2 .and. iprint .ge. 1) then
c  The data are admittances
      write(*, '(/t11,a,t24,a,t36,a,t58,a,t73,a)')
     $ fmode(ifreq),'c real','c imaginary','misfits','% chisq'
      do 1500 i=1,nfreq
        call forwrd(omega(i), c)
        if (obj.gt. 0.0) then
          dr=c(1) - cc(1,i)
          di=c(2) - cc(2,i)
          trib=100.0*((dr/err(1,i))**2 + (di/err(2,i))**2)/obj
        else
         dr=0.0
         di=0.0
         trib=0.0
        endif
        write(*,'(i4,1p,3g14.4,2g13.3,0p,f6.1)') i,zz(i,1),c,dr,di,trib
 1500 continue
c
      elseif (iprint .ge. 1) then
c  The data are apparent resistivities and phases.  Complications arise
c  from the fact that admittances are the working variables.
      write(*, '(/t11,a,t24,a,t39,a,t53,a,t62,a,t73,a)')
     $ fmode(ifreq),'rho-app','misfit','phase','misfit','% chisq'
      chisq=0.0
      do 2500 i=1, nfreq
        call forwrd(omega(i), c)
        d1=c(1) - cc(1,i)
        d2=c(2) - cc(2,i)
        p=atan2(cc(2,i), cc(1,i))
        ctilde=cexp(cmplx(0.0, -p))*cmplx(d1, d2)
        trib=(real(ctilde)/err(1,i))**2 + (aimag(ctilde)/err(2,i))**2
        chisq=trib + chisq
        rho =omega(i)*emuo*(c(1)**2 + c(2)**2)
        phi = 90.0 + atan2(c(2), c(1)) /radian
        if (obj.gt. 0.0) then
          drho=rho - omega(i)*emuo*(cc(1,i)**2 + cc(2,i)**2)
          dphi=phi -  90.0 - p/radian
          trib=100.0*trib/obj
        else
          drho=0.0
          dphi=0.0
          trib=0.0
        endif
        write(*,'(i5,1p,3g14.4,0p,2f10.3,f12.1)')
     $        i,zz(i,1),rho,drho,phi,dphi,trib
 2500 continue
      endif
c
c  Print chi-squared, close the plotfile and return
      if (iprint .ge. 1) write(*,'(/a, f13.4)')' Chi-squared ',chisq
      if (itype.eq. 1 .and. abs(chisq-obj).gt. 0.01) write(*,*)
     $'Differences between chi-squared and the quadratic',
     $'objective are to be expected with rho/phase data'
      return
      end
c_______________________________________________________________________
      subroutine sample(nlam, els, klmix)
      implicit double precision (a-h, o-z)
c$$$$ calls xsort getint
c  If  klmix.eq.1  Sets up lambda samples for quadratic program,
c      getting information from keyboard
c  If  klmix.gt.1  uses solution in /soln/ to densify net in regions of
c      the non-zero components in the solution vector.
c  If  klmix.eq.0  uses solution to supply a final lambda sampling in
c      which doublets (consecutive nonzero components in the solution
c      vector) have been removed.
      parameter (maxf=50, maxunk=500)
      common /datbas/ omega(maxf),cc(2,maxf),err(2,maxf),nfreq
      common /soln/ nonzer(maxunk),al(2,maxunk),nal
      common /option/ ifreq,itype,kount,iverb
      common /slabs/ dsq,bigh0,thick  /bottom/ czero
      common /prolix/ iprint
      dimension els(*),wt(10),range(2),augmen(2)
      character*1 answer
      data madd/2/, augmen/1.0d-2,5.0d00/, wt/10*1.0d00/
c
      if (klmix.eq.0) goto 3000
      if (klmix.gt.1) goto 2000
c
c  Get an initial lambda sampling based on stock responses
c  find frequency range.  Note if zero-freq occurs
      range(2)=0.0
      do 1100 j=1,nfreq
        range(2)=max(range(2), omega(j))
 1100 continue
      range(1)=range(2)
      lbegin=3
      do 1150 j=1,nfreq
        if (omega(j).gt.0.0) range(1)=min(range(1), omega(j))
        if (omega(j).eq.0.0 .or. czero.ne.0.0) lbegin=2
 1150 continue
c  Augment actual range by a factors in  augmen
      els(2)=0.0
      els(lbegin)=range(1)*augmen(1)
      range(2)=range(2)*augmen(2)
      if (dsq.ne.0.0) range(2)=min(range(2), 1.0/dsq)
      nlam=20 + 2*nfreq
c  Print parameters and if unsatisfactory get new ones from terminal
      write(*,120) nlam,els(2),range(2),madd
 120  format(' Initial number of unknowns',i9/' lambda interval',
     $1p,2g12.4/' Refined intervals subdivided by ',i3)
      call getint('intervene', ivn, nvn)
      if (nvn .ge. 0) then
        write(*, *)'Are these to your liking? (yes or no)'
        read (*,*) answer
        if (answer .ne. 'n') goto 1400
        write(*,*) 'Enter  nlam, range(1), range(2), madd'
        read (*,*) nlam,els(2),range(2),madd
c  If lower limit is not zero use it
        if (els(2).ne.0.0) lbegin=2
c  Put logarithmically spaced samples from els(lbegin) to range(2)
 1400   continue
      endif
c
      rfact=(range(2)/els(lbegin))**(1.0/(nlam-lbegin))
      do 1500 i=lbegin,nlam
        els(i)=els(lbegin)*rfact**(i-lbegin)
 1500 continue
      if (dsq.ne.0.0) els(nlam)=0.999999/dsq
      if (iprint.ge.2) write(*,150) nlam,(els(j),j=2,nlam)
 150  format(/' Samples of lambda, n =',i6/(1p,8g9.2))
c  Set interpolation weights
      nadd=2*madd
      do 1600 j=1,madd
        wt(j)=1.0 - (2*j-1)/(nadd + 1.0)
 1600 continue
      return
c
c  When  klmix  is greater than 1, use data in /soln/ to refine net
 2000 if (2*(nal+madd)+1 .gt. maxunk) then
        write(*,*) '>>>> Refinement requires too many unknowns',
     $  '     Process skipped.'
        return
      endif
      new=nlam
      do 2200 inter=2,nal
        mid=nonzer(inter)
        if (mid.eq.2 .or. mid.eq.nlam) goto 2200
        do 2100 j=1,madd
          new=2 + new
          els(new-1)=els(mid)*wt(j) + els(mid+1)*(1.0-wt(j))
          els(new)  =els(mid)*wt(j) + els(mid-1)*(1.0-wt(j))
 2100   continue
 2200 continue
      if (nonzer(nal).eq.nlam .and. dsq.eq. 0.0) then
        new=new+1
        els(new)=2.0*els(nlam)
      endif
c  Re-order the samples to be increasing
      call xsort(els(2), new-1)
      if (iprint.ge.2) write(*,150) new,(els(j),j=2,new)
      nlam=new
      return
c
c  When klmix.eq.0  this is final refinement stage where vanishing
c  components of the vector are removed and doublets merged
 3000 if ((nonzer(2).eq.2 .and. els(2).ne.0.0) .or.
     $(nonzer(nal).eq.nlam .and. dsq.eq.0.0)) write(*,*)
     $'>>>> Caution - Expanded lambda range may yield smaller objective'
      j=1
      do 3500 new=2,nal
        j=1+j
        if (j-nal) 3100,3200,3600
 3100   if (nonzer(j+1)-nonzer(j) .eq. 1) goto 3300
 3200   els(new)=al(2,j)
        goto 3500
 3300   els(new)=(al(2,j)*al(1,j) + al(2,j+1)*al(1,j+1))/
     $ (al(1,j) + al(1,j+1))
        j=1+j
 3500 continue
      new=nal + 1
 3600 nlam=new - 1
      return
      end
c_______________________________________________________________________
      subroutine rotate(me, m, n, e)
      implicit double precision (a-h, o-z)
c$$$$ calls nothing
c  Using /datbas/ rotates matrix  e  into a frame in which the magnitude
c  and phase of  c  are uncorrelated.  Then scales with inverse std err.
      parameter (maxf=50)
      common /datbas/ omega(maxf),cc(2,maxf),err(2,maxf),nfreq
      common /option/ ifreq,itype,kount,iverb
      dimension e(me, n)
c
      if (itype .eq. 2) return
      do 1500 i=1, nfreq
      ii=2*i - 1
      cmag=sqrt(cc(1,i)**2 + cc(2,i)**2)
      cosr=cc(1,i)/cmag
      sinr=cc(2,i)/cmag
      do 1400 j=1, n
        etemp    =e(ii+1,j)*cosr - e(ii,j)*sinr
        e(ii,  j)=e(ii+1,j)*sinr + e(ii,j)*cosr
        e(ii+1,j)=etemp*err(1,i)/err(2,i)
 1400 continue
 1500 continue
      return
      end
c_______________________________________________________________________
      subroutine xsort(x, no)
      implicit double precision (a-h, o-z)
c$$$$ calls nothing
c  In-place rapid sorter.  Rearranges  x  into ascending order
      dimension x(no)
c
      mo=no
 2    if (mo-15) 21,21,23
 21   if (mo.le.1) return
      mo=2*(mo/4) + 1
      goto 24
 23   mo=2*(mo/8) + 1
 24   ko=no - mo
      jo=1
 25   i=jo
 26   if (x(i) - x(i+mo)) 28,28,27
 27   temp=x(i)
      x(i)=x(i+mo)
      x(i+mo)=temp
      i=i - mo
      if (i-1) 28,26,26
 28   jo=1 + jo
      if (jo-ko) 25,25,2
      end
c_______________________________________________________________________
      subroutine qd(n, c, ans)
      implicit double precision (a-h,o-z)
c$$$$ calls solve
      dimension ans(*),c(2,n)
c
c  Program to convert from a partial to a 1/z-type continued fraction
c  based on the Rutishauser qd algorithm
c  in the notation of gragg's notes:
c         rho(1)    ith rho
c         rho(2)    ith rho primed
c         alamb(1)  (i-1)th lambda
c         alamb(2)  (i-1)th lambda primed
c         alamb(3)  ith lambda
c         alamb(4)  ith lambda primed
c  where i is i of "do 10" loop
c
c
c  c contains the coefficients of the partial fraction on entry
c  and the coefficients of the continued fraction on exit
c  so the partial fraction coefficients are lost
c
c  alpha array is stored in c(2,i),i=1,n
c  c0 is stored in c(1,1)
c  beta array is stored in c(1,i),i=2,n
c  where n is number of terms in partial fraction
c
c  internally declared arrays:
c    alamb(4),rho(2)
c
c  Calls subroutines solve and (for n>3) update
c
      if(n.gt.1) then
      do 10 i=2,n
        call solve(0, i, c)
 10   continue
      endif
c
c  Convert alpha,beta,c0 to rho,lambda,c0 to give fraction in
c  appropriate form. Ans contains ordered c0,rho,lambda arrays
c
      ans(1)=c(1,1)
      ans(2)=c(2,1)
      if(n.eq.1)return
      do 30 i=1,n-1
        j=2*i+1
        ans(j)=c(1,i+1)/ans(j-1)
        ans(j+1)=c(2,i+1)-ans(j)
   30 continue
      return
      end
c_______________________________________________________________________
      subroutine solve(n, i, c)
      implicit double precision (a-h,o-z)
c$$$$ calls update
      dimension alamb(4),c(2,*)
c
c  Adds in one more partial fraction and computes the new c0
c  and alpha,beta arrays
c
c
c  new values of c0,alpha(1),beta(1)
c
      c0dash = c(1,1)+c(1,i)
      rho1 = c(2,1)-c(2,i)
      rho2 = c(1,1)*rho1/c0dash
      c(1,1) = c0dash
      alamb(1) = 0.d0
      if(i.gt.2)alamb(1) = c(1,2)/rho1
      alamb(2) = rho1-rho2+alamb(1)
      c(2,1) = rho2+c(2,i)
      c(1,2) = alamb(2)*rho2
c
c  Special case when i = 2
      if(i.le.2) then
        c(2,2) = alamb(2)+c(2,i)
        return
      endif
c
c  Loop through updating successive members of alpha,beta arrays
c  Special case when i = 3
      if(i.gt.3) then
        call update(i-2, alamb, c, n, i)
      endif
c
c  Calculate final beta value and last two alphas
      rho1 = c(2,i-1)-c(2,i)-alamb(1)
      rho2 = alamb(1)*rho1/alamb(2)
      alamb(4) = rho1-rho2
      c(2,i-1) = alamb(2)+rho2+c(2,i)
      c(1,i) = alamb(4)*rho2
      c(2,i) = alamb(4)+c(2,i)
      return
      end
c_______________________________________________________________________
      subroutine update(n, alamb, c, m, i)
      implicit double precision (a-h,o-z)
c$$$$ calls nothing
      dimension alamb(4),c(2,*)
c
c  New values of alpha,beta for old when another partial fraction
c  is added in. Does all but the first and last values.
c  Enter subroutine with alamb(1),alamb(2) set from calculation
c  of alpha(1),beta(1)
c
      do 10 k=2,n
        rho1=c(2,k)-c(2,i)-alamb(1)
        rho2=alamb(1)*rho1/alamb(2)
        alamb(3)=c(1,k+1)/rho1
        alamb(4)=rho1-rho2+alamb(3)
        c(2,k)=alamb(2)+rho2+c(2,i)
        c(1,k+1)=rho2*alamb(4)
        alamb(1)=alamb(3)
        alamb(2)=alamb(4)
 10   continue
      return
      end
c_______________________________________________________________________
c=======================================================================
c        Unit K:    Decoding routines for command kit
c=======================================================================
      subroutine ascan
c$$$$ calls clear getchr icheck
      implicit double precision (a-h, o-z)
c  Input routine for commands 
c
c  Reads from the standard input until eof or code "execute " or "quit".
c  Saves lines in the input store  /store/ for later retrieval by
c  getvec, getone or getchr.
c
c  Prints a glossary upon request
c
      parameter (inmx=200)
      character *80 input(inmx),line, code*4
      common /ndict/ iecho,nin
      common /store/ input
c
      if (nin .eq. 0) write(*,*) ' ',
     $'Enter commands for D+ analysis (? for help)'
c
      do 1500 l=nin+1, inmx
        read(*,'(80a)', end=3000) line
        if (line(1:4).eq.'cont' .or. line(1:4).eq.'exec') goto 2000
c  List a glossary of codes
        if (line(1:1) .eq. '?') write(*, '(a/a/(2x,a))')
     $'Enter commands from the following list:',
     $'?:                  Remind me of the command list again',
     $'execute:            Quit reading commands & begin calculations',
     $'file filename:      Enter filename of data (mandatory)',
     $'                as:  f Re c Im c err;        with c in meters',
     $'                or:  f Re c Im c Re err Im err    c in meters',
     $'                or:  f rhoa drho phi dphi    with phi in deg',
     $'options T D:        T=freq/omega/period  - type of timing',
     $'                    D=admit/rhoa - type of data ',
     $'                    defaults: T=freq D=admit',
     $'scale s1,s2,...s5:  Scale each data row by factors s1, s2, etc',
     $'uncertainty f:               Assign +-f*|c| errors to Re, Im c',
     $'mu_0  muo:          Enter a nonstandard value for mu_0', 
     $'base h1,h2, ...:   Base depths (meters) in penetration analysis',
     $'print n:            Integer n=0,1,2,3 controls amount of output',
     $'plot:               Write a plotfile for model',
     $'intervene:          User may reset lambda sampling at terminal',
     $'clear command:      Delete all earlier occurrences of command',
     $'review:             List the current command stack',
     $' '
	if (line(1:1).eq.'?') goto 1500
c
c  Review the command stack
ce
        if (line(1:4) .eq. 'revi') then
          write(*,'(5x,a)')' ', '=================== ',
     $    (input(i)(1:60),i=1,nin),'=================== '
c
        elseif (line(1:1) .ne. ' ') then
          if (line(1:4) .eq. 'echo') iecho=-iecho
          if (line(1:4) .eq. 'data') line(1:4)='file'
          if (line(1:4) .eq. 'quit') then
            write(*,'(/a)') 'Normal termination'
            stop
          endif
          nin=nin + 1
          call icheck(line)
          input(nin)=line
c  Clear a command and clear clear itself
          if (line(1:4) .eq. 'clea') then
            call getchr('clear', code, klear)
            call clear (code)
            nin=nin - 1
          endif
        endif
 1500 continue
      write(*,*) '>>>> Maximum number of commands exceeded',
     $           '     dplus quits'
      stop
c
 2000 continue
      write(*,'(5x,a)')' ', '=================== ',
     $(input(i)(1:60),i=1,nin),'=================== '
      return
c
 3000 stop
      end
c______________________________________________________________________
      subroutine icheck(line)
c$$$$ calls nothing
c  Checks the command fragment com against the catalog; appends a
c  warning if com is not present in the list.
      character*80 line, com*4
ce
      com=line(1:4)
      if (0 .eq.
     $ index('file opti clea scal mu_0 base prin plot ',com)
     $+index('inte revi unce',com))
     $ line=line(1:20)//' %<<<<<<< Unrecognized command'
      return
      end
c______________________________________________________________________
      subroutine getvec(code, values, nwant, nfound)
c$$$$ calls ljust
      implicit double precision (a-h, o-z)
c
c  Extracts a vector of numbers from the input store.  That store is
c  a large character array in /store/ which has been filled earlier
c  by  ascan.
c
c  code    A 4-bye identifying code.  Only lines in the input store
c          beginning with this code are scanned for numbers.
c  values  the real output array of values found.
c  nwant   the maximum number of numbers expected .
c  nfound  the number of numbers actually found in the input.
c          If the line contains fewer than nwant  values, this is the
c          value returned in  nfound.  If an error is discovered
c          nfound=-n, where  n  is the number of numbers properly
c          decoded.  If there are no numbers after the codeword
c          nfound=0.  Finally, if the code is absent from the store
c          nfound=-99 and the array is left undisturbed.
c
      parameter (inmx=200)
      character *80 input(inmx),line,local,char, code*4
      common /store/ input
      common /ndict/ iecho,nin
      dimension values(*)
c
c  Read the store in reverse order (Thus last entry is obeyed)
ce
      do 1010 lin=nin, 1, -1
        line=input(lin)
c  Check for code 
        if (code .eq. line(1:4)) goto 1020
 1010 continue
c  Code word not found
      nfound=-99
      return
c
 1020 if (iecho .ge. 1) write(*,'(2a)')'==> ',line
      n1=index(line, ' ')+1
      n2=index(line, '%')
      n2=80 + min(n2,1)*(n2 - 81)
      char=line(n1:n2)
c
      do 1500 n=1, nwant
        call ljust(char, local)
        lbl=index(local, ' ')
        if (lbl .eq. 1) goto 1510
        read (local, *, err=2000) values(n) 
        char=local(lbl:80)
 1500 continue
      n=nwant+1
 1510 nfound=n - 1
      return
c
 2000 write(*,*)' ','>>>> Unreadable numbers in this input line:',line
      nfound=1 - n
      return
      end
c______________________________________________________________
      subroutine getone(code, value, nfound)
c$$$$ calls getvec
      implicit double precision (a-h, o-z)
c
c  Extracts a single number from the input store.  That store is
c  a large array in common /store/ which has been filled earlier.
c
c  code    A 4-bye identifying code.  Only lines in the input store
c          beginning with this code are scanned for numbers.
c  value   the real output variable containing the desired number.
c  nfound  is 1 if a number is successfully read in; it is zero
c          the number is absent or unreadable.  nfound = -99 if the
c          code is absent from the input store.
c
      character*4 code
      dimension v(1)
ce
      call getvec(code, v, 1, nfound)
      if (nfound .eq. 1) value=v(1)
      return
      end
c______________________________________________________________
      subroutine getint(code, number, nfound)
c$$$$ calls getvec
      implicit double precision (a-h, o-z)
c
c  Extracts a single integer from the input store.  
c  See getone for details.
c
      character*4 code
      dimension v(1)
ce
      call getvec(code, v, 1, nfound)
      if (nfound .eq. 1) then
        number=nint(v(1))
        if (number .ne. v(1)) print '(4a,1p,g16.8/a,i10)',
     $  '>>>> Warning: ',code,' expects an integer argument, but',
     $  ' found: ',v(1),'    Returns: ',number,' '
      endif
      return
      end
c______________________________________________________________
      subroutine getchr(code, char, nbytes)
c$$$$ calls ljust
      implicit double precision (a-h, o-z)
c  Extracts a character variable from the input store.  That 
c  store is a large array in common /store/ which has been filled 
c  earlier.
c
c  code    A 4-bye identifying code.  Only lines in the input store
c          beginning with this code are scanned.
c  char    the character output variable containing the desired string.
c  nbytes  is the length of the string read; it is zero if
c          the line is blank after the code.  nbytes = -99 if the
c          code is absent from the input store.
c
      parameter (inmx=200)
      character *80 input(inmx),line, char*(*),  code*4
      common /store/ input
      common /ndict/ iecho,nin
c  Inspect the store in reverse order (thus reading latest entry)
ce
      do 1010 lin=nin, 1, -1
        line=input(lin)
c  Check for code word  
        if (code .eq. line(1:4)) goto 1020
 1010 continue
c  Code word not found
      nbytes=-99
      return
c
 1020 if (iecho .ge. 1) write(*,'(2a)')'==> ',line
c  
      n1=index(line, ' ')+1
      n2=index(line, '%')
      n2=80 + min(n2,1)*(n2 - 81)
c  Check for blank string
      nbytes=0
      if (line(n1:n2) .eq. ' ') return
c  Save in char everything from 1st non-blank to last non-blank
      call ljust(line(n1:n2), char)
      nn=min(n2-n1+1, len(char))
      do 1200 k=nn, 1, -1
        if (char(k:k) .ne. ' ') goto 1210
 1200 continue
 1210 nbytes= k
      return
      end
c______________________________________________________________
      subroutine clear(code)
c$$$$ calls no other routines
      implicit double precision (a-h, o-z)
c
c  Removes the last command entry identified by code.
c  code    A 4-byte identifying code.  Only lines in the input store
c          beginning with this code are scanned.
c
      parameter (inmx=200)
      character *80 input(inmx),line, code*4
      common /store/ input
      common /ndict/ iecho,nin
c
c  Inspect the store in normal order
ce
      do 1010 lin= nin, 1, -1
        line=input(lin)
c  Check code  and clear the line if present, but only once
        if (code .eq. line(1:4)) then
          input(lin)=' --'
        endif
 1010 continue
      return
      end
c______________________________________________________________
      subroutine ljust(str1, str2)
c$$$$ calls nothing
c  Left justify: input character string str1 (up to 80 bytes in length)
c  is left justified, removing leading blanks, and returned in str2.
c  str1 and str2 may be set to the same variable in the call.
      character*(*) str1, str2, str3*80
c
ce
      str2=str1
      if (str1 .eq. ' ') return
c
      l1=len(str1)
      str3=str1
      do 1010 j=1, l1
        if (str1(j:j) .ne. ' ') goto 1011
 1010 continue
 1011 str2=str3(j:l1)
      return
      end
c______________________________________________________________
      blockdata iounit
      implicit double precision (a-h, o-z)
      parameter (inmx=200)
      common /prolix/ iprint
      common /ndict/ iecho,nin
      data iprint /1/
c  Don't echo command lines when read
      data iecho/-1/
c  Initial command line count
      data nin/0/
      end
c______________________________________________________________________
      subroutine nnls(a, mdim, m, n, b, x, rnorm, w, zz, index, mode)
      implicit double precision (a-h,o-z)
c$$$$ calls h12 g1 g2
c  Non-negative least-squares fitting routine.  Given by Lawson & Hanson
c  in  `Solving Least Squares Problems ,' Prentice-Hall 1974.
c  Given  matrix  a  (m by n), and an m-vector  b  , finds the n-vector
c  x  that solves the problem 2-norm(a*x - b) minimum with x.ge.0.
c
c  a(),mdim,m,n  a  is the array, with mrows,  n  cols.  mdim  is the
c     actual dimension given  a  in the calling program (mdim.ge.m)
c     a  is destroyed by program.
c  b()  on entry must contain m-vector  b.  Destroyed by program
c  x()  solution vector  x.  Need not be initialized.
c  rnorm  minimum two-norm.
c  w()  an array of working space.  Must be dimensioned at least n
c  zz()  another working array.  Must be dimensioned at least m.
c  index  another working array. Must be dimensioned at least  n.
c  mode  a flag indicateng outcome of subroutine.
c     mode=1   means routine worked properly.
c     mode=2   means dimensions of arrays were bad (m.le0 .or. n.le.0)
c     mode=3   means iteration count exceeded (more than 3*n iterations)
c
      dimension a(mdim,n),b(m),x(n),w(n),zz(m),index(n),dummy(1)
c
      zero=0.
      one=1.
      two=2.
      factor=.01
c
      mode=1
      if (m.gt.0 .and. n.gt.0) go to 10
      mode=2
      return
 10   iter=0
      itmax=3*n
c
c
c  Initialize arrays index  and  x
      do 20 i=1,n
      x(i)=zero
 20   index(i)=i
c
      iz2=n
      iz1=1
      nsetp=0
      npp1=1
c     Main loop begins here
 30   continue
c
c  Quit if all coeffs are already in solution. Or if m cols of a have
c  been triangularized.
      if (iz1.gt.iz2 .or. nsetp.ge.m) go to 350
c
c
c  Compute components of dual vector w
      do 50 iz=iz1,iz2
      j=index(iz)
      sm=zero
      do 40 l=npp1,m
 40   sm=a(l,j)*b(l) + sm
 50   w(j)=sm
c  Find largest +ve w(j)
 60   wmax=zero
      do 70 iz=iz1,iz2
      j=index(iz)
      if (w(j).le.wmax) go to 70
      wmax=w(j)
      izmax=iz
 70   continue
c
c
c
c  If wmax.le.0 go to termination (kuhn-tucker conditions ok)
      if (wmax) 350,350,80
 80   iz=izmax
      j=index(iz)
c
c
c
c  Sign of w(j) is ok for j to be moved to set  p.
c  Begin Householder trans. Check new diagonal element to avoid near ld.
      asave=a(npp1,j)
      call h12(1, npp1, npp1+1, m, a(1,j), 1, up, dummy, 1, 1, 0)
      unorm=zero
      if (nsetp.eq.0) go to 100
      do 90 l=1,nsetp
 90   unorm=a(l,j)**2 + unorm
 100  unorm=sqrt(unorm)
      temp=unorm + abs(a(npp1,j))*factor
      if (temp-unorm) 130,130,110
c  Col j is sufficiently indep. Copy  b  to zz. update  z  and solve
c  for  ztest (=proposed new val for x(j)).
c
 110  do 120 l=1,m
 120  zz(l)=b(l)
      call h12(2, npp1, npp1+1, m, a(1,j), 1, up, zz, 1, 1, 1)
      ztest=zz(npp1)/a(npp1,j)
c
c
c              See if  ztest is +ve
      if (ztest) 130,130,140
c
c  Reject  j  as a candidate to be moved from set  z  to set  p.
c  Restore a(npp1,j), set w(j)=0. And loop back to test dual
c  coeffs again.
c
 130  a(npp1,j)=asave
      w(j)=zero
      go to 60
c
c  The index j=index(iz) has been selected to be moved ftom set  z
c  to set  p .  Update b,  update indices, apply householder trans
c  to cols in new set z.  Zero subdiagonal  elts in col  j, set
c  w(j)=0.
c
 140  do 150 l=1,m
 150  b(l)=zz(l)
c
      index(iz)=index(iz1)
      index(iz1)=j
      iz1=1 + iz1
      nsetp=npp1
      npp1=1 + npp1
c
      if (iz1.gt.iz2) go to 170
      do 160 jz=iz1,iz2
      jj=index(jz)
 160  call h12(2, nsetp, npp1, m, a(1,j), 1, up, a(1,jj), 1, mdim, 1)
 170  continue
c
      if  (nsetp.eq.m) go to 190
      do 180 l=npp1,m
 180  a(l,j)=zero
 190  continue
c
      w(j)=zero
c  Solve triangulate system. Store sol in zz,temporarily
c
      next=1
      go to 400
 200  continue
c
c  Secondary loop begins here
 210  iter=1 + iter
      if (iter.le.itmax) go to 220
      mode=3
c  Write statement here deleted*********************
      go to 350
 220  continue
c
c
c  See if all new constrained coeffs are feasible. If not, find alpha.
      alpha=two
      do 240 ip=1,nsetp
      l=index(ip)
      if (zz(ip)) 230,230,240
c
 230  t=-x(l)/(zz(ip) - x(l))
      if (alpha.le.t) go to 240
      alpha=t
      jj=ip
 240  continue
c  If all new constrained coeffs are feasible, alpha is still 2. If so,
c  exit into main loop.
c
      if (alpha.eq.two) go to 330
c  Otherwise alpha will be in (0,1) to interpolate between old x
c  and new zz.
c
      do 250 ip=1,nsetp
      l=index(ip)
 250  x(l)=alpha*(zz(ip)-x(l)) + x(l)
c
c  Modify  a  bb  and the index arrays to move coeff  i  from  set  p
c   to set  z.
c
      i=index(jj)
 260  x(i)=zero
c
      if (jj.eq.nsetp) go to 290
      jj=1 + jj
      do 280 j=jj,nsetp
      ii=index(j)
      index(j-1)=ii
      call g1(a(j-1,ii), a(j,ii), cc, ss, a(j-1,ii))
      a(j,ii)=zero
      do 270 l=1,n
      if (l.ne.ii) call g2(cc, ss, a(j-1,l), a(j,l))
 270  continue
 280  call g2(cc, ss, b(j-1), b(j))
 290  npp1=nsetp
      nsetp=nsetp - 1
      iz1=iz1 -     1
      index(iz1)=i
c  All coeffs in set p should be feasible.  If they are not, this is due
c  to round-off.  Non-positive ones set to 0 and moved to set z.
c
c
      do 300 jj=1,nsetp
      i=index(jj)
      if (x(i)) 260,260,300
 300  continue
c
c
c  Copy  b  to  zz  and solve again and loop back
      do 310 i=1,m
 310  zz(i)=b(i)
      next=2
      go to 400
 320  continue
      go to 210
c
c  End of secondary loop
 330  do 340 ip=1,nsetp
      i=index(ip)
 340  x(i)=zz(ip)
c  All new coeffs are +ve.  Loop back to beginning.
      go to 30
c  End of main loop.  Terminating section.
c
c
 350  sm=zero
      if (npp1.gt.m) go to 370
      do 360 i=npp1,m
 360  sm=b(i)**2 + sm
      go to 390
 370  do 380 j=1,n
 380  w(j)=zero
 390  rnorm=sqrt(sm)
      return
c
c  Solution of linear triangular system
c
 400  do 430 l=1,nsetp
      ip=nsetp + 1 - l
      if (l.eq.1) go to 420
      do 410 ii=1,ip
 410  zz(ii)=zz(ii) - a(ii,jj)*zz(ip+1)
 420  jj=index(ip)
 430  zz(ip)=zz(ip)/a(ip,jj)
      go to (200, 320), next
c ******** format deleted for print statement here ****************
      end
c_______________________________________________________________________
      subroutine g2(cosa, sina, x, y)
      implicit double precision (a-h,o-z)
c$$$$ calls no other routines
c  Rotates a 2-vector (x,y) by angle a. Over-writes original vector.
c
      z=cosa*x + sina*y
      y=cosa*y - sina*x
      x=z
      return
      end
c_______________________________________________________________________
      subroutine g1(a, b, cosa, sina, sig)
      implicit double precision (a-h,o-z)
c$$$$ calls no other routines
c
c  From Lawson & Hanson `Solving Least Squares Problems' 1974 p 309.
c  Computes an orthogonal matrix to rotate (a,b) into (sqrt(a*a+b*b),0)
c  matrix in form (cosa, sina/-sina, cosa) and sig holds magnitude of
c  (a,b).  Sig is allowed to overwrite a or b in calling routine.
c
      zero=0.0
      one=1.0
      if (abs(a) .le. abs(b)) go to 10
      xr=b/a
      yr=sqrt(one + xr**2)
      cosa=sign(one/yr, a)
      sina=cosa*xr
      sig=abs(a)*yr
      return
 10   if (b) 20,30,20
 20   xr=a/b
      yr=sqrt(one + xr**2)
      sina=sign(one/yr, b)
      cosa=sina*xr
      sig=abs(b)*yr
      return
 30   sig=zero
      cosa=zero
      sina=one
      return
      end
c_______________________________________________________________________
      subroutine h12(mode, lpivot, l1, m, u, iue, up, c, ice, icv, ncv)
      implicit double precision (a-h,o-z)
c$$$$ calls no other routines
c  Construction and application of householder transformation for nnls
c  from  h12  in Lawson & Hanson `Solving Least Squares Problems'
c  - Prentice-Hall 1974 (pp308,309).  Double precision throughout.
      dimension u(iue,m),c(*)
      if (0.ge.lpivot .or. lpivot.ge.l1 .or. l1.gt.m) return
      cl=abs(u(1,lpivot))
      if (mode.eq.2) go to 60
c  Construct transformation
      do 10 j=l1,m
 10   cl=max(abs(u(1,j)),cl)
      if (cl) 130,130,20
 20   clinv=1.0/cl
      sm=(u(1,lpivot)*clinv)**2
      do 30 j=l1,m
 30   sm=(u(1,j)*clinv)**2 + sm
      sm1=sm
      cl=sqrt(sm1)*cl
      if (u(1,lpivot)) 50,50,40
 40   cl=-cl
 50   up=u(1,lpivot) - cl
      u(1,lpivot)=cl
      go to 70
c  Apply transformation
 60   if (cl) 130,130,70
 70   if (ncv.le.0) return
      b=up*u(1,lpivot)
      if (b) 80,130,130
 80   b=1.0/b
      i2=1-icv+ice*(lpivot-1)
      incr=ice*(l1-lpivot)
      do 120 j=1,ncv
      i2=icv + i2
      i3=i2 + incr
      i4=i3
      sm=c(i2)*up
      do 90 i=l1,m
      sm=c(i3)*u(1,i) + sm
 90   i3=ice + i3
      if (sm) 100,120,100
 100  sm=b*sm
      c(i2)=sm*up + c(i2)
      do 110 i=l1,m
      c(i4)=sm*u(1,i) + c(i4)
 110  i4=ice + i4
 120  continue
 130  return
      end 
c_______________________________________________________________________
      subroutine plot(nal, d,tau,bot, skilo)
c$$$$ calls getchr getint forwrd
c  Writes a plot file to fort.19 
      implicit double precision (a-h,o-z)
      parameter (maxf=50)
      character*20 frmt,bk*1,sqchi*13,idfile*60
      common /datbas/ omega(maxf),cc(2,maxf),err(2,maxf),nfreq
      common /cases/ plist(25), chisq(25), kase
      common /prolix/ iprint
      dimension d(*),tau(*), ctheo(2,5*maxf),om(5*maxf)
      data ngraph/0/
c
c  Write a plotfile of this model. Its name is fort.19
      call getint('plot', ignor, notice)
      if (notice.eq. -99) return
      call getchr('file', idfile, ignor)
c
      taulo=tau(1)
      tauhi=taulo
      do 2100 j=2, nal-1
*       if (tau(j).le.1d-10 .or.  tau(j) .gt.1d10) goto 2200
        taulo=min(taulo, tau(j))
        tauhi=max(tauhi, tau(j))
 2100 continue
      iall=nal-1
 2200 taulo=10d0**(int(log10(taulo)+1000.7) - 1001)
c
      write(*,*) 'Plotxy file written to fort.19'
c
c  Write instructions for lin-lin plot of model
c
      write(19,'(a)') 'weight 6 12','char 0.08','xlim 3.5','ylim 2.5',
     $'xlab Depth, km','ylab Conductance, S','frame gridsol'
      write(19,'(a,i6)')'read ',3*iall+1
 210  format(2g15.5)
      write(19,210)
     $ (skilo*d(j),taulo,skilo*d(j),tau(j),skilo*d(j),taulo, j=1,iall)
      write(19,210) skilo*bot,taulo
c  Perfect conductor
      if (tau(nal).lt. 0.0d0 ) then
        z1=skilo*1.08*bot
        write(19,'(a/a/a/(2g12.4))')'fill','color pale gray','read 4',
     $  bot*skilo,taulo, bot*skilo,0.9*tauhi,z1,0.9*tauhi,z1,taulo
        write(19,'(a)')
     $ 'color black','char 0.08 90',
     $ 'note (3.3 0.3 i) Perfect conductor','char 0.08 0'
      endif
      write(19,'(a)')'plot 1.2 1.7',' '
c
c  Next log-log plot of model
c  First, perfect conductor`
      if (tau(nal).lt. 0.0d0 ) then
        zmin=d(1)+d(2)*(0.5d0+sign(0.5d0, 1.0d-8 -d(1)))
        z1= skilo*zmin*exp(3.5*log(bot/zmin)/3.2)
        write(19,'(a/a/a/(2g12.4))')'fill','color pale gray','read 4',
     $  bot*skilo,taulo, bot*skilo,0.9*tauhi,z1,0.9*tauhi,z1,taulo
        write(19,'(a)')
     $ 'color black','char 0.08 90',
     $ 'note (3.3 0.3 i) Perfect conductor','char 0.08 0'
      endif
      j1=1
      if (d(1).eq. 0.0) then
        skd1=skilo*d(2)/10.0
        write(19,'(a,2g12.4,a)')'note (',skd1,taulo/1.3,' c)z=0'
        write(19,'(a/(2g15.5))')  'read 2',skd1,taulo,skd1,tau(1)
        j1=2
      endif
      write(19,'(a,i6)')'read ',3*iall+1-3*(j1-1)
      write(19,210)
     $ (skilo*d(j),taulo,skilo*d(j),tau(j),skilo*d(j),taulo, j=j1,iall)
      write(19,210) skilo*bot,taulo
      write(19,'(a)')'logxy loglog','plot',' '
c
c  Plot the data and model response
      write(19,'(a)') 'xlab Frequency, Hz','ylab Admittance c*, meters',
     $'file *','logxy loglog', 'mode 3','symbol circle 0.02'
      write(19,'(a,i6/(3g12.4))') 'read ',nfreq,
     $(omega(j)/6.283,cc(1,j),err(1,j),j=1,nfreq)
      write(19,'(a,i6/(3g12.4))') 'read ',nfreq,
     $(omega(j)/6.283, -cc(2,j),err(2,j),j=1,nfreq)
      i=0
      do 3200 j=1, nfreq-1
        do 3100 k=1, 5
          i=i + 1
          om(i) = 0.25*((k-0.96)*omega(j+1) + (5-k)*omega(j))
          call forwrd(om(i), ctheo(1,i))
 3100   continue
 3200 continue
c  Backslash is char(92)
      bk = char(92)
      sqchi=bk//'chi'//bk//bk//'sup{2}='
      ndig=4 + max(0, int(log10(chisq(kase)+0.01)))
      write(frmt,'(a,i5,a)') '(a/2a,f',ndig,'.1)'
      write(19,frmt)'note','note (3.3 2.2ir)',sqchi,chisq(kase)
c
      write(19,'(a)')'background white notes'
      write(19,'(a)')'smooth','mode 2','affine 0.1592 0 1 0'
      write(19,'(a/a,i6/(2g12.4))') 'color red',
     $'read ',i,(om(k),ctheo(1,k),k=1,i)
      write(19,'(a)')'note (vr) Re c'
      write(19,'(a/a,i6/(2g12.4))') 'color blue',
     $'read ',i,(om(k),-ctheo(2,k),k=1,i)
      write(19,'(a)')'note (vr) -Im c','color blac'
      write(19,'(2a)')'title From ',idfile
      write(19,'(a)')'plot',' ','stop'
c
      ngraph = ngraph + 1
      write(19,'(i5,a)') ngraph,' --------------------------'

      return
      end
c_______________________________________________________________________