surfdisp96.f90

subroutine surfdisp96(thkm, vpm, vsm, rhom, nlayer, iflsph, iwave, mode, igr, kmax, t, cg, err)
parameter(ler = 0, lin = 5, lot = 6)
integer :: nl, nl2, nlay
parameter(nl = 100, nlay = 100, nl2 = nl + nl)
integer :: np
parameter(np = 60)
real(4) :: thkm(nlay), vpm(nlay), vsm(nlay), rhom(nlay)
integer :: nlayer, iflsph, iwave, mode, igr, kmax, err
double precision :: twopi, one, onea
double precision :: cc, c1, clow, cm, dc, t1
double precision :: t(np), c(np), cb(np), cg(np)
real(4) :: d(nl), a(nl), b(nl), rho(nl), rtp(nl), dtp(nl), btp(nl)
integer(4) :: iverb(2)
integer(4) :: llw
integer(4) :: nsph, ifunc, idispl, idispr, is, ie
real(4) :: sone0, ddc0, h0, sone, ddc, h
mmax = nlayer
nsph = iflsph
err = 0
do i = 1, mmax
    b(i) = vsm(i)
    a(i) = vpm(i)
    d(i) = thkm(i)
    rho(i) = rhom(i)
    39 continue
end do
if (iwave == 1) then
    idispl = kmax
    idispr = 0
else
    if (iwave == 2) then
        idispl = 0
        idispr = kmax
    end if
end if
iverb(1) = 0
iverb(2) = 0
sone0 = 1.500
ddc0 = 0.005
h0 = 0.005
llw = 1
if (b(1) <= 0.0) then
    llw = 2
end if
twopi = 2.d0*3.141592653589793d0
one = 1.0d-2
if (nsph == 1) then
    call sphere(0, 0, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
end if
jmn = 1
betmx = -1.e20
betmn = 1.e20
do i = 1, mmax
    if (b(i) > 0.01 .and. b(i) < betmn) then
        betmn = b(i)
        jmn = i
        jsol = 1
    else
        if (b(i) <= 0.01 .and. a(i) < betmn) then
            betmn = a(i)
            jmn = i
            jsol = 0
        end if
    end if
    if (b(i) > betmx) then
        betmx = b(i)
    end if
    20 continue
end do
do ifunc = 1, 2
    if (ifunc == 1 .and. idispl <= 0) then
        go to 2000
    end if
    if (ifunc == 2 .and. idispr <= 0) then
        go to 2000
    end if
    if (nsph == 1) then
        call sphere(ifunc, 1, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
    end if
    ddc = ddc0
    sone = sone0
    h = h0
    if (sone < 0.01) then
        sone = 2.0
    end if
    onea = dble(sone)
    if (jsol == 0) then
        cc1 = betmn
    else
        call gtsolh(a(jmn), b(jmn), cc1)
    end if
    cc1 = .95*cc1
    cc1 = .90*cc1
    cc = dble(cc1)
    dc = dble(ddc)
    dc = dabs(dc)
    c1 = cc
    cm = cc
    do i = 1, kmax
        cb(i) = 0.0d0
        c(i) = 0.0d0
        450 continue
    end do
    ift = 999
    do iq = 1, mode
        is = 1
        ie = kmax
        itst = ifunc
        do k = is, ie
            if (k >= ift) then
                go to 1700
            end if
            t1 = dble(t(k))
            if (igr > 0) then
                t1a = t1/(1. + h)
                t1b = t1/(1. - h)
                t1 = dble(t1a)
            else
                t1a = sngl(t1)
                tlb = 0.0
            end if
            if (k == is .and. iq == 1) then
                c1 = cc
                clow = cc
                ifirst = 1
            else
                if (k == is .and. iq > 1) then
                    c1 = c(is) + one*dc
                    clow = c1
                    ifirst = 1
                else
                    if (k > is .and. iq > 1) then
                        ifirst = 0
                        clow = c(k) + one*dc
                        c1 = c(k - 1)
                        if (c1 < clow) then
                            c1 = clow
                        end if
                    else
                        if (k > is .and. iq == 1) then
                            ifirst = 0
                            c1 = c(k - 1) - onea*dc
                            clow = cm
                        end if
                    end if
                end if
            end if
            call getsol(t1, c1, clow, dc, cm, betmx, iret, ifunc, ifirst, d, a, b, rho, rtp, dtp, btp, mmax, llw)
            if (iret == -1) then
                go to 1700
            end if
            c(k) = c1
            if (igr > 0) then
                t1 = dble(t1b)
                ifirst = 0
                clow = cb(k) + one*dc
                c1 = c1 - onea*dc
                call getsol(t1, c1, clow, dc, cm, betmx, iret, ifunc, ifirst, d, a, b, rho, rtp, dtp, btp, mmax, llw)
                if (iret == -1) then
                    c1 = c(k)
                end if
                cb(k) = c1
            else
                c1 = 0.0d+00
            end if
            cc0 = sngl(c(k))
            cc1 = sngl(c1)
            if (igr == 0) then
                cg(k) = cc0
            else
                gvel = (1/t1a - 1/t1b)/(1/t1a*cc0 - 1/t1b*cc1)
                cg(k) = gvel
            end if
            1600 continue
        end do
        go to 1800
        1700 if (iq > 1) then
            go to 1750
        end if
        if (iverb(ifunc) == 0) then
            iverb(ifunc) = 1
            err = 1
        end if
        1750 ift = k
        itst = 0
        do i = k, ie
            t1a = t(i)
            cg(i) = 0.0
            1770 continue
        end do
        1800 continue
    end do
    2000 continue
end do
end subroutine surfdisp96

subroutine gtsolh(a, b, c)
real(4) :: kappa, k2, gk2
c = 0.95*b
do i = 1, 5
    gamma = b/a
    kappa = c/b
    k2 = kappa**2
    gk2 = (gamma*kappa)**2
    fac1 = sqrt(1.0 - gk2)
    fac2 = sqrt(1.0 - k2)
    fr = (2.0 - k2)**2 - 4.0*fac1*fac2
    frp = (-4.0)*(2.0 - k2)*kappa + 4.0*fac2*gamma*gamma*kappa/fac1 + 4.0*fac1*kappa/fac2
    frp = frp/b
    c = c - fr/frp
    100 continue
end do
return
end subroutine gtsolh

subroutine getsol(t1, c1, clow, dc, cm, betmx, iret, ifunc, ifirst, d, a, b, rho, rtp, dtp, btp, mmax, llw)
parameter(nl = 100)
real(8) :: wvno, omega, twopi
real(8) :: c1, c2, cn, cm, dc, t1, clow
real(8) :: dltar, del1, del2, del1st, plmn
save :: del1st
real(4) :: d(nl), a(nl), b(nl), rho(nl), rtp(nl), dtp(nl), btp(nl)
integer :: llw, mmax
twopi = 2.d0*3.141592653589793d0
omega = twopi/t1
wvno = omega/c1
del1 = dltar(wvno, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
if (ifirst == 1) then
    del1st = del1
end if
plmn = dsign(1.0d+00, del1st)*dsign(1.0d+00, del1)
if (ifirst == 1) then
    idir = 1
else
    if (ifirst /= 1 .and. plmn >= 0.0d+00) then
        idir = 1
    else
        if (ifirst /= 1 .and. plmn < 0.0d+00) then
            idir = -1
        end if
    end if
end if
1000 continue
if (idir > 0) then
    c2 = c1 + dc
else
    c2 = c1 - dc
end if
if (c2 <= clow) then
    idir = 1
    c1 = clow
end if
if (c2 <= clow) then
    go to 1000
end if
omega = twopi/t1
wvno = omega/c2
del2 = dltar(wvno, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
if (dsign(1.0d+00, del1) /= dsign(1.0d+00, del2)) then
    go to 1300
end if
c1 = c2
del1 = del2
if (c1 < cm) then
    go to 1700
end if
if (c1 >= betmx + dc) then
    go to 1700
end if
go to 1000
1300 call nevill(t1, c1, c2, del1, del2, ifunc, cn, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
c1 = cn
if (c1 > (betmx)) then
    go to 1700
end if
iret = 1
return
1700 continue
iret = -1
return
end subroutine getsol

subroutine sphere(ifunc, iflag, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
parameter(nl = 100, np = 60)
real(4) :: d(nl), a(nl), b(nl), rho(nl), rtp(nl), dtp(nl), btp(nl)
integer :: mmax, llw
double precision :: z0, z1, r0, r1, dr, ar, tmp, twopi
save :: dhalf
ar = 6370.0d0
dr = 0.0d0
r0 = ar
d(mmax) = 1.0
if (iflag == 0) then
    do i = 1, mmax
        dtp(i) = d(i)
        rtp(i) = rho(i)
        5 continue
    end do
    do i = 1, mmax
        dr = dr + dble(d(i))
        r1 = ar - dr
        z0 = ar*dlog(ar/r0)
        z1 = ar*dlog(ar/r1)
        d(i) = z1 - z0
        tmp = (ar + ar)/(r0 + r1)
        a(i) = a(i)*tmp
        b(i) = b(i)*tmp
        btp(i) = tmp
        r0 = r1
        10 continue
    end do
    dhalf = d(mmax)
else
    d(mmax) = dhalf
    do i = 1, mmax
        if (ifunc == 1) then
            rho(i) = rtp(i)*btp(i)**(-5)
        else
            if (ifunc == 2) then
                rho(i) = rtp(i)*btp(i)**(-2.275)
            end if
        end if
        30 continue
    end do
end if
d(mmax) = 0.0
return
end subroutine sphere

subroutine nevill(t, c1, c2, del1, del2, ifunc, cc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
implicit double precision (a-h,o-z)
parameter(nl = 100, np = 60)
real(4) :: d(nl), a(nl), b(nl), rho(nl), rtp(nl), dtp(nl), btp(nl)
dimension :: x(20), y(20)
integer :: llw, mmax
omega = twopi/t
call half(c1, c2, c3, del3, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
nev = 1
nctrl = 1
100 continue
nctrl = nctrl + 1
if (nctrl >= 100) then
    go to 1000
end if
if (c3 < dmin1(c1, c2) .or. c3 > dmax1(c1, c2)) then
    nev = 0
    call half(c1, c2, c3, del3, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
end if
s13 = del1 - del3
s32 = del3 - del2
if (dsign(1.d+00, del3)*dsign(1.d+00, del1) < 0.0d+00) then
    c2 = c3
    del2 = del3
else
    c1 = c3
    del1 = del3
end if
if (dabs(c1 - c2) <= 1.d-6*c1) then
    go to 1000
end if
if (dsign(1.0d+00, s13) /= dsign(1.0d+00, s32)) then
    nev = 0
end if
ss1 = dabs(del1)
s1 = 0.01*ss1
ss2 = dabs(del2)
s2 = 0.01*ss2
if (s1 > ss2 .or. s2 > ss1 .or. nev == 0) then
    call half(c1, c2, c3, del3, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
    nev = 1
    m = 1
else
    if (nev == 2) then
        x(m + 1) = c3
        y(m + 1) = del3
    else
        x(1) = c1
        y(1) = del1
        x(2) = c2
        y(2) = del2
        m = 1
    end if
    do kk = 1, m
        j = m - kk + 1
        denom = y(m + 1) - y(j)
        if (dabs(denom) < 1.0d-10*abs(y(m + 1))) then
            go to 950
        end if
        x(j) = ((-y(j))*x(j + 1) + y(m + 1)*x(j))/denom
        900 continue
    end do
    c3 = x(1)
    wvno = omega/c3
    del3 = dltar(wvno, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
    nev = 2
    m = m + 1
    if (m > 10) then
        m = 10
    end if
    go to 951
    950 continue
    call half(c1, c2, c3, del3, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
    nev = 1
    m = 1
    951 continue
end if
go to 100
1000 continue
cc = c3
return
end subroutine nevill

subroutine half(c1, c2, c3, del3, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
implicit double precision (a-h,o-z)
parameter(nl = 100)
real(4) :: d(nl), a(nl), b(nl), rho(nl), rtp(nl), dtp(nl), btp(nl)
c3 = 0.5*(c1 + c2)
wvno = omega/c3
del3 = dltar(wvno, omega, ifunc, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
return
end subroutine half

function dltar(wvno, omega, kk, d, a, b, rho, rtp, dtp, btp, mmax, llw, twop)
implicit double precision (a-h,o-z)
parameter(nl = 100)
real(4) :: d(nl), a(nl), b(nl), rho(nl), rtp(nl), dtp(nl), btp(nl)
if (kk == 1) then
    dltar = dltar1(wvno, omega, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
else
    if (kk == 2) then
        dltar = dltar4(wvno, omega, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
    end if
end if
end function dltar

function dltar1(wvno, omega, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
implicit double precision (a-h,o-z)
parameter(nl = 100, np = 60)
real(4) :: d(nl), a(nl), b(nl), rho(nl), rtp(nl), dtp(nl), btp(nl)
integer :: llw, mmax
beta1 = dble(b(mmax))
rho1 = dble(rho(mmax))
xkb = omega/beta1
wvnop = wvno + xkb
wvnom = dabs(wvno - xkb)
rb = dsqrt(wvnop*wvnom)
e1 = rho1*rb
e2 = 1.d+00/beta1*beta1
mmm1 = mmax - 1
do m = mmm1, llw, -1
    beta1 = dble(b(m))
    rho1 = dble(rho(m))
    xmu = rho1*beta1*beta1
    xkb = omega/beta1
    wvnop = wvno + xkb
    wvnom = dabs(wvno - xkb)
    rb = dsqrt(wvnop*wvnom)
    q = dble(d(m))*rb
    if (wvno < xkb) then
        sinq = dsin(q)
        y = sinq/rb
        z = (-rb)*sinq
        cosq = dcos(q)
    else
        if (wvno == xkb) then
            cosq = 1.0d+00
            y = dble(d(m))
            z = 0.0d+00
        else
            fac = 0.0d+00
            if (q < 16) then
                fac = dexp((-2.0d+0)*q)
            end if
            cosq = (1.0d+00 + fac)*0.5d+00
            sinq = (1.0d+00 - fac)*0.5d+00
            y = sinq/rb
            z = rb*sinq
        end if
    end if
    e10 = e1*cosq + e2*xmu*z
    e20 = e1*y/xmu + e2*cosq
    xnor = dabs(e10)
    ynor = dabs(e20)
    if (ynor > xnor) then
        xnor = ynor
    end if
    if (xnor < 1.d-40) then
        xnor = 1.0d+00
    end if
    e1 = e10/xnor
    e2 = e20/xnor
    600 continue
end do
dltar1 = e1
return
end function dltar1

function dltar4(wvno, omga, d, a, b, rho, rtp, dtp, btp, mmax, llw, twopi)
implicit double precision (a-h,o-z)
parameter(nl = 100, np = 60)
dimension :: e(5), ee(5), ca(5,5)
real(4) :: d(nl), a(nl), b(nl), rho(nl), rtp(nl), dtp(nl), btp(nl)
omega = omga
if (omega < 1.0d-4) then
    omega = 1.0d-4
end if
wvno2 = wvno*wvno
xka = omega/dble(a(mmax))
xkb = omega/dble(b(mmax))
wvnop = wvno + xka
wvnom = dabs(wvno - xka)
ra = dsqrt(wvnop*wvnom)
wvnop = wvno + xkb
wvnom = dabs(wvno - xkb)
rb = dsqrt(wvnop*wvnom)
t = dble(b(mmax))/omega
gammk = 2.d+00*t*t
gam = gammk*wvno2
gamm1 = gam - 1.d+00
rho1 = dble(rho(mmax))
e(1) = rho1*rho1*(gamm1*gamm1 - gam*gammk*ra*rb)
e(2) = (-rho1)*ra
e(3) = rho1*(gamm1 - gammk*ra*rb)
e(4) = rho1*rb
e(5) = wvno2 - ra*rb
mmm1 = mmax - 1
do m = mmm1, llw, -1
    xka = omega/dble(a(m))
    xkb = omega/dble(b(m))
    t = dble(b(m))/omega
    gammk = 2.d+00*t*t
    gam = gammk*wvno2
    wvnop = wvno + xka
    wvnom = dabs(wvno - xka)
    ra = dsqrt(wvnop*wvnom)
    wvnop = wvno + xkb
    wvnom = dabs(wvno - xkb)
    rb = dsqrt(wvnop*wvnom)
    dpth = dble(d(m))
    rho1 = dble(rho(m))
    p = ra*dpth
    q = rb*dpth
    beta = dble(b(m))
    call var(p, q, ra, rb, wvno, xka, xkb, dpth, w, cosp, exa, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
    call dnka(ca, wvno2, gam, gammk, rho1, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
    do i = 1, 5
        cr = 0.0d+00
        do j = 1, 5
            cr = cr + e(j)*ca(j, i)
            100 continue
        end do
        ee(i) = cr
        200 continue
    end do
    call normc(ee, exa)
    do i = 1, 5
        e(i) = ee(i)
        300 continue
    end do
    500 continue
end do
if (llw /= 1) then
    xka = omega/dble(a(1))
    wvnop = wvno + xka
    wvnom = dabs(wvno - xka)
    ra = dsqrt(wvnop*wvnom)
    dpth = dble(d(1))
    rho1 = dble(rho(1))
    p = ra*dpth
    beta = dble(b(1))
    znul = 1.0d-05
    call var(p, znul, ra, znul, wvno, xka, znul, dpth, w, cosp, exa, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
    w0 = (-rho1)*w
    dltar4 = cosp*e(1) + w0*e(2)
else
    dltar4 = e(1)
end if
return
end function dltar4

subroutine var(p, q, ra, rb, wvno, xka, xkb, dpth, w, cosp, exa, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
implicit double precision (a-h,o-z)
exa = 0.0d+00
a0 = 1.0d+00
pex = 0.0d+00
sex = 0.0d+00
if (wvno < xka) then
    sinp = dsin(p)
    w = sinp/ra
    x = (-ra)*sinp
    cosp = dcos(p)
else
    if (wvno == xka) then
        cosp = 1.0d+00
        w = dpth
        x = 0.0d+00
    else
        if (wvno > xka) then
            pex = p
            fac = 0.0d+00
            if (p < 16) then
                fac = dexp((-2.0d+00)*p)
            end if
            cosp = (1.0d+00 + fac)*0.5d+00
            sinp = (1.0d+00 - fac)*0.5d+00
            w = sinp/ra
            x = ra*sinp
        end if
    end if
end if
if (wvno < xkb) then
    sinq = dsin(q)
    y = sinq/rb
    z = (-rb)*sinq
    cosq = dcos(q)
else
    if (wvno == xkb) then
        cosq = 1.0d+00
        y = dpth
        z = 0.0d+00
    else
        if (wvno > xkb) then
            sex = q
            fac = 0.0d+00
            if (q < 16) then
                fac = dexp((-2.0d+0)*q)
            end if
            cosq = (1.0d+00 + fac)*0.5d+00
            sinq = (1.0d+00 - fac)*0.5d+00
            y = sinq/rb
            z = rb*sinq
        end if
    end if
end if
exa = pex + sex
a0 = 0.0d+00
if (exa < 60.0d+00) then
    a0 = dexp(-exa)
end if
cpcq = cosp*cosq
cpy = cosp*y
cpz = cosp*z
cqw = cosq*w
cqx = cosq*x
xy = x*y
xz = x*z
wy = w*y
wz = w*z
qmp = sex - pex
fac = 0.0d+00
if (qmp > -40.0d+00) then
    fac = dexp(qmp)
end if
cosq = cosq*fac
y = fac*y
z = fac*z
return
end subroutine var

subroutine normc(ee, ex)
implicit double precision (a-h,o-z)
dimension :: ee(5)
ex = 0.0d+00
t1 = 0.0d+00
do i = 1, 5
    if (dabs(ee(i)) > t1) then
        t1 = dabs(ee(i))
    end if
    10 continue
end do
if (t1 < 1.d-40) then
    t1 = 1.d+00
end if
do i = 1, 5
    t2 = ee(i)
    t2 = t2/t1
    ee(i) = t2
    20 continue
end do
ex = dlog(t1)
return
end subroutine normc

subroutine dnka(ca, wvno2, gam, gammk, rho, a0, cpcq, cpy, cpz, cqw, cqx, xy, xz, wy, wz)
implicit double precision (a-h,o-z)
dimension :: ca(5,5)
data one, two/1.d+00, 2.d+00/
gamm1 = gam - one
twgm1 = gam + gamm1
gmgmk = gam*gammk
gmgm1 = gam*gamm1
gm1sq = gamm1*gamm1
rho2 = rho*rho
a0pq = a0 - cpcq
ca(1, 1) = cpcq - two*gmgm1*a0pq - gmgmk*xz - wvno2*gm1sq*wy
ca(1, 2) = (wvno2*cpy - cqx)/rho
ca(1, 3) = (-(twgm1*a0pq + gammk*xz + wvno2*gamm1*wy))/rho
ca(1, 4) = (cpz - wvno2*cqw)/rho
ca(1, 5) = (-(two*wvno2*a0pq + xz + wvno2*wvno2*wy))/rho2
ca(2, 1) = (gmgmk*cpz - gm1sq*cqw)*rho
ca(2, 2) = cpcq
ca(2, 3) = gammk*cpz - gamm1*cqw
ca(2, 4) = -wz
ca(2, 5) = ca(1, 4)
ca(4, 1) = (gm1sq*cpy - gmgmk*cqx)*rho
ca(4, 2) = -xy
ca(4, 3) = gamm1*cpy - gammk*cqx
ca(4, 4) = ca(2, 2)
ca(4, 5) = ca(1, 2)
ca(5, 1) = (-(two*gmgmk*gm1sq*a0pq + gmgmk*gmgmk*xz + gm1sq*gm1sq*wy))*rho2
ca(5, 2) = ca(4, 1)
ca(5, 3) = (-(gammk*gamm1*twgm1*a0pq + gam*gammk*gammk*xz + gamm1*gm1sq*wy))*rho
ca(5, 4) = ca(2, 1)
ca(5, 5) = ca(1, 1)
t = (-two)*wvno2
ca(3, 1) = t*ca(5, 3)
ca(3, 2) = t*ca(4, 3)
ca(3, 3) = a0 + two*(cpcq - ca(1, 1))
ca(3, 4) = t*ca(2, 3)
ca(3, 5) = t*ca(1, 3)
return
end subroutine dnka