      subroutine solve(meflag)
c
c ...... solve equations of condition using weighted least-squares
c
      implicit real*8 (a-h,o-z)
      character*1 polar
      character*72 id
      character*4 stn
      integer*4 meflag
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common/obs/pt(1000),st(1000),nsreq
      common/obs1/stn(1000),polar(1000)
      common/stn/syf(1000),sxf(1000),psaf(1000),ssaf(1000),nstns
      common/epid/id
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common/condn/cn
      common/dim/m,n,nmax,nit
      common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd
      common/c4/cov44
      common/mmax/mmax
c
      if(nit.gt.1) go to 2
      sum=0.d0
      do 1 i=1,m
      sum=sum+b(i)*b(i)
    1 continue
      if((sum.gt.0.d0).and.(m.gt.0)) then
        nnn=15
        sem1=dsqrtf(nnn,sum/(dreal(m)-1.d0))
      else
        sem1=0.d0
      endif
    2 continue
c
c      write(*,*) ' m & n = ',m,n
c
      do 4 i=1,m
      if(w(i).gt.0.d0) then
        nnn=16
        w2=dsqrtf(nnn,w(i))
      else
        w2=0.
      endif
c ...... construct weighted equations of condition -> w A x = w b
      do 3 j=1,n
      a(i,j)=w2*a(i,j)
c
c      write(*,*) 'i,j,w2,a(i,j) = ',i,j,w2,a(i,j)
c
    3 continue
      b(i)=w2*b(i)
c
c      write(*,*) 'i,w2,b(i) = ',i,w2,b(i)
c
    4 continue
c
      call svdcom(a,m,n,mmax,nmax,q,u,v)
c
c      write(*,*) 'n,q(i) = ',n,(q(kl),kl=1,n)
c
      call damper(a,b,w,m,n,mmax,q,u,v,qi,x,bmax,var,cn,meflag,nit)
      do 45 ij=1,n
c      write(*,*) ',st_slv - i & x = ',ij,x(ij)
   45 continue
      call infoco(m,n,u)
      if(var.gt.0.d0) then
        nnn=17
        se=dsqrtf(nnn,var)
      else
        se=0.d0
      endif
      do 7 i=1,n
      do 6 j=1,n
      sum=0.d0
      do 5 k=1,n
      sum=sum+v(i,k)*v(j,k)*var*qi(k)*qi(k)
    5 continue
      cov(i,j)=sum
    6 continue
    7 continue
      if(n.eq.4) then
        cov44=cov(4,4)
      else
        cov44=0.d0
      endif
      do 30 ii=1,n
c      write(*,*) 'solve - cov = ',ii,(cov(ii,jj),jj=1,n)
   30 continue
c
      return
      end

      subroutine conver
c
c ...... test for convergence
c
      implicit real*8 (a-h,o-z)
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common/dim/m,n,nmax,nit
      integer*4 swhdf,swhdx,swnit,swnmx,swcon
      common/swtch/swhdf,swhdx,swnit,swnmx,swcon
      common/step/step
      sum=0.d0
      do 1 i=1,n
      sum=sum+x(i)*x(i)
    1 continue
c ...... step = magnitude of change vector
      if(sum.gt.0.d0) then
        nnn=18
        step=dsqrtf(nnn,sum)
      else
        step=0.d0
      endif
c ...... set convergence switch if step size < p_step
c
c ...... modified convergence criteria for historical earthquakes
c
	if(nit.le.10) then
	  p_step=0.025d0+0.0075d0*dble(nit)
	elseif(nit.le.20) then
	  p_step=0.1d0+0.03d0*dble(nit-10)
	elseif(nit.le.30) then
	  p_step=0.4d0+0.12d0*dble(nit-20)
	else
	  p_step=1.6d0+0.48d0*dble(nit-30)
	endif
c
      if(step.lt.p_step) then
        swcon=1
      endif
c      write(*,*) 'conver: nit & step = ',nit,'   ',step
c
      return
      end

      subroutine pertur(n_zdse)
c
c ...... perturb solution
c
c        perturbations -
c            x(1) - latitude step (km) (+=north)
c            x(2) - longitude step (km) (+=east)
c            x(3) - origin time step (sec) (+=forward)
c            x(4) - hypocentral depth step (km) (+=down)
c
      implicit real*8 (a-h,o-z)
      character*72 id
      common/epid/id
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd
      common/dim/m,n,nmax,nit
      common/c4/cov44
      integer*4 swhdf,swhdx,swnit,swnmx,swcon,n_zdse
      common/swtch/swhdf,swhdx,swnit,swnmx,swcon
      real*8 kpd
      kpd=111.195d0
      hstp2=x(1)*x(1)+x(2)*x(2)
      if(hstp2.gt.0.d0) then
        nnn=19
        hstp=dsqrtf(nnn,hstp2)
      else
        hstp=0.d0
      endif
c ...... set switch if horizontal step < 1.0 km
      if(n_zdse.eq.0) then
        if(hstp.le.1.0d0) then
          swhdx=1
        endif
      endif
c
c ...... inhibit large perturbations
c
      ext=ex
      if(dabs(x(1)).gt.8.d0) then
        dex=dsign(8.d0,x(1))
      else
        dex=x(1)
      endif
      ex=ex+dex/(kpd*dcosd(ey))
      eyt=ey
      if(dabs(x(2)).gt.8.d0) then
        dey=dsign(8.d0,x(2))
      else
        dey=x(2)
      endif
      ey=ey+dey/kpd
      ett=et
      if(dabs(x(3)).gt.1.d0) then
        det=dsign(1.d0,x(3))
      else
        det=x(3)
      endif
      et=et+det
      dez=0.d0
      if(n.eq.4) then
c
c        dez=x(4)
c        ez=ez+dez
c
        dez_max=0.1d0*ez
        if(dez_max.lt.0.1d0) then
          dez_max=0.1d0
        elseif(dez_max.gt.4.d0) then
          dez_max=4.d0
        endif
        dez=x(4)
        if(dabs(dez).gt.dez_max) dez=dsign(1.d0,dez)*dez_max
        ez=ez+dez
c
c ...... restrict hypocentral depth if:
c            1 - ez.lt.3
c            2 - ez.gt.15
c
        if(ez.le.3.d0) then
          ez=3.d0
          edepth=ez
          x(4) = 0.d0
          cov(4,4) = 0.d0
          n = 3
          swnmx=3
          ex=ext
          ey=eyt
          et=ett
          dez=0.d0
        elseif(ez.gt.15.d0) then
          if(m.gt.12) then
            if(ez.gt.15.d0) then
              ez=15.d0
              edepth=ez
              x(4) = 0.d0
              cov(4,4) = 0.d0
              n = 3
              swnmx=3
              ex=ext
              ey=eyt
              et=ett
              dez=0.d0
            endif
          else
            ez=15.d0
            edepth=ez
            x(4) = 0.d0
            cov(4,4) = 0.d0
            n = 3
            swnmx=3
            ex=ext
            ey=eyt
            et=ett
            dez=0.d0
          endif
        endif
      endif
c      if(n.lt.4) then
c        write(*,*) 'pertur - ex, ey & et = ',ex,ey,et
c      else
c        write(*,*) 'pertur - ex, ey, ez & et = ',ex,ey,ez,et
c      endif
c
      return
      end

      subroutine partia (ex,ey,ez,sx,sy,va,vb,vc,vd,dtdx,dtdy,dtdt,
     1                   dtdz)
c
c ...... compute partial derivatives
c
      implicit real*8 (a-h,o-z)
c
c ...... computes partial derivatives with respect to model parameters
c        (dtda etc.) which are not used in this program
c
c        step = ray shooting radius of convergence
c	 06 mar 95 - Changed zhgtdpf routine to expedite convergence.
c
      call dazd(ex,ey,sx,sy,d,az)
      step=1.d-4
      z0=0.d0
      call ptdcal (va,vb,vc,vd,z0,ez,step,d,p,t,dtda,dtdb,dtdc,dtdd,
     1             dtdz,dtdw,dtdt,hemi,aoi)
      dtdx=-1.0d0*dtdw*dsind(az)
      dtdy=-1.0d0*dtdw*dcosd(az)
c
      return
      end

      subroutine abcomp
c
c ...... construct equations of condition
c
      implicit real*8 (a-h,o-z)
      character*72 id
      character*1 polar,iwtp,iwts,iwt
      character*4 stn
      common/files/ ncr,nsp,ns1,ns2,ns3,ns4
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common/obs/pt(1000),st(1000),nsreq
      common/obs1/stn(1000),polar(1000)
      common/stn/syf(1000),sxf(1000),psaf(1000),ssaf(1000),nstns
      common/epid/id
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common/dim/m,n,nmax,nit
      common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd
      common/dstnc/dist(1000)
      common/wght/iwtp(1000),iwts(1000),iwt(1000)
      m = 0
      do 1 i=1,nsreq
      if (pt(i).ne.0.d0) then
        m = m + 1
        sy=syf(i)
        sx=sxf(i)
        psa=psaf(i)
        ssa=ssaf(i)
        call partia(ex,ey,ez,sx,sy,vpa,vpb,vpc,vd,a(m,1),a(m,2),
     1              a(m,3),a(m,4))
        call tcal(ex,ey,ez,sx,sy,vpa,vpb,vpc,vd,et,cpt,hemi,aoi)
        b(m)=pt(i)-cpt-psa
c        if(b(m).lt.-600.d0) then
c          b(m)=b(m)+86400.d0
c        elseif(b(m).gt.+600.d0) then
c          b(m)=b(m)-86400.d0
c        endif
c        write(*,*) 'abcomp - ex, ey, ez & et = ',ex,ey,ez,et
c        write(*,*) 'abcomp - sx & sy = ',sx,sy
c        write(*,*) 'abcomp - P - ',m,a(m,1),a(m,2),a(m,3),a(m,4),
c     1             b(m)
        iwt(m)=iwtp(i)
        call dazd(ex,ey,sx,sy,dist(m),dmy)
      endif
      if (st(i).ne.0.d0) then
        m = m + 1
        sy=syf(i)
        sx=sxf(i)
        psa=psaf(i)
        ssa=ssaf(i)
        call partia(ex,ey,ez,sx,sy,vsa,vsb,vsc,vd,a(m,1),a(m,2),
     1              a(m,3),a(m,4))
        call tcal(ex,ey,ez,sx,sy,vsa,vsb,vsc,vd,et,cst,hemi,aoi)
        b(m)=st(i)-cst-ssa
c        if(b(m).lt.-600.d0) then
c          b(m)=b(m)+86400.d0
c        elseif(b(m).gt.+600.d0) then
c          b(m)=b(m)-86400.d0
c        endif
c        write(*,*) 'abcomp - ex, ey, ez & et = ',ex,ey,ez,et
c        write(*,*) 'abcomp - sx & sy = ',sx,sy
c        write(*,*) 'abcomp - S - ',m,a(m,1),a(m,2),a(m,3),a(m,4),
c     1             b(m)
        iwt(m)=iwts(i)
        call dazd(ex,ey,sx,sy,dist(m),dmy)
      endif
    1 continue
c      write(*,*) 'm = ',m
c
      return
      end

      subroutine ordst(m,qd1,qd3)
c
c ...... order by distance
c
      implicit real*8 (a-h,o-z)
      common/dstnc/dist(1000)
      real*8 ord(1000)
      do 1 i=1,m
      ord(i)=dist(i)
    1 continue
      m1=m-1
      do 3 j=1,m1
      j1=j+1
      do 2 k=j1,m
      if(ord(k).le.ord(j)) then
        ot=ord(j)
        ord(j)=ord(k)
        ord(k)=ot
      endif
    2 continue
    3 continue
      m1=m/4+1
      m3=(3*m)/4+1
      qd1=ord(m1)
      qd3=ord(m3)
      if(qd1.lt.50.d0) qd1=50.d0
      if(qd3.lt.200.d0) qd3=200.d0
c
      return
      end

      subroutine ordrs(m,nrt,qb1,qb3)
c
c ...... order by residual
c
      implicit real*8 (a-h,o-z)
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      real*8 ord(1000)
      if(nrt.gt.1) go to 2
      do 1 i=1,m
      ord(i)=dabs(b(i))
    1 continue
      go to 4
    2 continue
      do 3 i=1,m
      ord(i)=dabs(bmax(i))
    3 continue
    4 continue
      m1=m-1
      do 6 j=1,m1
      j1=j+1
      do 5 k=j1,m
      if(ord(k).le.ord(j)) then
        ot=ord(j)
        ord(j)=ord(k)
        ord(k)=ot
      endif
    5 continue
    6 continue
      m1=m/4+1
      m3=(3*m)/4+1
      qb1=ord(m1)
      qb3=ord(m3)
      if(qb1.lt.0.5d0) qb1=0.5d0
      if(qb3.lt.2.0d0) qb3=2.0d0
c
      return
      end

      subroutine weight(meflag,mwflag)
c
c ...... calculate weight for each equation of condition
c
      implicit real*8 (a-h,o-z)
      character*1 iblnk,izero,ione,itwo,ithree,ifour,iwtp,iwts,iwt
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common/dim/m,n,nmax,nit
      common/dstnc/dist(1000)
      common/wght/iwtp(1000),iwts(1000),iwt(1000)
      integer*4 meflag,mwflag
      data iblnk,izero,ione,itwo,ithree,ifour/' ','0','1','2','3','4'/
c
      i1=1
      i2=2
      if(nit.eq.1) then
        sem1=10.d0
        se=10.d0
      endif
c
c ...... fix weighting after 8 interations
c
	if(nit.ge.8) return
c
      wmin=1.e-8
      fm=m
      if(mwflag.eq.1) then
        do 16 i=1,m
        w(i)=1.d0
   16   continue
        return
      endif
      if(meflag.eq.1) go to 8
      if (nit.gt.1) go to 8
      sum = 0.d0
      do 1 i=1,m
      sum = sum + b(i)*b(i)
    1 continue
      if((sum.gt.0.d0).and.(m.gt.0)) then
        nnn=20
        se = dsqrtf(nnn,sum/(dble(m)-1.d0))
      else
        se=0.d0
      endif
c
      do 2 i=1,m
      w(i)=1.0
    2 continue
      call ordst(m,qd1,qd3)
      call ordrs(m,i1,qb1,qb3)
      if(nit.lt.100) then
        return
      endif
c
      abd=1.d0/(qd1-qd3)
      aad=1.d0-abd*qd1
      abb=0.5d0/(qb1-qb3)
      aab=1.d0-abb*qb1
      do 5 i=1,m
      abi=dabs(b(i))
      wbi =1.0d0
      if(m.le.10) go to 3
      if(abi.le.qb1) go to 3
      wbi=aab+abb*abi
      if(wbi .lt.1.e-9) wbi =1.e-9
    3 continue
      wdi=1.0d0
      if(m.le.10) go to 4
      if(dist(i).le.qd1) go to 4
      wdi=aad+abd*dist(i)
      if(wdi.lt.0.1d0) wdi=0.1d0
    4 continue
      w(i)=wbi*wdi
      if(w(i).lt.1.e-9) w(i)=1.e-9
    5 continue
      wmax=w(1)
      do 6 i=2,m
      if(w(i).gt.wmax) wmax=w(i)
    6 continue
      do 7 i=1,m
      if(iwt(i).eq.iblnk) w(i)=w(i)/wmax
      if(iwt(i).eq.izero) w(i)=0.d0
      if(iwt(i).eq.ione) w(i)=1.d0
      if(iwt(i).eq.itwo) w(i)=0.5d0
      if(iwt(i).eq.ithree) w(i)=0.25d0
      if(iwt(i).eq.ifour) w(i)=0.1d0
    7 continue
      go to 15
    8 continue
      if(se.eq.0.d0) go to 13
      if(sem1.gt.se) sem1=se
      call ordst(m,qd1,qd3)
      call ordrs(m,i2,qb1,qb3)
      abd=1.d0/(qd1-qd3)
      aad=1.d0-abd*qd1
      abb=0.5d0/(qb1-qb3)
      aab=1.d0-abb*qb1
      do 11 i=1,m
      abi=dabs(bmax(i))
      wbi =1.0d0
      if(m.le.6) go to 9
      if(abi.le.qb1) go to 9
      wbi=aab+abb*abi
      if(wbi .lt.1.e-9) wbi =1.e-9
    9 continue
      wdi=1.0d0
      if(m.le.6) go to 10
      if(dist(i).le.qd1) go to 10
      wdi=aad+abd*dist(i)
      if(wdi.le.0.1d0) wdi=0.1d0
   10 continue
      w(i)=wdi*wbi
      if(w(i).le.1.e-9) w(i)=1.e-9
   11 continue
      wmax=w(1)
      do 12 i=2,m
      if(w(i).gt.wmax) wmax=w(i)
   12 continue
      if(iwt(i).eq.iblnk) w(i)=w(i)/wmax
      if(iwt(i).eq.izero) w(i)=0.d0
      if(iwt(i).eq.ione) w(i)=1.d0
      if(iwt(i).eq.itwo) w(i)=0.5d0
      if(iwt(i).eq.ithree) w(i)=0.25d0
      if(iwt(i).eq.ifour) w(i)=0.1d0
      go to 15
   13 continue
      do 14 i=1,m
      if(iwt(i).eq.izero) w(i)=0.d0
      if(iwt(i).eq.ione) w(i)=1.d0
      if(iwt(i).eq.itwo) w(i)=0.5d0
      if(iwt(i).eq.ithree) w(i)=0.25d0
      if(iwt(i).eq.ifour) w(i)=0.1d0
   14 continue
   15 continue
c
      return
      end

      subroutine tcal(ex,ey,ez,sx,sy,va,vb,vc,vd,et,ct,hemi,aoi)
c
c ...... calculate travel time
c
      implicit real*8 (a-h,o-z)
      call dazd(ex,ey,sx,sy,d,az)
      err=1.e-3
      z0=0.
      call ptcalc(va,vb,vc,vd,z0,ez,err,d,p,t,hemi,aoi,dtdz)
      ct=et+t
c
      return
      end

      subroutine infoco(m,n,u)
c
c ...... compute information matrix
c
c        indicates relative importance of the observed onset times
c
c        Reference: Uhrhammer, R. A., "The Optimal Estimation of
c        Earthquake Parameters", PEPI, Vol. 30, p. 105-118, 1982.
c
      implicit real*8 (a-h,o-z)
      common/files/ ncr,nsp,ns1,ns2,ns3,ns4
      common/info/sinfo(1000),sclinf
      real*8 u(1000,4)
      do 2 i=1,m
      sum=0.d0
      do 1 j=1,n
      sum=sum+u(i,j)*u(i,j)
    1 continue
      sinfo(i)=sum
    2 continue
      smax=sinfo(1)
      do 3 i=1,m
      if(sinfo(i).gt.smax) smax=sinfo(i)
    3 continue
      sclinf=smax
      if(sclinf.ne.0.d0) then
        do 4 i=1,m
        sinfo(i)=sinfo(i)/sclinf
    4   continue
      else
        do 5 i=1,m
        sinfo(i)=1.d0
    5   continue
      endif
c
      return
      end

      subroutine damper(a,b,w,m,n,mmax,q,u,v,qi,x,bmax,var,cn,
     1           meflag,nit)
c
c ...... damp equations of condition
c
      implicit real*8 (a-h,o-z)
      real*8 a(mmax,4),b(mmax),w(mmax),q(4),qi(4),u(mmax,4),
     1       v(4,4),x(4),bmax(mmax),ub(4)
      integer*4 meflag
      real*8 qils(4)
      data qmin,zero,one /0.00001d0,0.d0,1.d0/
c
      call nsfig(q,n)
      do 2 i=1,n
      sum=0.d0
      do 1 j=1,m
      sum=sum+u(j,i)*b(j)
    1 continue
      ub(i)=sum
    2 continue
c                       -1
c ...... limit size of q
c
      if(n.eq.3) then
        q_max=dmax1(q(1),q(2),q(3))
      else
        q_max=dmax1(q(1),q(2),q(3),q(4))
      endif
      if(nit.le.20) then
        q_min=qmin
      else
        q_min=qmin*1.35d0**dble(nit-20)
      endif
      do 3 i=1,n
      if(q(i).gt.0.d0) then
        qils(i)=one/q(i)
        if(qils(i).gt.100.d0) qils(i)=100.d0
      else
        qils(i)=zero
      endif
      if(q(i).eq.q_max) then
        qi(i)=one/q(i)
      else
        if(q(i).gt.q_min) then
          qi(i)=one/q(i)
        else
          qi(i)=zero
        endif
      endif
    3 continue
c                 -1    T
c ...... x = V * Q   * U  * b
c
      do 5 i=1,n
      sum=0.d0
      do 4 j=1,n
      sum=sum+v(i,j)*qi(j)*ub(j)
    4 continue
      x(i)=sum
    5 continue
      do 7 i=1,m
      sum=0.d0
      do 6 j=1,n
      sum=sum+a(i,j)*x(j)
    6 continue
      if(meflag.lt.1) then
        nnn=21
        bmax(i)=(b(i)-sum)/dsqrtf(nnn,w(i))
      else
        bmax(i)=b(i)
      endif
    7 continue
      sn=0.d0
      sd=0.d0
      do 8 i=1,m
      wi=w(i)
      sn=sn+wi*bmax(i)*bmax(i)
      sd=sd+wi
    8 continue
c
      do 9 i=1,n
      qi(i)=qils(i)
    9 continue
c
      var=sn/sd
      call psicom(m,n,u,qi,var)
c
      return
      end

      subroutine svdcom(a,m,n,mmax,nmax,q,u,v)
c
c ...... singular value decomposition procedure
c
c        Computes the singular value decomposition of a m by n matrix.
c        The procedure is described in detail in chapter 18 of
c        "Solving Least Squares Problems" by C. L. Lawson and
c        R. J. Hanson, Prentice-Hall, 1974.
c
c        PROFORT Compiler Version - 08 September 1987
c
      implicit real*8 (a-h,o-z)
      real*8 a(mmax,nmax),u(mmax,nmax),v(nmax,nmax),q(nmax),w(100)
      zero=0.0d0
      one=1.0d0
      tol=1.0e-290
      eps=1.0e-16
         do 1 i=1,m
            do 1 j=1,n
            u(i,j)=a(i,j)
    1    continue
      x=zero
      g=zero
         do  16 i=1,n
         w(i)=g
         s=zero
         l=i+1
            do 2 j=i,m
            s=s+u(j,i)**2
    2       continue
         if(s.ge.tol) go to 3
         g=zero
         go to 7
    3    f=u(i,i)
         nnn=22
         g=-dsqrtf(nnn,s)
         if(f.lt.zero) g=-g
         h=f*g-s
         hr=one/h
         u(i,i)=f-g
         if(l.gt.n) go to 7
            do 6 j=l,n
            s=zero
               do 4 k=i,m
               s=s+u(k,i)*u(k,j)
    4          continue
            f=s*hr
               do 5 k=i,m
               u(k,j)=u(k,j)+f*u(k,i)
    5          continue
    6       continue
    7    q(i)=g
         s=zero
         if(l.gt.n) go to 9
            do 8 j=l,n
            s=s+u(i,j)**2
    8       continue
    9    if(s.ge.tol) go to 10
         g=zero
         go to 15
   10    if(i.lt.n) f=u(i,i+1)
         nnn=23
         g=-dsqrtf(nnn,s)
         if(f.lt.zero) g=-g
         h=f*g-s
         hr=one/h
         if(i.lt.n) u(i,i+1)=f-g
         if(l.le.n) then
            do 11  j=l,n
            w(j)=u(i,j)*hr
   11       continue
         endif
         if(l.le.m) then
            do 14  j=l,m
            s=zero
            if(l.le.n) then
               do 12  k=l,n
               s=s+u(j,k)*u(i,k)
   12          continue
               do 13  k=l,n
               u(j,k)=u(j,k)+s*w(k)
   13          continue
            endif
   14       continue
         endif
   15       y=dabs(q(i))+dabs(w(i))
            if(y.gt.x) x=y
   16       continue
         do 108 i=1,n
         ii=n-i+1
         if(g.eq.zero) go to 105
         if(l.gt.n) go to 107
         h=u(ii,ii+1)*g
         hr=one/h
            do 101 j=l,n
            v(j,ii)=u(ii,j)*hr
  101       continue
            do 104 j=l,n
            s=zero
               do 102 k=l,n
               s=s+u(ii,k)*v(k,j)
  102          continue
               do 103 k=l,n
               v(k,j)=v(k,j)+s*v(k,ii)
  103          continue
  104       continue
  105    if(l.le.n) then
            do 106 j=l,n
            v(j,ii)=zero
            v(ii,j)=zero
  106       continue
         endif
  107    v(ii,ii)=one
         g=w(ii)
         l=ii
  108    continue
         do 211 i=1,n
         ii=n-i+1
         ll=ii+1
         g=q(ii)
         if(ll.le.n) then
            do 201 j=ll,n
            u(ii,j)=zero
  201       continue
         endif
         if(g.eq.zero) go to 208
         h=u(ii,ii)*g
         hr=one/h
         if(ll.le.n) then
            do 205 j=ll,n
            s=zero
               do 203 k=ll,m
               s=s+u(k,ii)*u(k,j)
  203          continue
            f=s*hr
               do 204 k=ii,m
               u(k,j)=u(k,j)+f*u(k,ii)
  204          continue
  205       continue
         endif
            gr=one/g
            do 207 j=ii,m
            u(j,ii)=u(j,ii)*gr
  207       continue
            go to 210
  208       do 209 j=ii,m
            u(j,ii)=zero
  209       continue
  210    u(ii,ii)=u(ii,ii)+one
  211    continue
      eps=eps*x
         do 265  k=1,n
         kk=n-k+1
  212       do 213  l=1,kk
            ll=kk-l+1
            if(dabs(w(ll)).le.eps) go to 230
            if(ll.eq.1) go to 230
            if(dabs(q(ll-1)).le.eps) go to 214
  213       continue
  214    c=zero
         s=one
         l1=ll-1
         if(kk.lt.ll) go to 230
            do 225  i=ll,kk
            f=s*w(i)
            w(i)=c*w(i)
            if(dabs(f) .le. eps) go to 230
            g=q(i)
            nnn=24
            q(i)=dsqrtf(nnn,f*f+g*g)
            h=q(i)
            hr=one/h
            c=g*hr
            s=-f*hr
               do 220  j=1,m
               y=u(j,l1)
               z=u(j,i)
               u(j,l1)=y*c+z*s
               u(j,i)=-y*s+z*c
  220          continue
  225       continue
  230    z=q(kk)
         if(ll.eq.kk) go to 255
         x=q(ll)
         if(kk.gt.1) y=q(kk-1)
         if(kk.gt.1) g=w(kk-1)
         h=w(kk)
         f=((y-z)*(y+z)+(g-h)*(g+h))/(2.*h*y)
         nnn=25
         g=dsqrtf(nnn,f*f+one)
         if(f.lt.zero) then
            f=((x-z)*(x+z)+h*(y/(f-g)-h))/x
         else
            f=((x-z)*(x+z)+h*(y/(f+g)-h))/x
         endif
         c=1
         s=1
         l2=ll+1
         if(kk.lt.l2) go to 250
            do 245  i=l2,kk
            g=w(i)
            y=q(i)
            h=s*g
            g=c*g
            nnn=26
            z=dsqrtf(nnn,f*f+h*h)
            zr=one/z
            w(i-1)=z
            c=f*zr
            s=h*zr
            f=x*c+g*s
            g=-x*s+g*c
            h=y*s
            y=y*c
               do 235  j=1,n
               x=v(j,i-1)
               z=v(j,i)
               v(j,i-1)=x*c+z*s
               v(j,i)=-x*s+z*c
  235          continue
            nnn=27
            z=dsqrtf(nnn,f*f+h*h)
            zr=one/z
            q(i-1)=z
            c=f*zr
            s=h*zr
            f=c*g+s*y
            x=-s*g+c*y
               do 240  j=1,m
               y=u(j,i-1)
               z=u(j,i)
               u(j,i-1)=y*c+z*s
               u(j,i)=-y*s+z*c
  240          continue
  245       continue
  250    w(ll)=zero
         w(kk)=f
         q(kk)=x
         go to 212
  255    if(z.lt.zero) then
            q(kk)=-z
            do 260  j=1,n
            v(j,kk)=-v(j,kk)
  260       continue
         endif
  265    continue
c
      return
      end

      subroutine corr(cov,n,cor,nit)
c
c ...... calculate correlation matrix for the hypocental parameters
c
      implicit real*8 (a-h,o-z)
      real*8 cov(4,4),cor(4,4)
      do 2 i=1,4
      do 1 j=1,4
      cor(i,j)=0.d0
    1 continue
    2 continue
      do 4 i=1,n
      do 3 j=1,n
      arg=cov(i,i)*cov(j,j)
      if(arg.gt.0.d0) then
        nnn=28
        cor(i,j)=cov(i,j)/dsqrtf(nnn,arg)
      else
        if(cov(i,j).lt.0.d0) then
          cor(i,j)=-1.0d0
        elseif(cov(i,j).gt.0.d0) then
          cor(i,j)=+1.0d0
        else
          if(i.eq.j) then
            cor(i,j)=1.d0
          else
            cor(i,j)=0.d0
          endif
        endif
      endif
    3 continue
    4 continue
c
      return
      end

      subroutine psicom(m,n,u,qi,var)
c
c ...... calculate information density matrix
c
c        For details see above Uhrhammer reference.
c
      implicit real*8 (a-h,o-z)
      common/psiv/psi(1000),psiscl
      real*8 u(1000,4),qi(4)
      do 2 i=1,m
      sum=0.d0
      do 1 j=1,n
      sum=sum+u(i,j)*u(i,j)*qi(j)*qi(j)*var
    1 continue
      if(sum.gt.0.d0) then
        nnn=29
        psi(i)=dsqrtf(nnn,sum)
      else
        psi(i)=0.d0
      endif
    2 continue
      psimax=psi(1)
      do 3 i=1,m
      if(psi(i).gt.psimax) psimax=psi(i)
    3 continue
      psiscl=psimax
      do 4 i=1,m
      if(psimax.ne.0.d0) then
        psi(i)=psi(i)/psimax
      else
        psi(i)=0.d0
      endif
    4 continue
c
      return
      end

      subroutine ptdcal(a,b,c,d,z0,zh,err,x,p,t,dtda,dtdb,dtdc,dtdd,
     1                  dtdz,dtdx,dtdt,hemi,aoi)
c
c ...... calculate travel time for a gradient layer over a halfspace
c        velocity model
c
c        02 SEPTEMBER 1987
c
c        Modified 16 March 1994 to work with hypocenters in halfspace
c
c        input -
c           velocity model parameters - a, b, c & d
c           distance in km            - x
c           source and station depths - zh & z0
c
      implicit real*8 (a-h,o-z)
      dtda=0.d0
      dtdb=0.d0
      dtdc=0.d0
      dtdd=0.d0
      call ptcalc(a,b,c,d,z0,zh,err,x,p,t,hemi,aoi,dtdz)
      dtdx=p
      dtdt=1.d0
c
      return
      end

      subroutine ptcalc(a,b,c,d,z0,zh,err,x,p,t,hemi,aoi,dtdz)
c
c ...... calculate slowness and travel time
c
      implicit real*8 (a-h,o-z)
      arg=45.d0/datan(1.d0)
      dtdz=0.d0
      im1=-1
      ip1=+1
      if(zh.gt.d) then
        p=pcalc(a,b,c,d,z0,zh,x,err,pc1,pc2,pc3,xc3,im1)
        t=tcalc(a,b,c,d,z0,zh,p,pc1,pc2,pc3,xc3,x,im1)
        hemi=-1.d0
        q=p*c
        if(q.lt.1.d0) then
          aoi=arg*dasin(q)
        else
          aoi=90.d0
        endif
        pc3=1.0d0/c
        temp=pc3**2-p**2
        if(temp.gt.0.d0) then
          nnn=30
          dtdz=dsqrtf(nnn,temp)
        else
          dtdz=0.d0
        endif
        return
      endif
      pc1=1.d0/(a+b*zh)
      pc2=1.d0/(a+b*d)
      pc3=1.d0/c
      xc1=xcalc(a,b,c,d,z0,zh,pc1,pc1,pc2,pc3,im1)
      xc2=xcalc(a,b,c,d,z0,zh,pc2,pc1,pc2,pc3,ip1)
      xc3=xcalc(a,b,c,d,z0,zh,pc3,pc1,pc2,pc3,ip1)
      if(x.le.xc1) then
        p=pcalc(a,b,c,d,z0,zh,x,err,pc1,pc2,pc3,xc3,im1)
        t=tcalc(a,b,c,d,z0,zh,p,pc1,pc2,pc3,xc3,x,im1)
        hemi=-1.d0
        q=p*(a+b*zh)
        if(q.lt.1.d0) then
          aoi=arg*dasin(q)
        else
          aoi=90.d0
        endif
        temp=pc1**2-p**2
        if(temp.gt.0.d0) then
          nnn=31
          dtdz=dsqrtf(nnn,temp)
        else
          dtdz=0.d0
        endif
      elseif(x.le.xc2) then
        p=pcalc(a,b,c,d,z0,zh,x,err,pc1,pc2,pc3,xc3,ip1)
        t=tcalc(a,b,c,d,z0,zh,p,pc1,pc2,pc3,xc3,x,ip1)
        hemi=+1.d0
        q=p*(a+b*zh)
        if(q.lt.1.d0) then
          aoi=arg*dasin(q)
        else
          aoi=90.d0
        endif
        temp=pc1**2-p**2
        if(temp.gt.0.d0) then
          nnn=32
          dtdz=-1.0d0*dsqrtf(nnn,temp)
        else
          dtdz=0.d0
        endif
        if(x.lt.xc3) then
          return
        endif
        t1=tcalc(a,b,c,d,z0,zh,pc3,pc1,pc2,pc3,xc3,x,ip1)
        if(t1.gt.t) then
          return
        endif
        p=pc3
        t=t1
        q=p*(a+b*zh)
        if(q.lt.1.d0) then
          aoi=arg*dasin(q)
        else
          aoi=90.d0
        endif
        temp=pc1**2-p**2
        if(temp.gt.0.d0) then
          nnn=33
          dtdz=-1.0d0*dsqrtf(nnn,temp)
        else
          dtdz=0.d0
        endif
      else
        p=pc3
        t=tcalc(a,b,c,d,z0,zh,pc3,pc1,pc2,pc3,xc3,x,ip1)
        hemi=+1.d0
        q=p*(a+b*zh)
        if(q.lt.1.d0) then
          aoi=arg*dasin(q)
        else
          aoi=90.d0
        endif
        temp=pc1**2-p**2
        if(temp.gt.0.d0) then
          nnn=34
          dtdz=-1.0d0*dsqrtf(nnn,temp)
        else
          dtdz=0.d0
        endif
      endif
c
      return
      end

      real*8 function xcalc(a,b,c,d,z0,zh,p,pc1,pc2,pc3,jd)
c
c ...... calculate distance for a given slowness, model & path
c
      implicit real*8 (a-h,o-z)
c
c      write(*,*) 'xcalc - a,b,c,d = ',a,b,c,d
c      write(*,*) 'xcalc - z0,zh,p = ',z0,zh,p
c
      q=-1.d0/(b*p)
      if(jd.eq.+1) then
        if(p.ge.pc1) then
          xcalc=q*(r5(a,b,p,zh)-r5(a,b,p,z0))
        elseif(p.ge.pc2) then
          xcalc=-1.d0*q*(r5(a,b,p,z0)+r5(a,b,p,zh))
        elseif(p.ge.pc3) then
          xcalc=q*(2.*r5(a,b,p,d)-r5(a,b,p,z0)-r5(a,b,p,zh))
        else
          nnn=36
          xcalc=q*(r5(a,b,p,d)-r5(a,b,p,z0))+p*c*(zh-d)/
     1          dsqrtf(nnn,1.d0-(p*c)**2)
        endif
      elseif(zh.gt.d) then
        nnn=37
        xcalc=q*(r5(a,b,p,d)-r5(a,b,p,z0))
     1        +p*c*(zh-d)/dsqrtf(nnn,1.d0-(p*c)**2)
      else
        xcalc=q*(r5(a,b,p,zh)-r5(a,b,p,z0))
      endif
c
      return
      end

      real*8 function r5(a,b,p,z)
c
c ...... function called by xcalc
c
      implicit real*8 (a-h,o-z)
      arg=1.d0-(p*(a+b*z))**2
      if(arg.le.0.d0) then
        r5=0.d0
      else
        nnn=36
        r5=dsqrtf(nnn,arg)
      endif
c
      return
      end

      real*8 function tcalc(a,b,c,d,z0,zh,p,pc1,pc2,pc3,xc3,x,jd)
c
c ...... calculate travel time for a given slowness, model and path
c
      implicit real*8 (a-h,o-z)
      q=-1.d0/b
      t0=1.d0/(z0+a/b)
      th=1.d0/(zh+a/b)
      td=1.d0/(d+a/b)
      if(jd.eq.+1) then
        if(p.ge.pc1) then
          tcalc=q*(r4(b,p,th)-r4(b,p,t0))
        elseif(p.ge.pc2) then
          tp=b*p
          tcalc=q*(2.*r4(b,p,tp)-r4(b,p,t0)-r4(b,p,th))
        elseif(p.gt.pc3) then
          tcalc=q*(2.*r4(b,p,td)-r4(b,p,t0)-r4(b,p,th))
        elseif(p.eq.pc3) then
          tcalc=q*(2.*r4(b,p,td)-r4(b,p,t0)-r4(b,p,th))+(x-xc3)/c
        else
          nnn=37
          tcalc=q*(r4(b,p,td)-r4(b,p,t0)-r4(b,p,th))+
     1          p*(zh-d)/dsqrtf(nnn,1.-(p*c)**2)
        endif
      elseif(zh.gt.d) then
        nnn=38
        tcalc=q*(r4(b,p,td)-r4(b,p,t0))+p*(zh-d)/
     1        dsqrtf(nnn,1.-(p*c)**2)
      else
        tcalc=q*(r4(b,p,th)-r4(b,p,t0))
      endif
c
      return
      end

      real*8 function r4(b,p,t)
c
c ...... function called by tcalc
c
      implicit real*8 (a-h,o-z)
      arg=t**2-(b*p)**2
      if(arg.le.0.d0) then
        sarg=0.d0
      else
        nnn=39
        sarg=dsqrtf(nnn,arg)
      endif
      r4=dlog(2.d0*(sarg+t))
c
      return
      end

      real*8 function pcalc(a,b,c,d,z0,zh,x,err,pc1,pc2,pc3,xc3,jd)
c
c ...... calculate slowness
c
      implicit real*8 (a-h,o-z)
c
      if(zh.le.d) then
        if(x.eq.0.0d0) then
          pcalc=0.0d0
        else
          nnn=40
          pcalc=2.d0*x*b/dsqrtf(nnn,((x*b)**2+(a+b*zh)**2+
     1          (a+b*z0)**2)**2-4.d0*((a+b*zh)**2)*((a+b*z0)**2))
        endif
      else
        if(x.eq.0.0d0) then
          pcalc=0.0d0
        else
          call zhgtdpf(a,b,c,d,z0,zh,x,err,p)
          pcalc=p
c
c          nnn=41
c          ps=x/(dsqrtf(nnn,x**2+(zh-d)**2)*c)
c          xs=xcalc(a,b,c,d,z0,zh,ps,pc1,pc2,pc3,jd)
c          ps=dsin(x/zh)/c
c          xs=xcalc(a,b,c,d,z0,zh,ps,pc1,pc2,pc3,jd)
c          p=ps
c          dx=x-xs
c          ij=0
c    1     dpdx=dpdxc(a,b,c,d,z0,zh,p,pc1,pc2,pc3,jd)
c          if(dabs(dx).gt.1.d0) then
c            p=p+dpdx*dsign(1.d0,dx)
c          else
c            p=p+dpdx*dx
c          endif
c          write(*,*) 'x, zh, p, dpdx & dx = ',x,zh,p,dpdx,dx
c          if(p.ge.(1.d0/c)) then
c            write(*,*) 'pcalc - p >= 1/c - ',p,1.d0/c
c            write(*,*) 'dpdx & dx = ',dpdx,dx
c            stop
c          endif
c          xc=xcalc(a,b,c,d,z0,zh,p,pc1,pc2,pc3,jd)
c          dx=x-xc
c          ij=ij+1
c          if(dabs(dx).gt.err) go to 1
c          pcalc=p
        endif
      endif
c
      return
      end

      subroutine zhgtdpf(a,b,c,d,z0,zh,x,err,p)
c
c ...... numerically search for the slowness when the zh > d
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 pt(1000),xt(1000)
      data ncall /0/
c
c	write(*,*) 'zhgtdpf - a,b,c,d,z0,zh,x,err'
c	write(*,'(a,4f10.4)') 'a,b,c,d     = ',a,b,c,d
c	write(*,'(a,4f10.6)') 'z0,zh,x,err = ',z0,zh,x,err
      if(ncall.eq.0) then
        call pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
        at=a
        bt=b
        ct=c
        dt=d
        z0t=z0
        zht=zh
        xtt=x
        ncall=1
      endif
      if(    a.ne.at) then
        call pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
      elseif(b.ne.bt) then
        call pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
      elseif(c.ne.ct) then
        call pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
      elseif(d.ne.dt) then
        call pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
      elseif(z0.ne.z0t) then
        call pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
      elseif(zh.ne.zht) then
        call pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
      elseif(x.ne.xtt) then
        call pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
      endif
      at=a
      bt=b
      ct=c
      dt=d
      z0t=z0
      zht=zh
      xtt=x
c
c      write(*,*) 'zhgtdpf - x = ',x
      i=0
    1 i=i+1
      if(x.gt.xt(i)) go to 1
      dpdx=(pt(i)-pt(i-1))/(xt(i)-xt(i-1))
      pp=pt(i-1)+dpdx*(x-xt(i-1))
      nnn=53
      q=-1.d0/(b*pp)
      xx=q*(r5(a,b,pp,d)-r5(a,b,pp,z0))+
     1   pp*c*(zh-d)/dsqrtf(nnn,1.d0-(pp*c)**2)
c
      p0=pt(i-1)
      p1=pt(i)
      x0=xt(i-1)
      x1=xt(i)
	loop=0
    2 continue
	loop=loop+1
c	write(*,'(a,i10)') 'zhgtdpf - loop = ',loop
c      write(*,*) 'zhgtdpf - p0, pp & p1 = ',p0,pp,p1
c      write(*,*) 'zhgtdpf - x0, xx & x1 = ',x0,xx,x1
      dx=x-xx
c      write(*,*) 'zhgtdpf - x, dx & err = ',x,dx,err
      if(dabs(x-xx).lt.err) then
        p=pp
        return
      endif
      if(loop.gt.100) then
	p=pp
	return
      endif
      if(x.lt.xx) then
        p1=pp
        x1=xx
        dpdx=(p1-p0)/(x1-x0)
        pp=p0+dpdx*(x-x0)
        q=-1.d0/(b*pp)
        nnn=54
        xx=q*(r5(a,b,pp,d)-r5(a,b,pp,z0))+
     1     pp*c*(zh-d)/dsqrtf(nnn,1.d0-(pp*c)**2)
      else
        p0=pp
        x0=xx
        dpdx=(p1-p0)/(x1-x0)
        pp=p0+dpdx*(x-x0)
        q=-1.d0/(b*pp)
        nnn=55
        xx=q*(r5(a,b,pp,d)-r5(a,b,pp,z0))+
     1     pp*c*(zh-d)/dsqrtf(nnn,1.d0-(pp*c)**2)
      endif
      go to 2
c
      end

      subroutine pxfan(a,b,c,d,z0,zh,x,pt,xt,nt)
c
c ...... generate slowness and distance fan
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 pt(1000),xt(1000)
c
c      write(*,*) 'pxfan - a, b, c & d = ',a,b,c,d
c      write(*,*) 'pxfan - z0 & zh = ',z0,zh
      nnn=52
      nt=20
      p0=(0.0001d0/dsqrt(1.d-8+zh*zh))/c
      p1=(x/dsqrt(x*x+(zh-d)*(zh-d)))/c
c      write(*,*) 'pxfan - p0 & p1 = ',p0,p1
      con=2.d0*datan(1.d0)/dble(nt-1)
      do 1 i=1,nt
      p=p0+(p1-p0)*dsin(con*dble(i-1))
      pt(i)=p
      q=-1.d0/(b*p)
      xt(i)=q*(r5(a,b,p,d)-r5(a,b,p,z0))+
     1      p*c*(zh-d)/dsqrtf(nnn,1.d0-(p*c)**2)
    1 continue
c      do 2 i=1,nt
c      write(*,*) 'pxfan - i, pt(i) & xt(i) = ',i,' ',pt(i),xt(i)
c    2 continue
c
      return
      end

      real*8 function dpdxc(a,b,c,d,z0,zh,p,pc1,pc2,pc3,jd)
c
c ...... calculate dp/dx for ray shooting in halfspace
c
      implicit real*8 (a-h,o-z)
      if(zh.gt.d) then
        dxdp=r1(a,b,p,d)-r1(a,b,p,z0)+r2(a,b,p,d)-r2(a,b,p,z0)
     1      +r3(c,p,zh)-r3(c,p,d)+r6(c,p,zh)-r6(c,p,d)
        dpdxc=1.d0/dxdp
      elseif(jd.eq.-1) then
        dxdp=r1(a,b,p,zh)-r1(a,b,p,z0)+r2(a,b,p,zh)-r2(a,b,p,z0)
        dpdxc=1.d0/dxdp
      elseif(p.lt.pc1) then
        if(p.ge.pc2) then
          zp=(1.d0-a*p)/(b*p)
          dxdp=2.d0*(r1(a,b,p,zp)+r2(a,b,p,zp))-r1(a,b,p,z0)
     1         -r1(a,b,p,zh)-r2(a,b,p,z0)-r2(a,b,p,zh)
          dpdxc=1.d0/dxdp
        elseif(p.gt.pc3) then
          dxdp=2.d0*(r1(a,b,p,d)+r2(a,b,p,d))-r1(a,b,p,z0)
     1         -r1(a,b,p,zh)-r2(a,b,p,z0)-r2(a,b,p,zh)
          dpdxc=1./dxdp
        elseif(p.lt.pc3) then
          dxdp=r1(a,b,p,d)-r1(a,b,p,z0)+r2(a,b,p,d)-r2(a,b,p,z0)
     1        +r3(c,p,zh)-r3(c,p,d)+r6(c,p,zh)-r6(c,p,d)
          dpdxc=1.d0/dxdp
        else
          dpdxc=0.d0
        endif
      else
        dxdp=r1(a,b,p,zh)-r1(a,b,p,z0)+r2(a,b,p,zh)-r2(a,b,p,z0)
        dpdxc=1.d0/dxdp
      endif
c
      return
      end

      real*8 function r1(a,b,p,z)
c
c ...... function called by dxdp
c
      implicit real*8 (a-h,o-z)
      arg=1.d0-(p*(a+b*z))**2
      if(arg.le.0.d0) then
        r1=0.d0
      else
        nnn=42
        r1=dsqrtf(nnn,arg)/(b*p*p)
      endif
      return
    1 r1=0.d0
c
      return
      end

      real*8 function r2(a,b,p,z)
c
c ...... function called by dxdp
c
      implicit real*8 (a-h,o-z)
      arg=1.d0-(p*(a+b*z))**2
      if(arg.le.1.e-6) then
        r2=0.d0
      else
        nnn=43
        r2=((a+b*z)**2)/(b*dsqrtf(nnn,arg))
      endif
c
      return
      end

      real*8 function r3(c,p,z)
c
c ...... function called by dxdp
c
      implicit real*8 (a-h,o-z)
      arg=1.d0-(p*c)**2
      nnn=44
      r3=c*z/dsqrtf(nnn,arg)
c
      return
      end

      real*8 function r6(c,p,z)
c
c ...... function called by dxdp
c
      implicit real*8 (a-h,o-z)
      nnn=45
      sarg=dsqrtf(nnn,1.d0-(p*c)**2)
      r6=(c*z*(c*p)**2)/sarg**3
c
      return
      end
