37subroutine ylmscale(p, vscales, v4pi2lp1, vscales_rel)
 
   39    integer, 
intent(in) :: p
 
   41    real(dp), 
intent(out) :: vscales((p+1)**2), v4pi2lp1(p+1), &
 
   42        & vscales_rel((p+1)**2)
 
   44    real(dp) :: tmp, twolp1
 
   53        vscales_rel(ind) = tmp
 
   54        vscales(ind) = one / tmp
 
   56        tmp = vscales(ind) * sqrt2
 
   59            tmp = -tmp / sqrt(dble((l-m+1)*(l+m)))
 
   62            vscales_rel(ind+m) = tmp * v4pi2lp1(l+1)
 
   63            vscales_rel(ind-m) = vscales_rel(ind+m)
 
   81        & m2l_ztranslate_adj_coef)
 
   83    integer, 
intent(in) :: dmax, pm, pl
 
   85    real(dp), 
intent(out) :: vcnk((2*dmax+1)*(dmax+1)), &
 
   86        & m2l_ztranslate_coef(pm+1, pl+1, pl+1), &
 
   87        & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1)
 
   89    integer :: i, indi, j, k, n, indjn
 
  104            vcnk(indi+j) = vcnk(indi-i+j+1) + vcnk(indi-i+j)
 
  112            vcnk(indi+j) = sqrt(vcnk(indi+j))
 
  120                indjn = (j+n)*(j+n+1)/2 + 1
 
  121                m2l_ztranslate_coef(n-k+1, k+1, j-k+1) = &
 
  122                    & tmp1 * vcnk(indjn+j-k) * vcnk(indjn+j+k)
 
  123                m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1) = &
 
  124                    & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)
 
  144    real(dp), 
intent(in) :: x(3)
 
  146    real(dp), 
intent(out) :: rho, ctheta, stheta, cphi, sphi
 
  148    real(dp) :: max12, ssq12
 
  150    if ((x(1) .eq. zero) .and. (x(2) .eq. zero)) 
then 
  152        ctheta = sign(one, x(3))
 
  160    if (abs(x(2)) .gt. abs(x(1))) 
then 
  162        ssq12 = one + (x(1)/x(2))**2
 
  165        ssq12 = one + (x(2)/x(1))**2
 
  167    stheta = max12 * sqrt(ssq12)
 
  171    if (abs(x(3)) .gt. max12) 
then 
  172        rho = one + ssq12*(max12/x(3))**2
 
  173        rho = abs(x(3)) * sqrt(rho)
 
  174        stheta = stheta / rho
 
  177        rho = ssq12 + (x(3)/max12)**2
 
  178        rho = max12 * sqrt(rho)
 
  179        stheta = stheta / rho
 
  214subroutine ylmbas(x, rho, ctheta, stheta, cphi, sphi, p, vscales, vylm, vplm, &
 
  217    integer, 
intent(in) :: p
 
  218    real(dp), 
intent(in) :: x(3)
 
  219    real(dp), 
intent(in) :: vscales((p+1)**2)
 
  221    real(dp), 
intent(out) :: rho, ctheta, stheta, cphi, sphi
 
  222    real(dp), 
intent(out) :: vylm((p+1)**2), vplm((p+1)**2)
 
  223    real(dp), 
intent(out) :: vcos(p+1), vsin(p+1)
 
  226    real(dp) :: max12, ssq12, tmp
 
  231    if (x(1) .eq. zero) 
then 
  234    else if (abs(x(2)) .gt. abs(x(1))) 
then 
  236        ssq12 = one + (x(1)/x(2))**2
 
  239        ssq12 = one + (x(2)/x(1))**2
 
  242    if (x(3) .eq. zero) 
then 
  243        rho = max12 * sqrt(ssq12)
 
  244    else if (abs(x(3)) .gt. max12) 
then 
  245        rho = one + ssq12 *(max12/x(3))**2
 
  246        rho = abs(x(3)) * sqrt(rho)
 
  248        rho = ssq12 + (x(3)/max12)**2
 
  249        rho = max12 * sqrt(rho)
 
  252    if (rho .eq. zero) 
then 
  256    stheta = max12 * sqrt(ssq12)
 
  258    if (stheta .ne. zero) 
then 
  264        stheta = stheta / rho
 
  272                vylm(2) = - vscales(4) * stheta * sphi
 
  273                vylm(3) = vscales(3) * ctheta
 
  274                vylm(4) = - vscales(4) * stheta * cphi
 
  281        call trgev(cphi, sphi, p, vcos, vsin)
 
  286        vylm(2) = - vscales(4) * stheta * sphi
 
  287        vylm(3) = vscales(3) * ctheta
 
  288        vylm(4) = - vscales(4) * stheta * cphi
 
  295            vylm(ind) = vscales(ind) * vplm(ind)
 
  298                tmp = vplm(ind+m) * vscales(ind+m)
 
  300                vylm(ind+m) = tmp * vcos(m+1)
 
  302                vylm(ind-m) = tmp * vsin(m+1)
 
  310        ctheta = sign(one, x(3))
 
  324            vylm(ind) = vscales(ind)
 
  330            vylm(ind) = ctheta * vscales(ind)
 
  348subroutine trgev(cphi, sphi, p, vcos, vsin)
 
  350    integer, 
intent(in) :: p
 
  351    real(dp), 
intent(in) :: cphi, sphi
 
  353    real(dp), 
intent(out) :: vcos(p+1), vsin(p+1)
 
  356    real(dp) :: c4phi, s4phi
 
  380            vcos(3) = cphi**2 - sphi**2
 
  381            vsin(3) = 2 * cphi * sphi
 
  389            vcos(3) = cphi**2 - sphi**2
 
  390            vsin(3) = 2 * cphi * sphi
 
  391            vcos(4) = vcos(3)*cphi - vsin(3)*sphi
 
  392            vsin(4) = vcos(3)*sphi + vsin(3)*cphi
 
  400            vcos(3) = cphi**2 - sphi**2
 
  401            vsin(3) = 2 * cphi * sphi
 
  402            vcos(4) = vcos(3)*cphi - vsin(3)*sphi
 
  403            vsin(4) = vcos(3)*sphi + vsin(3)*cphi
 
  404            vcos(5) = vcos(3)**2 - vsin(3)**2
 
  405            vsin(5) = 2 * vcos(3) * vsin(3)
 
  411        vcos(m:m+3) = vcos(m-4:m-1)*c4phi - vsin(m-4:m-1)*s4phi
 
  412        vsin(m:m+3) = vcos(m-4:m-1)*s4phi + vsin(m-4:m-1)*c4phi
 
  417            vcos(p-1) = vcos(p-2)*cphi - vsin(p-2)*sphi
 
  418            vsin(p-1) = vcos(p-2)*sphi + vsin(p-2)*cphi
 
  419            vcos(p) = vcos(p-2)*vcos(3) - vsin(p-2)*vsin(3)
 
  420            vsin(p) = vcos(p-2)*vsin(3) + vsin(p-2)*vcos(3)
 
  421            vcos(p+1) = vcos(p-2)*vcos(4) - vsin(p-2)*vsin(4)
 
  422            vsin(p+1) = vcos(p-2)*vsin(4) + vsin(p-2)*vcos(4)
 
  424            vcos(p) = vcos(p-1)*cphi - vsin(p-1)*sphi
 
  425            vsin(p) = vcos(p-1)*sphi + vsin(p-1)*cphi
 
  426            vcos(p+1) = vcos(p-1)*vcos(3) - vsin(p-1)*vsin(3)
 
  427            vsin(p+1) = vcos(p-1)*vsin(3) + vsin(p-1)*vcos(3)
 
  429            vcos(p+1) = vcos(p)*cphi - vsin(p)*sphi
 
  430            vsin(p+1) = vcos(p)*sphi + vsin(p)*cphi
 
  452subroutine polleg(ctheta, stheta, p, vplm)
 
  454    integer, 
intent(in) :: p
 
  455    real(dp), 
intent(in) :: ctheta, stheta
 
  457    real(dp), 
intent(out) :: vplm((p+1)**2)
 
  459    real(dp) :: work(p+1)
 
  484    integer, 
intent(in) :: p
 
  485    real(dp), 
intent(in) :: ctheta, stheta
 
  487    real(dp), 
intent(out) :: vplm((p+1)**2)
 
  489    real(dp), 
intent(out) :: work(p+1)
 
  492    real(dp) :: tmp1, tmp2, tmp3, tmp4, pl2m, pl1m, plm, pmm
 
  507            tmp1 = three * stheta
 
  508            tmp2 = ctheta * ctheta
 
  509            vplm(7) = 1.5d0*tmp2 - pt5
 
  510            vplm(8) = -tmp1 * ctheta
 
  511            vplm(9) = tmp1 * stheta
 
  517            tmp1 = three * stheta
 
  518            tmp2 = ctheta * ctheta
 
  519            vplm(7) = 1.5d0*tmp2 - pt5
 
  520            vplm(8) = -tmp1 * ctheta
 
  521            vplm(9) = tmp1 * stheta
 
  522            tmp3 = 2.5d0 * tmp2 - pt5
 
  523            vplm(13) = tmp3*ctheta - ctheta
 
  524            vplm(14) = -tmp1 * tmp3
 
  526            vplm(15) = tmp4 * vplm(8)
 
  527            vplm(16) = tmp4 * vplm(9)
 
  533        work(l) = dble(2*l-1) * ctheta
 
  544        tmp1 = dble(2*m+1) * pmm
 
  550        vplm((m+1)*(m+3)) = pl1m
 
  553            plm = work(l)*pl1m - dble(l+m-1)*pl2m
 
  554            plm = plm / dble(l-m)
 
  555            vplm(l*l+l+1+m) = plm
 
  572subroutine modified_spherical_bessel_first_kind(lmax, argument, SI, DI, work)
 
  575    integer, 
intent(in) :: lmax
 
  576    real(dp), 
intent(in) :: argument
 
  578    real(dp), 
dimension(0:lmax), 
intent(out) :: SI, DI
 
  580    complex(dp), 
dimension(0:max(1, lmax)), 
intent(out) :: work
 
  589    complex(dp) :: complex_argument
 
  590    real(dp) :: scaling_factor
 
  591    real(dp), 
parameter :: fnu = 1.5d0
 
  592    integer :: NZ, ierr, l
 
  594    if (argument .eq. zero) 
then 
  603    scaling_factor = sqrt(pi/(2*argument))
 
  605    complex_argument = argument
 
  607    call cbesi(complex_argument, pt5, 1, 1, work(0), nz, ierr)
 
  608    if (ierr .ne. 0) stop 
'Error in computing Bessel function of first kind' 
  610    if (lmax .gt. 0) 
then 
  611        call cbesi(complex_argument, fnu, 1, lmax, work(1:lmax), nz, ierr)
 
  612        if (ierr .ne. 0) stop 
'Error in computing Bessel function of first kind' 
  615        call cbesi(complex_argument, fnu, 1, 1, work(1), nz, ierr)
 
  616        if (ierr .ne. 0) stop 
'Error in computing Bessel function of first kind' 
  619    si = real(work(0:lmax))
 
  621    si = scaling_factor * si
 
  623    di(0) = scaling_factor * real(work(1))
 
  625        di(l)= si(l-1) - ((l+1.0d0)*si(l))/argument
 
  627end subroutine modified_spherical_bessel_first_kind
 
  635subroutine modified_spherical_bessel_second_kind(lmax, argument, SK, DK, work)
 
  638    integer , 
intent(in) :: lmax
 
  639    real(dp), 
intent(in) :: argument
 
  641    real(dp), 
dimension(0:lmax), 
intent(out) :: SK, DK
 
  643    complex(dp), 
dimension(0:max(1, lmax)), 
intent(out) :: work
 
  652    complex(dp) :: complex_argument
 
  653    real(dp) :: scaling_factor
 
  654    real(dp), 
parameter :: fnu=1.5d0
 
  655    integer :: NZ, ierr, l
 
  658    scaling_factor = sqrt(2/(pi*argument))
 
  660    complex_argument = argument
 
  662    call cbesk(complex_argument, pt5, 1, 1, work(0), nz, ierr)
 
  663    if (ierr .ne. 0) stop 
'Error in computing Bessel function of second kind' 
  665    if (lmax .gt. 0) 
then 
  666        call cbesk(complex_argument, fnu, 1, lmax, work(1:lmax), nz, ierr)
 
  667        if (ierr .ne. 0) stop 
'Error in computing Bessel function of second kind' 
  670        call cbesk(complex_argument, fnu, 1, 1, work(1), nz, ierr)
 
  671        if (ierr .ne. 0) stop 
'Error in computing Bessel function of second kind' 
  674    sk = real(work(0:lmax))
 
  676    sk = scaling_factor * sk
 
  678    dk(0) = -scaling_factor * real(work(1))
 
  680        dk(l) = -sk(l-1) - ((l+1.0d0)*sk(l))/argument
 
  682end subroutine modified_spherical_bessel_second_kind
 
  710subroutine fmm_p2m(c, src_q, dst_r, p, vscales, beta, dst_m)
 
  712    integer, 
intent(in) :: p
 
  713    real(dp), 
intent(in) :: c(3), src_q, dst_r, vscales((p+1)**2), beta
 
  715    real(dp), 
intent(inout) :: dst_m((p+1)**2)
 
  717    real(dp) :: work(2*(p+1)*(p+2))
 
  719    call fmm_p2m_work(c, src_q, dst_r, p, vscales, beta, dst_m, work)
 
  751    integer, 
intent(in) :: p
 
  752    real(dp), 
intent(in) :: c(3), src_q, dst_r, vscales((p+1)**2), beta
 
  754    real(dp), 
intent(inout) :: dst_m((p+1)**2)
 
  756    real(dp), 
intent(out), 
target :: work(2*(p+1)*(p+2))
 
  758    real(dp) :: rho, ctheta, stheta, cphi, sphi, t, rcoef
 
  759    integer :: n, ind, vylm1, vylm2, vplm1, vplm2, vcos1, vcos2, vsin1, vsin2
 
  761    if (src_q .eq. zero) 
then 
  763        if (beta .eq. zero) 
then 
  783    call ylmbas(c, rho, ctheta, stheta, cphi, sphi, p, vscales, &
 
  784        & work(vylm1:vylm2), work(vplm1:vplm2), work(vcos1:vcos2), &
 
  787    if (rho .ne. zero) 
then 
  791        if (beta .eq. zero) 
then 
  795                dst_m(ind-n:ind+n) = t * work(ind-n:ind+n)
 
  803                dst_m(ind-n:ind+n) = beta*dst_m(ind-n:ind+n) + &
 
  804                    & t*work(ind-n:ind+n)
 
  811        if (beta .eq. zero) 
then 
  812            dst_m(1) = src_q / dst_r / sqrt4pi
 
  816            dst_m(1) = beta*dst_m(1) + src_q/dst_r/sqrt4pi
 
  817            dst_m(2:) = beta * dst_m(2:)
 
  848subroutine fmm_m2p(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v)
 
  850    integer, 
intent(in) :: p
 
  851    real(dp), 
intent(in) :: c(3), src_r, vscales_rel((p+1)*(p+1)), alpha, &
 
  852        & src_m((p+1)*(p+1)), beta
 
  854    real(dp), 
intent(inout) :: dst_v
 
  856    real(dp) :: work(p+1)
 
  858    call fmm_m2p_work(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v, &
 
  894    integer, 
intent(in) :: p
 
  895    real(dp), 
intent(in) :: c(3), src_r, vscales_rel((p+1)*(p+1)), alpha, &
 
  896        & src_m((p+1)*(p+1)), beta
 
  898    real(dp), 
intent(inout) :: dst_v
 
  900    real(dp), 
intent(out), 
target :: work(p+1)
 
  902    real(dp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
 
  903        & tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
 
  904    integer :: l, m, indl
 
  906    if (beta .eq. zero) 
then 
  913    if (alpha .eq. zero) 
then 
  917    if (c(1) .eq. zero) 
then 
  920    else if (abs(c(2)) .gt. abs(c(1))) 
then 
  922        ssq12 = one + (c(1)/c(2))**2
 
  925        ssq12 = one + (c(2)/c(1))**2
 
  928    if (c(3) .eq. zero) 
then 
  929        rho = max12 * sqrt(ssq12)
 
  930    else if (abs(c(3)) .gt. max12) 
then 
  931        rho = one + ssq12 *(max12/c(3))**2
 
  932        rho = abs(c(3)) * sqrt(rho)
 
  934        rho = ssq12 + (c(3)/max12)**2
 
  935        rho = max12 * sqrt(rho)
 
  939    if (rho .eq. zero) 
then 
  945    stheta = max12 * sqrt(ssq12)
 
  947    if (stheta .ne. zero) 
then 
  953        stheta = stheta / rho
 
  958                dst_v = dst_v + alpha*rcoef*vscales_rel(1)*src_m(1)
 
  962                tmp = src_m(1) * vscales_rel(1)
 
  964                tmp2 = ctheta * src_m(3) * vscales_rel(3)
 
  965                tmp3 = vscales_rel(4) * stheta
 
  966                tmp2 = tmp2 - tmp3*sphi*src_m(2)
 
  967                tmp2 = tmp2 - tmp3*cphi*src_m(4)
 
  968                dst_v = dst_v + alpha*rcoef*(tmp+rcoef*tmp2)
 
  973        work(1) = alpha * rcoef
 
  975            work(l+1) = rcoef * work(l)
 
  980        dst_v = dst_v + work(1)*src_m(1)*vscales_rel(1)
 
  985        ylm = pl1m * vscales_rel(3)
 
  986        dst_v = dst_v + work(2)*src_m(3)*ylm
 
  989            plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
  991            ylm = plm * vscales_rel(l*l+l+1)
 
  992            dst_v = dst_v + work(l+1)*src_m(l*l+l+1)*ylm
 
 1004            ylm = pmm * vscales_rel((m+1)**2)
 
 1005            tmp1 = cmphi*src_m((m+1)**2) + smphi*src_m(m*m+1)
 
 1006            dst_v = dst_v + work(m+1)*ylm*tmp1
 
 1008            tmp1 = dble(2*m+1) * pmm
 
 1010            pmm = -stheta * tmp1
 
 1012            pl1m = ctheta * tmp1
 
 1013            ylm = pl1m * vscales_rel((m+1)*(m+3))
 
 1014            tmp1 = cmphi*src_m((m+1)*(m+3)) + smphi*src_m((m+1)*(m+2)+1-m)
 
 1015            dst_v = dst_v + work(m+2)*ylm*tmp1
 
 1018                plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 1019                plm = plm / dble(l-m)
 
 1020                ylm = plm * vscales_rel(l*l+l+1+m)
 
 1021                tmp1 = cmphi*src_m(l*l+l+1+m) + smphi*src_m(l*l+l+1-m)
 
 1022                dst_v = dst_v + work(l+1)*ylm*tmp1
 
 1029            cmphi = cmphi*cphi - smphi*sphi
 
 1030            smphi = tmp1*sphi + smphi*cphi
 
 1033        ylm = pmm * vscales_rel((p+1)**2)
 
 1034        tmp1 = cmphi*src_m((p+1)**2) + smphi*src_m(p*p+1)
 
 1035        dst_v = dst_v + work(p+1)*ylm*tmp1
 
 1044        if (c(3) .lt. zero) 
then 
 1053            dst_v = dst_v + t*src_m(indl)*vscales_rel(indl)
 
 1089    integer, 
intent(in) :: p
 
 1090    real(dp), 
intent(in) :: c(3), src_r, vscales((p+1)*(p+1)), alpha, &
 
 1091        & src_m((p+1)*(p+1)), beta
 
 1093    real(dp), 
intent(inout) :: dst_v
 
 1095    complex(dp) :: work_complex(max(2, p+1))
 
 1096    real(dp) :: work(p+1)
 
 1098    real(dp) :: src_sk(p+1)
 
 1100    call modified_spherical_bessel_second_kind(p, src_r, src_sk, work, &
 
 1103        & dst_v, work_complex, work)
 
 1133        & dst_v, work_complex, work)
 
 1136    integer, 
intent(in) :: p
 
 1137    real(dp), 
intent(in) :: c(3), vscales((p+1)*(p+1)), alpha, &
 
 1138        & src_m((p+1)*(p+1)), beta, SK_ri(p+1)
 
 1140    real(dp), 
intent(inout) :: dst_v
 
 1142    complex(dp), 
intent(out) :: work_complex(p+1)
 
 1143    real(dp), 
intent(out) :: work(p+1)
 
 1145    real(dp) :: rho, ctheta, stheta, cphi, sphi, tmp, tmp1, tmp2, &
 
 1146        & tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm, t
 
 1147    real(dp) :: scaling_factor
 
 1148    complex(dp) :: complex_argument
 
 1149    integer :: l, m, indl, ierr, NZ
 
 1150    real(dp), 
parameter :: fnu=1.5d0
 
 1152    if (beta .eq. zero) 
then 
 1155        dst_v = beta * dst_v
 
 1159    if (alpha .eq. zero) 
then 
 1163    if (c(1) .eq. zero) 
then 
 1166    else if (abs(c(2)) .gt. abs(c(1))) 
then 
 1168        ssq12 = one + (c(1)/c(2))**2
 
 1171        ssq12 = one + (c(2)/c(1))**2
 
 1174    if (c(3) .eq. zero) 
then 
 1175        rho = max12 * sqrt(ssq12)
 
 1176    else if (abs(c(3)) .gt. max12) 
then 
 1177        rho = one + ssq12 *(max12/c(3))**2
 
 1178        rho = abs(c(3)) * sqrt(rho)
 
 1180        rho = ssq12 + (c(3)/max12)**2
 
 1181        rho = max12 * sqrt(rho)
 
 1185    if (rho .eq. zero) 
then 
 1189    scaling_factor = alpha * sqrt(two / (pi*rho))
 
 1191    complex_argument = rho
 
 1193    call cbesk(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 1196        call cbesk(complex_argument, fnu, 1, p, work_complex(2:p+1), nz, ierr)
 
 1198    work = scaling_factor * real(work_complex) / sk_ri
 
 1201    stheta = max12 * sqrt(ssq12)
 
 1203    if (stheta .ne. zero) 
then 
 1205        cphi = c(1) / stheta
 
 1206        sphi = c(2) / stheta
 
 1209        stheta = stheta / rho
 
 1214                dst_v = dst_v + work(1)*vscales(1)*src_m(1)
 
 1218                tmp = work(1) * src_m(1) * vscales(1)
 
 1220                tmp2 = ctheta * src_m(3) * vscales(3)
 
 1221                tmp3 = stheta * vscales(4)
 
 1222                tmp2 = tmp2 - tmp3*sphi*src_m(2)
 
 1223                tmp2 = tmp2 - tmp3*cphi*src_m(4)
 
 1224                dst_v = dst_v + tmp + work(2)*tmp2
 
 1231        dst_v = dst_v + work(1)*src_m(1)*vscales(1)
 
 1236        ylm = pl1m * vscales(3)
 
 1237        dst_v = dst_v + work(2)*src_m(3)*ylm
 
 1240            plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 1242            ylm = plm * vscales(l*l+l+1)
 
 1243            dst_v = dst_v + work(l+1)*src_m(l*l+l+1)*ylm
 
 1255            ylm = pmm * vscales((m+1)**2)
 
 1256            tmp1 = cmphi*src_m((m+1)**2) + smphi*src_m(m*m+1)
 
 1257            dst_v = dst_v + work(m+1)*ylm*tmp1
 
 1259            tmp1 = dble(2*m+1) * pmm
 
 1261            pmm = -stheta * tmp1
 
 1263            pl1m = ctheta * tmp1
 
 1264            ylm = pl1m * vscales((m+1)*(m+3))
 
 1265            tmp1 = cmphi*src_m((m+1)*(m+3)) + smphi*src_m((m+1)*(m+2)+1-m)
 
 1266            dst_v = dst_v + work(m+2)*ylm*tmp1
 
 1269                plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 1270                plm = plm / dble(l-m)
 
 1271                ylm = plm * vscales(l*l+l+1+m)
 
 1272                tmp1 = cmphi*src_m(l*l+l+1+m) + smphi*src_m(l*l+l+1-m)
 
 1273                dst_v = dst_v + work(l+1)*ylm*tmp1
 
 1280            cmphi = cmphi*cphi - smphi*sphi
 
 1281            smphi = tmp1*sphi + smphi*cphi
 
 1284        ylm = pmm * vscales((p+1)**2)
 
 1285        tmp1 = cmphi*src_m((p+1)**2) + smphi*src_m(p*p+1)
 
 1286        dst_v = dst_v + work(p+1)*ylm*tmp1
 
 1295        if (c(3) .lt. zero) 
then 
 1302                dst_v = dst_v + t*work(l+1)*src_m(indl)*vscales(indl)
 
 1312                dst_v = dst_v + work(l+1)*src_m(indl)*vscales(indl)
 
 1318subroutine fmm_m2p_bessel_grad(c, src_r, p, vscales, alpha, src_m, beta, dst_g)
 
 1321    integer, 
intent(in) :: p
 
 1322    real(dp), 
intent(in) :: c(3), src_r, vscales((p+2)*(p+2)), alpha, &
 
 1323        & src_m((p+1)*(p+1)), beta
 
 1325    real(dp), 
intent(inout) :: dst_g(3)
 
 1327    complex(dp) :: work_complex(p+2)
 
 1328    real(dp) :: work(p+2), src_sk(p+2), src_m_grad((p+2)**2, 3)
 
 1330    call modified_spherical_bessel_second_kind(p+1, src_r, src_sk, work, &
 
 1332    call fmm_m2m_bessel_grad(p, src_sk, vscales, src_m, src_m_grad)
 
 1334        & beta, dst_g(1), work_complex, work)
 
 1336        & beta, dst_g(2), work_complex, work)
 
 1338        & beta, dst_g(3), work_complex, work)
 
 1339end subroutine fmm_m2p_bessel_grad
 
 1368    integer, 
intent(in) :: p
 
 1369    real(dp), 
intent(in) :: c(3), src_q, dst_r, vscales_rel((p+1)*(p+1)), beta
 
 1371    real(dp), 
intent(inout) :: dst_m((p+1)**2)
 
 1373    real(dp) :: work(p+1)
 
 1406    integer, 
intent(in) :: p
 
 1407    real(dp), 
intent(in) :: c(3), src_q, dst_r, vscales_rel((p+1)*(p+1)), beta
 
 1409    real(dp), 
intent(inout) :: dst_m((p+1)**2)
 
 1411    real(dp), 
intent(out), 
target :: work(p+1)
 
 1413    real(dp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
 
 1414        & max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
 
 1415    integer :: l, m, indl
 
 1417    if (src_q .eq. zero) 
then 
 1419        if (beta .eq. zero) 
then 
 1423            dst_m = beta * dst_m
 
 1430    if (c(1) .eq. zero) 
then 
 1433    else if (abs(c(2)) .gt. abs(c(1))) 
then 
 1435        ssq12 = one + (c(1)/c(2))**2
 
 1438        ssq12 = one + (c(2)/c(1))**2
 
 1441    if (c(3) .eq. zero) 
then 
 1442        rho = max12 * sqrt(ssq12)
 
 1443    else if (abs(c(3)) .gt. max12) 
then 
 1444        rho = one + ssq12 *(max12/c(3))**2
 
 1445        rho = abs(c(3)) * sqrt(rho)
 
 1447        rho = ssq12 + (c(3)/max12)**2
 
 1448        rho = max12 * sqrt(rho)
 
 1452    if (rho .eq. zero) 
then 
 1458    stheta = max12 * sqrt(ssq12)
 
 1460    if (stheta .ne. zero) 
then 
 1462        cphi = c(1) / stheta
 
 1463        sphi = c(2) / stheta
 
 1466        stheta = stheta / rho
 
 1471                if (beta .eq. zero) 
then 
 1472                    dst_m(1) = t * vscales_rel(1)
 
 1474                    dst_m(1) = beta*dst_m(1) + t*vscales_rel(1)
 
 1478                if (beta .eq. zero) 
then 
 1480                    dst_m(1) = t * vscales_rel(1)
 
 1483                    tmp = -vscales_rel(4) * stheta
 
 1485                    dst_m(2) = tmp2 * sphi
 
 1486                    dst_m(3) = t * vscales_rel(3) * ctheta
 
 1487                    dst_m(4) = tmp2 * cphi
 
 1490                    dst_m(1) = beta*dst_m(1) + t*vscales_rel(1)
 
 1493                    tmp = -vscales_rel(4) * stheta
 
 1495                    dst_m(2) = beta*dst_m(2) + tmp2*sphi
 
 1496                    dst_m(3) = beta*dst_m(3) + t*vscales_rel(3)*ctheta
 
 1497                    dst_m(4) = beta*dst_m(4) + tmp2*cphi
 
 1503        work(1) = src_q * rcoef
 
 1505            work(l+1) = rcoef * work(l)
 
 1508        if (beta .eq. zero) 
then 
 1512            dst_m(1) = work(1) * vscales_rel(1)
 
 1517            ylm = pl1m * vscales_rel(3)
 
 1518            dst_m(3) = work(2) * ylm
 
 1521                plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 1523                ylm = plm * vscales_rel(l*l+l+1)
 
 1524                dst_m(l*l+l+1) = work(l+1) * ylm
 
 1536                ylm = pmm * vscales_rel((m+1)**2)
 
 1537                ylm = work(m+1) * ylm
 
 1538                dst_m((m+1)**2) = cmphi * ylm
 
 1539                dst_m(m*m+1) = smphi * ylm
 
 1541                tmp1 = dble(2*m+1) * pmm
 
 1543                pmm = -stheta * tmp1
 
 1545                pl1m = ctheta * tmp1
 
 1546                ylm = pl1m * vscales_rel((m+1)*(m+3))
 
 1547                ylm = work(m+2) * ylm
 
 1548                dst_m((m+1)*(m+3)) = cmphi * ylm
 
 1549                dst_m((m+1)*(m+2)+1-m) = smphi * ylm
 
 1552                    plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 1553                    plm = plm / dble(l-m)
 
 1554                    ylm = plm * vscales_rel(l*l+l+1+m)
 
 1555                    ylm = work(l+1) * ylm
 
 1556                    dst_m(l*l+l+1+m) = cmphi * ylm
 
 1557                    dst_m(l*l+l+1-m) = smphi * ylm
 
 1564                cmphi = cmphi*cphi - smphi*sphi
 
 1565                smphi = tmp1*sphi + smphi*cphi
 
 1568            ylm = pmm * vscales_rel((p+1)**2)
 
 1569            ylm = work(p+1) * ylm
 
 1570            dst_m((p+1)**2) = cmphi * ylm
 
 1571            dst_m(p*p+1) = smphi * ylm
 
 1577            dst_m(1) = beta*dst_m(1) + work(1)*vscales_rel(1)
 
 1582            ylm = pl1m * vscales_rel(3)
 
 1583            dst_m(3) = beta*dst_m(3) + work(2)*ylm
 
 1586                plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 1588                ylm = plm * vscales_rel(l*l+l+1)
 
 1589                dst_m(l*l+l+1) = beta*dst_m(l*l+l+1) + work(l+1)*ylm
 
 1601                ylm = pmm * vscales_rel((m+1)**2)
 
 1602                ylm = work(m+1) * ylm
 
 1603                dst_m((m+1)**2) = beta*dst_m((m+1)**2) + cmphi*ylm
 
 1604                dst_m(m*m+1) = beta*dst_m(m*m+1) + smphi*ylm
 
 1606                tmp1 = dble(2*m+1) * pmm
 
 1608                pmm = -stheta * tmp1
 
 1610                pl1m = ctheta * tmp1
 
 1611                ylm = pl1m * vscales_rel((m+1)*(m+3))
 
 1612                ylm = work(m+2) * ylm
 
 1613                dst_m((m+1)*(m+3)) = beta*dst_m((m+1)*(m+3)) + cmphi*ylm
 
 1614                dst_m((m+1)*(m+2)+1-m) = beta*dst_m((m+1)*(m+2)+1-m) + &
 
 1618                    plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 1619                    plm = plm / dble(l-m)
 
 1620                    ylm = plm * vscales_rel(l*l+l+1+m)
 
 1621                    ylm = work(l+1) * ylm
 
 1622                    dst_m(l*l+l+1+m) = beta*dst_m(l*l+l+1+m) + cmphi*ylm
 
 1623                    dst_m(l*l+l+1-m) = beta*dst_m(l*l+l+1-m) + smphi*ylm
 
 1630                cmphi = cmphi*cphi - smphi*sphi
 
 1631                smphi = tmp1*sphi + smphi*cphi
 
 1634            ylm = pmm * vscales_rel((p+1)**2)
 
 1635            ylm = work(p+1) * ylm
 
 1636            dst_m((p+1)**2) = beta*dst_m((p+1)**2) + cmphi*ylm
 
 1637            dst_m(p*p+1) = beta*dst_m(p*p+1) + smphi*ylm
 
 1647        if (c(3) .lt. zero) 
then 
 1652        if (beta .eq. zero) 
then 
 1657                dst_m(indl) = t * vscales_rel(indl)
 
 1658                dst_m(indl-l:indl-1) = zero
 
 1659                dst_m(indl+1:indl+l) = zero
 
 1668                dst_m(indl) = beta*dst_m(indl) + t*vscales_rel(indl)
 
 1669                dst_m(indl-l:indl-1) = beta * dst_m(indl-l:indl-1)
 
 1670                dst_m(indl+1:indl+l) = beta * dst_m(indl+1:indl+l)
 
 1705    integer, 
intent(in) :: p
 
 1706    real(dp), 
intent(in) :: c(3), src_q, dst_r, vscales((p+1)*(p+1)), beta, &
 
 1709    real(dp), 
intent(inout) :: dst_m((p+1)**2)
 
 1711    complex(dp) :: work_complex(max(2, p+1))
 
 1712    real(dp) :: work(p+1)
 
 1714    real(dp) :: dst_sk(p+1), ck(3)
 
 1716    call modified_spherical_bessel_second_kind(p, kappa*dst_r, dst_sk, work, &
 
 1720        & dst_m, work_complex, work)
 
 1750        & work_complex, work)
 
 1753    integer, 
intent(in) :: p
 
 1754    real(dp), 
intent(in) :: c(3), src_q, dst_sk(p+1), vscales((p+1)*(p+1)), beta
 
 1756    real(dp), 
intent(inout) :: dst_m((p+1)**2)
 
 1758    complex(dp), 
intent(out) :: work_complex(p+1)
 
 1759    real(dp), 
intent(out) :: work(p+1)
 
 1761    real(dp) :: rho, ctheta, stheta, cphi, sphi, t, tmp, tmp1, tmp2, &
 
 1762        & max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm, &
 
 1764    integer :: l, m, indl, NZ, ierr
 
 1765    complex(dp) :: complex_argument
 
 1767    if (src_q .eq. zero) 
then 
 1769        if (beta .eq. zero) 
then 
 1773            dst_m = beta * dst_m
 
 1780    if (c(1) .eq. zero) 
then 
 1783    else if (abs(c(2)) .gt. abs(c(1))) 
then 
 1785        ssq12 = one + (c(1)/c(2))**2
 
 1788        ssq12 = one + (c(2)/c(1))**2
 
 1791    if (c(3) .eq. zero) 
then 
 1792        rho = max12 * sqrt(ssq12)
 
 1793    else if (abs(c(3)) .gt. max12) 
then 
 1794        rho = one + ssq12 *(max12/c(3))**2
 
 1795        rho = abs(c(3)) * sqrt(rho)
 
 1797        rho = ssq12 + (c(3)/max12)**2
 
 1798        rho = max12 * sqrt(rho)
 
 1802    if (rho .eq. zero) 
then 
 1806    scaling_factor = src_q * sqrt(two / (pi*rho))
 
 1808    complex_argument = rho
 
 1810    call cbesk(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 1813        call cbesk(complex_argument, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 1815    work = scaling_factor * real(work_complex) / dst_sk
 
 1817    stheta = max12 * sqrt(ssq12)
 
 1819    if (stheta .ne. zero) 
then 
 1821        cphi = c(1) / stheta
 
 1822        sphi = c(2) / stheta
 
 1825        stheta = stheta / rho
 
 1829                if (beta .eq. zero) 
then 
 1830                    dst_m(1) = work(1) * vscales(1)
 
 1832                    dst_m(1) = beta*dst_m(1) + work(1)*vscales(1)
 
 1836                if (beta .eq. zero) 
then 
 1838                    dst_m(1) = work(1) * vscales(1)
 
 1840                    tmp = -vscales(4) * stheta
 
 1841                    tmp2 = work(2) * tmp
 
 1842                    dst_m(2) = tmp2 * sphi
 
 1843                    dst_m(3) = work(2) * vscales(3) * ctheta
 
 1844                    dst_m(4) = tmp2 * cphi
 
 1847                    dst_m(1) = beta*dst_m(1) + work(1)*vscales(1)
 
 1849                    tmp = -vscales(4) * stheta
 
 1850                    tmp2 = work(2) * tmp
 
 1851                    dst_m(2) = beta*dst_m(2) + tmp2*sphi
 
 1852                    dst_m(3) = beta*dst_m(3) + work(2)*vscales(3)*ctheta
 
 1853                    dst_m(4) = beta*dst_m(4) + tmp2*cphi
 
 1859        if (beta .eq. zero) 
then 
 1863            dst_m(1) = work(1) * vscales(1)
 
 1868            ylm = pl1m * vscales(3)
 
 1869            dst_m(3) = work(2) * ylm
 
 1872                plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 1874                ylm = plm * vscales(l*l+l+1)
 
 1875                dst_m(l*l+l+1) = work(l+1) * ylm
 
 1887                ylm = pmm * vscales((m+1)**2)
 
 1888                ylm = work(m+1) * ylm
 
 1889                dst_m((m+1)**2) = cmphi * ylm
 
 1890                dst_m(m*m+1) = smphi * ylm
 
 1892                tmp1 = dble(2*m+1) * pmm
 
 1894                pmm = -stheta * tmp1
 
 1896                pl1m = ctheta * tmp1
 
 1897                ylm = pl1m * vscales((m+1)*(m+3))
 
 1898                ylm = work(m+2) * ylm
 
 1899                dst_m((m+1)*(m+3)) = cmphi * ylm
 
 1900                dst_m((m+1)*(m+2)+1-m) = smphi * ylm
 
 1903                    plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 1904                    plm = plm / dble(l-m)
 
 1905                    ylm = plm * vscales(l*l+l+1+m)
 
 1906                    ylm = work(l+1) * ylm
 
 1907                    dst_m(l*l+l+1+m) = cmphi * ylm
 
 1908                    dst_m(l*l+l+1-m) = smphi * ylm
 
 1915                cmphi = cmphi*cphi - smphi*sphi
 
 1916                smphi = tmp1*sphi + smphi*cphi
 
 1919            ylm = pmm * vscales((p+1)**2)
 
 1920            ylm = work(p+1) * ylm
 
 1921            dst_m((p+1)**2) = cmphi * ylm
 
 1922            dst_m(p*p+1) = smphi * ylm
 
 1928            dst_m(1) = beta*dst_m(1) + work(1)*vscales(1)
 
 1933            ylm = pl1m * vscales(3)
 
 1934            dst_m(3) = beta*dst_m(3) + work(2)*ylm
 
 1937                plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 1939                ylm = plm * vscales(l*l+l+1)
 
 1940                dst_m(l*l+l+1) = beta*dst_m(l*l+l+1) + work(l+1)*ylm
 
 1952                ylm = pmm * vscales((m+1)**2)
 
 1953                ylm = work(m+1) * ylm
 
 1954                dst_m((m+1)**2) = beta*dst_m((m+1)**2) + cmphi*ylm
 
 1955                dst_m(m*m+1) = beta*dst_m(m*m+1) + smphi*ylm
 
 1957                tmp1 = dble(2*m+1) * pmm
 
 1959                pmm = -stheta * tmp1
 
 1961                pl1m = ctheta * tmp1
 
 1962                ylm = pl1m * vscales((m+1)*(m+3))
 
 1963                ylm = work(m+2) * ylm
 
 1964                dst_m((m+1)*(m+3)) = beta*dst_m((m+1)*(m+3)) + cmphi*ylm
 
 1965                dst_m((m+1)*(m+2)+1-m) = beta*dst_m((m+1)*(m+2)+1-m) + &
 
 1969                    plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 1970                    plm = plm / dble(l-m)
 
 1971                    ylm = plm * vscales(l*l+l+1+m)
 
 1972                    ylm = work(l+1) * ylm
 
 1973                    dst_m(l*l+l+1+m) = beta*dst_m(l*l+l+1+m) + cmphi*ylm
 
 1974                    dst_m(l*l+l+1-m) = beta*dst_m(l*l+l+1-m) + smphi*ylm
 
 1981                cmphi = cmphi*cphi - smphi*sphi
 
 1982                smphi = tmp1*sphi + smphi*cphi
 
 1985            ylm = pmm * vscales((p+1)**2)
 
 1986            ylm = work(p+1) * ylm
 
 1987            dst_m((p+1)**2) = beta*dst_m((p+1)**2) + cmphi*ylm
 
 1988            dst_m(p*p+1) = beta*dst_m(p*p+1) + smphi*ylm
 
 1998        if (c(3) .lt. zero) 
then 
 2001            if (beta .eq. zero) 
then 
 2006                    dst_m(indl) = t * work(l+1) * vscales(indl)
 
 2007                    dst_m(indl-l:indl-1) = zero
 
 2008                    dst_m(indl+1:indl+l) = zero
 
 2017                    dst_m(indl) = beta*dst_m(indl) + t*work(l+1)*vscales(indl)
 
 2018                    dst_m(indl-l:indl-1) = beta * dst_m(indl-l:indl-1)
 
 2019                    dst_m(indl+1:indl+l) = beta * dst_m(indl+1:indl+l)
 
 2027            if (beta .eq. zero) 
then 
 2032                    dst_m(indl) = work(l+1) * vscales(indl)
 
 2033                    dst_m(indl-l:indl-1) = zero
 
 2034                    dst_m(indl+1:indl+l) = zero
 
 2041                    dst_m(indl) = beta*dst_m(indl) + work(l+1)*vscales(indl)
 
 2042                    dst_m(indl-l:indl-1) = beta * dst_m(indl-l:indl-1)
 
 2043                    dst_m(indl+1:indl+l) = beta * dst_m(indl+1:indl+l)
 
 2075    integer, 
intent(in) :: p
 
 2076    real(dp), 
intent(in) :: c(3), r, vscales((p+1)**2)
 
 2078    real(dp), 
intent(out) :: mat((p+1)**2)
 
 2080    real(dp) :: rho, ctheta, stheta, cphi, sphi, vcos(p+1), vsin(p+1)
 
 2081    real(dp) :: vylm((p+1)**2), vplm((p+1)**2), rcoef, t
 
 2084    call ylmbas(c, rho, ctheta, stheta, cphi, sphi, p, vscales, vylm, vplm, &
 
 2088    if (rho .eq. zero) 
then 
 2097        mat(ind-n:ind+n) = t / vscales(ind)**2 * vylm(ind-n:ind+n)
 
 2126subroutine fmm_l2p(c, src_r, p, vscales, alpha, src_l, beta, dst_v)
 
 2128    integer, 
intent(in) :: p
 
 2129    real(dp), 
intent(in) :: c(3), src_r, vscales((p+1)*(p+1)), alpha, &
 
 2130        & src_l((p+1)*(p+1)), beta
 
 2132    real(dp), 
intent(inout) :: dst_v
 
 2134    real(dp) :: work(p+1)
 
 2136    call fmm_l2p_work(c, src_r, p, vscales, alpha, src_l, beta, dst_v, work)
 
 2165subroutine fmm_l2p_work(c, src_r, p, vscales_rel, alpha, src_l, beta, dst_v, &
 
 2168    integer, 
intent(in) :: p
 
 2169    real(dp), 
intent(in) :: c(3), src_r, vscales_rel((p+1)*(p+1)), alpha, &
 
 2170        & src_l((p+1)*(p+1)), beta
 
 2172    real(dp), 
intent(inout) :: dst_v
 
 2174    real(dp), 
intent(out), 
target :: work(p+1)
 
 2176    real(dp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
 
 2177        & tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
 
 2178    integer :: l, m, indl
 
 2180    if (beta .eq. zero) 
then 
 2183        dst_v = beta * dst_v
 
 2187    if (alpha .eq. zero) 
then 
 2191    if (c(1) .eq. zero) 
then 
 2194    else if (abs(c(2)) .gt. abs(c(1))) 
then 
 2196        ssq12 = one + (c(1)/c(2))**2
 
 2199        ssq12 = one + (c(2)/c(1))**2
 
 2202    if (c(3) .eq. zero) 
then 
 2203        rho = max12 * sqrt(ssq12)
 
 2204    else if (abs(c(3)) .gt. max12) 
then 
 2205        rho = one + ssq12 *(max12/c(3))**2
 
 2206        rho = abs(c(3)) * sqrt(rho)
 
 2208        rho = ssq12 + (c(3)/max12)**2
 
 2209        rho = max12 * sqrt(rho)
 
 2213    if (rho .eq. zero) 
then 
 2214        dst_v = dst_v + alpha*src_l(1)*vscales_rel(1)
 
 2220    stheta = max12 * sqrt(ssq12)
 
 2222    if (stheta .ne. zero) 
then 
 2224        cphi = c(1) / stheta
 
 2225        sphi = c(2) / stheta
 
 2228        stheta = stheta / rho
 
 2233                dst_v = dst_v + alpha*vscales_rel(1)*src_l(1)
 
 2237                tmp = src_l(1) * vscales_rel(1)
 
 2239                tmp2 = ctheta * src_l(3) * vscales_rel(3)
 
 2240                tmp3 = vscales_rel(4) * stheta
 
 2241                tmp2 = tmp2 - tmp3*sphi*src_l(2)
 
 2242                tmp2 = tmp2 - tmp3*cphi*src_l(4)
 
 2243                dst_v = dst_v + alpha*(tmp+rcoef*tmp2)
 
 2250            work(l+1) = rcoef * work(l)
 
 2255        dst_v = dst_v + work(1)*src_l(1)*vscales_rel(1)
 
 2260        ylm = pl1m * vscales_rel(3)
 
 2261        dst_v = dst_v + work(2)*src_l(3)*ylm
 
 2264            plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 2266            ylm = plm * vscales_rel(l*l+l+1)
 
 2267            dst_v = dst_v + work(l+1)*src_l(l*l+l+1)*ylm
 
 2279            ylm = pmm * vscales_rel((m+1)**2)
 
 2280            tmp1 = cmphi*src_l((m+1)**2) + smphi*src_l(m*m+1)
 
 2281            dst_v = dst_v + work(m+1)*ylm*tmp1
 
 2283            tmp1 = dble(2*m+1) * pmm
 
 2285            pmm = -stheta * tmp1
 
 2287            pl1m = ctheta * tmp1
 
 2288            ylm = pl1m * vscales_rel((m+1)*(m+3))
 
 2289            tmp1 = cmphi*src_l((m+1)*(m+3)) + smphi*src_l((m+1)*(m+2)+1-m)
 
 2290            dst_v = dst_v + work(m+2)*ylm*tmp1
 
 2293                plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 2294                plm = plm / dble(l-m)
 
 2295                ylm = plm * vscales_rel(l*l+l+1+m)
 
 2296                tmp1 = cmphi*src_l(l*l+l+1+m) + smphi*src_l(l*l+l+1-m)
 
 2297                dst_v = dst_v + work(l+1)*ylm*tmp1
 
 2304            cmphi = cmphi*cphi - smphi*sphi
 
 2305            smphi = tmp1*sphi + smphi*cphi
 
 2308        ylm = pmm * vscales_rel((p+1)**2)
 
 2309        tmp1 = cmphi*src_l((p+1)**2) + smphi*src_l(p*p+1)
 
 2310        dst_v = dst_v + work(p+1)*ylm*tmp1
 
 2316        rcoef = sign(rcoef, c(3))
 
 2324            dst_v = dst_v + t*src_l(indl)*vscales_rel(indl)
 
 2358    integer, 
intent(in) :: p
 
 2359    real(dp), 
intent(in) :: c(3), src_r, vscales((p+1)*(p+1)), alpha, &
 
 2360        & src_l((p+1)*(p+1)), beta
 
 2362    real(dp), 
intent(inout) :: dst_v
 
 2364    complex(dp) :: work_complex(max(2, p+1))
 
 2365    real(dp) :: work(p+1)
 
 2367    real(dp) :: src_si(p+1)
 
 2369    call modified_spherical_bessel_first_kind(p, src_r, src_si, work, &
 
 2372        & dst_v, work_complex, work)
 
 2402        & dst_v, work_complex, work)
 
 2405    integer, 
intent(in) :: p
 
 2406    real(dp), 
intent(in) :: c(3), vscales((p+1)*(p+1)), alpha, &
 
 2407        & src_l((p+1)*(p+1)), beta, SI_ri(p+1)
 
 2409    real(dp), 
intent(inout) :: dst_v
 
 2411    complex(dp), 
intent(out) :: work_complex(p+1)
 
 2412    real(dp), 
intent(out) :: work(p+1)
 
 2414    real(dp) :: rho, ctheta, stheta, cphi, sphi, t, tmp, tmp1, tmp2, &
 
 2415        & tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
 
 2416    real(dp) :: scaling_factor
 
 2417    complex(dp) :: complex_argument
 
 2418    integer :: l, m, indl, ierr, NZ
 
 2419    real(dp), 
parameter :: fnu=1.5d0
 
 2421    if (beta .eq. zero) 
then 
 2424        dst_v = beta * dst_v
 
 2428    if (alpha .eq. zero) 
then 
 2432    if (c(1) .eq. zero) 
then 
 2435    else if (abs(c(2)) .gt. abs(c(1))) 
then 
 2437        ssq12 = one + (c(1)/c(2))**2
 
 2440        ssq12 = one + (c(2)/c(1))**2
 
 2443    if (c(3) .eq. zero) 
then 
 2444        rho = max12 * sqrt(ssq12)
 
 2445    else if (abs(c(3)) .gt. max12) 
then 
 2446        rho = one + ssq12 *(max12/c(3))**2
 
 2447        rho = abs(c(3)) * sqrt(rho)
 
 2449        rho = ssq12 + (c(3)/max12)**2
 
 2450        rho = max12 * sqrt(rho)
 
 2454    if (rho .eq. zero) 
then 
 2455        dst_v = dst_v + alpha*src_l(1)*vscales(1)/si_ri(1)
 
 2459    scaling_factor = alpha * sqrt(pi / (2*rho))
 
 2461    complex_argument = rho
 
 2463    call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 2466        call cbesi(complex_argument, fnu, 1, p, work_complex(2:p+1), nz, ierr)
 
 2468    work = scaling_factor * real(work_complex) / si_ri
 
 2471    stheta = max12 * sqrt(ssq12)
 
 2473    if (stheta .ne. zero) 
then 
 2475        cphi = c(1) / stheta
 
 2476        sphi = c(2) / stheta
 
 2479        stheta = stheta / rho
 
 2484                dst_v = dst_v + work(1)*vscales(1)*src_l(1)
 
 2488                tmp = work(1) * src_l(1) * vscales(1)
 
 2490                tmp2 = ctheta * src_l(3) * vscales(3)
 
 2491                tmp3 = stheta * vscales(4)
 
 2492                tmp2 = tmp2 - tmp3*sphi*src_l(2)
 
 2493                tmp2 = tmp2 - tmp3*cphi*src_l(4)
 
 2494                dst_v = dst_v + tmp + work(2)*tmp2
 
 2501        dst_v = dst_v + work(1)*src_l(1)*vscales(1)
 
 2506        ylm = pl1m * vscales(3)
 
 2507        dst_v = dst_v + work(2)*src_l(3)*ylm
 
 2510            plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 2512            ylm = plm * vscales(l*l+l+1)
 
 2513            dst_v = dst_v + work(l+1)*src_l(l*l+l+1)*ylm
 
 2525            ylm = pmm * vscales((m+1)**2)
 
 2526            tmp1 = cmphi*src_l((m+1)**2) + smphi*src_l(m*m+1)
 
 2527            dst_v = dst_v + work(m+1)*ylm*tmp1
 
 2529            tmp1 = dble(2*m+1) * pmm
 
 2531            pmm = -stheta * tmp1
 
 2533            pl1m = ctheta * tmp1
 
 2534            ylm = pl1m * vscales((m+1)*(m+3))
 
 2535            tmp1 = cmphi*src_l((m+1)*(m+3)) + smphi*src_l((m+1)*(m+2)+1-m)
 
 2536            dst_v = dst_v + work(m+2)*ylm*tmp1
 
 2539                plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 2540                plm = plm / dble(l-m)
 
 2541                ylm = plm * vscales(l*l+l+1+m)
 
 2542                tmp1 = cmphi*src_l(l*l+l+1+m) + smphi*src_l(l*l+l+1-m)
 
 2543                dst_v = dst_v + work(l+1)*ylm*tmp1
 
 2550            cmphi = cmphi*cphi - smphi*sphi
 
 2551            smphi = tmp1*sphi + smphi*cphi
 
 2554        ylm = pmm * vscales((p+1)**2)
 
 2555        tmp1 = cmphi*src_l((p+1)**2) + smphi*src_l(p*p+1)
 
 2556        dst_v = dst_v + work(p+1)*ylm*tmp1
 
 2565        if (c(3) .lt. zero) 
then 
 2572                dst_v = dst_v + t*work(l+1)*src_l(indl)*vscales(indl)
 
 2582                dst_v = dst_v + work(l+1)*src_l(indl)*vscales(indl)
 
 2588subroutine fmm_l2p_bessel_grad(c, src_r, p, vscales, alpha, src_l, beta, dst_g)
 
 2591    integer, 
intent(in) :: p
 
 2592    real(dp), 
intent(in) :: c(3), src_r, vscales((p+2)*(p+2)), alpha, &
 
 2593        & src_l((p+1)*(p+1)), beta
 
 2595    real(dp), 
intent(inout) :: dst_g(3)
 
 2597    complex(dp) :: work_complex(p+2)
 
 2598    real(dp) :: work(p+2), src_si(p+2), src_l_grad((p+2)**2, 3)
 
 2600    call modified_spherical_bessel_first_kind(p+1, src_r, src_si, work, &
 
 2602    call fmm_l2l_bessel_grad(p, src_si, vscales, src_l, src_l_grad)
 
 2604        & beta, dst_g(1), work_complex, work)
 
 2606        & beta, dst_g(2), work_complex, work)
 
 2608        & beta, dst_g(3), work_complex, work)
 
 2609end subroutine fmm_l2p_bessel_grad
 
 2637    integer, 
intent(in) :: p
 
 2638    real(dp), 
intent(in) :: c(3), src_q, dst_r, vscales_rel((p+1)*(p+1)), beta
 
 2640    real(dp), 
intent(inout) :: dst_l((p+1)**2)
 
 2642    real(dp) :: work(p+1)
 
 2674    integer, 
intent(in) :: p
 
 2675    real(dp), 
intent(in) :: c(3), src_q, dst_r, vscales_rel((p+1)*(p+1)), beta
 
 2677    real(dp), 
intent(inout) :: dst_l((p+1)**2)
 
 2679    real(dp), 
intent(out), 
target :: work(p+1)
 
 2681    real(dp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
 
 2682        & max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
 
 2683    integer :: l, m, indl
 
 2685    if (src_q .eq. zero) 
then 
 2687        if (beta .eq. zero) 
then 
 2691            dst_l = beta * dst_l
 
 2698    if (c(1) .eq. zero) 
then 
 2701    else if (abs(c(2)) .gt. abs(c(1))) 
then 
 2703        ssq12 = one + (c(1)/c(2))**2
 
 2706        ssq12 = one + (c(2)/c(1))**2
 
 2709    if (c(3) .eq. zero) 
then 
 2710        rho = max12 * sqrt(ssq12)
 
 2711    else if (abs(c(3)) .gt. max12) 
then 
 2712        rho = one + ssq12 *(max12/c(3))**2
 
 2713        rho = abs(c(3)) * sqrt(rho)
 
 2715        rho = ssq12 + (c(3)/max12)**2
 
 2716        rho = max12 * sqrt(rho)
 
 2720    if (rho .eq. zero) 
then 
 2721        if (beta .eq. zero) 
then 
 2722            dst_l(1) = src_q * vscales_rel(1)
 
 2725            dst_l(1) = beta*dst_l(1) + src_q*vscales_rel(1)
 
 2726            dst_l(2:) = beta * dst_l(2:)
 
 2733    stheta = max12 * sqrt(ssq12)
 
 2735    if (stheta .ne. zero) 
then 
 2737        cphi = c(1) / stheta
 
 2738        sphi = c(2) / stheta
 
 2741        stheta = stheta / rho
 
 2746                if (beta .eq. zero) 
then 
 2747                    dst_l(1) = src_q * vscales_rel(1)
 
 2749                    dst_l(1) = beta*dst_l(1) + src_q*vscales_rel(1)
 
 2755                tmp2 = -tmp * vscales_rel(4) * stheta
 
 2756                if (beta .eq. zero) 
then 
 2757                    dst_l(1) = src_q * vscales_rel(1)
 
 2758                    dst_l(2) = tmp2 * sphi
 
 2759                    dst_l(3) = tmp * vscales_rel(3) * ctheta
 
 2760                    dst_l(4) = tmp2 * cphi
 
 2762                    dst_l(1) = beta*dst_l(1) + src_q*vscales_rel(1)
 
 2763                    dst_l(2) = beta*dst_l(2) + tmp2*sphi
 
 2764                    dst_l(3) = beta*dst_l(3) + tmp*vscales_rel(3)*ctheta
 
 2765                    dst_l(4) = beta*dst_l(4) + tmp2*cphi
 
 2773            work(l+1) = rcoef * work(l)
 
 2776        if (beta .eq. zero) 
then 
 2780            dst_l(1) = work(1) * vscales_rel(1)
 
 2785            ylm = pl1m * vscales_rel(3)
 
 2786            dst_l(3) = work(2) * ylm
 
 2789                plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 2791                ylm = plm * vscales_rel(l*l+l+1)
 
 2792                dst_l(l*l+l+1) = work(l+1) * ylm
 
 2804                ylm = pmm * vscales_rel((m+1)**2)
 
 2805                ylm = work(m+1) * ylm
 
 2806                dst_l((m+1)**2) = cmphi * ylm
 
 2807                dst_l(m*m+1) = smphi * ylm
 
 2809                tmp1 = dble(2*m+1) * pmm
 
 2811                pmm = -stheta * tmp1
 
 2813                pl1m = ctheta * tmp1
 
 2814                ylm = pl1m * vscales_rel((m+1)*(m+3))
 
 2815                ylm = work(m+2) * ylm
 
 2816                dst_l((m+1)*(m+3)) = cmphi * ylm
 
 2817                dst_l((m+1)*(m+2)+1-m) = smphi * ylm
 
 2820                    plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 2821                    plm = plm / dble(l-m)
 
 2822                    ylm = plm * vscales_rel(l*l+l+1+m)
 
 2823                    ylm = work(l+1) * ylm
 
 2824                    dst_l(l*l+l+1+m) = cmphi * ylm
 
 2825                    dst_l(l*l+l+1-m) = smphi * ylm
 
 2832                cmphi = cmphi*cphi - smphi*sphi
 
 2833                smphi = tmp1*sphi + smphi*cphi
 
 2836            ylm = pmm * vscales_rel((p+1)**2)
 
 2837            ylm = work(p+1) * ylm
 
 2838            dst_l((p+1)**2) = cmphi * ylm
 
 2839            dst_l(p*p+1) = smphi * ylm
 
 2845            dst_l(1) = beta*dst_l(1) + work(1)*vscales_rel(1)
 
 2850            ylm = pl1m * vscales_rel(3)
 
 2851            dst_l(3) = beta*dst_l(3) + work(2)*ylm
 
 2854                plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
 
 2856                ylm = plm * vscales_rel(l*l+l+1)
 
 2857                dst_l(l*l+l+1) = beta*dst_l(l*l+l+1) + work(l+1)*ylm
 
 2869                ylm = pmm * vscales_rel((m+1)**2)
 
 2870                ylm = work(m+1) * ylm
 
 2871                dst_l((m+1)**2) = beta*dst_l((m+1)**2) + cmphi*ylm
 
 2872                dst_l(m*m+1) = beta*dst_l(m*m+1) + smphi*ylm
 
 2874                tmp1 = dble(2*m+1) * pmm
 
 2876                pmm = -stheta * tmp1
 
 2878                pl1m = ctheta * tmp1
 
 2879                ylm = pl1m * vscales_rel((m+1)*(m+3))
 
 2880                ylm = work(m+2) * ylm
 
 2881                dst_l((m+1)*(m+3)) = beta*dst_l((m+1)*(m+3)) + cmphi*ylm
 
 2882                dst_l((m+1)*(m+2)+1-m) = beta*dst_l((m+1)*(m+2)+1-m) + &
 
 2886                    plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
 
 2887                    plm = plm / dble(l-m)
 
 2888                    ylm = plm * vscales_rel(l*l+l+1+m)
 
 2889                    ylm = work(l+1) * ylm
 
 2890                    dst_l(l*l+l+1+m) = beta*dst_l(l*l+l+1+m) + cmphi*ylm
 
 2891                    dst_l(l*l+l+1-m) = beta*dst_l(l*l+l+1-m) + smphi*ylm
 
 2898                cmphi = cmphi*cphi - smphi*sphi
 
 2899                smphi = tmp1*sphi + smphi*cphi
 
 2902            ylm = pmm * vscales_rel((p+1)**2)
 
 2903            ylm = work(p+1) * ylm
 
 2904            dst_l((p+1)**2) = beta*dst_l((p+1)**2) + cmphi*ylm
 
 2905            dst_l(p*p+1) = beta*dst_l(p*p+1) + smphi*ylm
 
 2912        rcoef = sign(rcoef, c(3))
 
 2917        if (beta .eq. zero) 
then 
 2922                dst_l(indl) = t * vscales_rel(indl)
 
 2923                dst_l(indl-l:indl-1) = zero
 
 2924                dst_l(indl+1:indl+l) = zero
 
 2934                dst_l(indl) = beta*dst_l(indl) + t*vscales_rel(indl)
 
 2935                dst_l(indl-l:indl-1) = beta * dst_l(indl-l:indl-1)
 
 2936                dst_l(indl+1:indl+l) = beta * dst_l(indl+1:indl+l)
 
 3000    integer, 
intent(in) :: p
 
 3001    real(dp), 
intent(in) :: r1(-1:1, -1:1), alpha, src((p+1)*(p+1)), beta
 
 3003    real(dp), 
intent(inout) :: dst((p+1)*(p+1))
 
 3005    real(dp) :: work(2*(2*p+1)*(2*p+3))
 
 3025    integer, 
intent(in) :: p
 
 3026    real(dp), 
intent(in) :: r1(-1:1, -1:1), alpha, src((p+1)*(p+1)), beta
 
 3028    real(dp), 
intent(out) :: dst((p+1)*(p+1))
 
 3030    real(dp), 
intent(out), 
target :: work(2*(2*p+1)*(2*p+3))
 
 3033    integer :: l, m, n, ind
 
 3035    real(dp), 
pointer :: r_prev(:, :), r(:, :), scal_uvw_m(:), scal_u_n(:), &
 
 3036        & scal_v_n(:), scal_w_n(:), r_swap(:, :)
 
 3038    if (alpha .eq. zero) 
then 
 3040        if (beta .eq. zero) 
then 
 3050    if (beta .eq. zero) 
then 
 3053        dst(1) = alpha * src(1)
 
 3058        dst(2) = alpha * dot_product(r1(-1:1,-1), src(2:4))
 
 3059        dst(3) = alpha * dot_product(r1(-1:1,0), src(2:4))
 
 3060        dst(4) = alpha * dot_product(r1(-1:1,1), src(2:4))
 
 3067        r_prev(-p:p, -p:p) => work(1:l)
 
 3069        r(-p:p, -p:p) => work(l+1:m)
 
 3071        scal_uvw_m(-p:p) => work(m+1:l)
 
 3073        scal_u_n(-p:p) => work(l+1:m)
 
 3075        scal_v_n(-p:p) => work(m+1:l)
 
 3077        scal_w_n(-p:p) => work(l+1:m)
 
 3079        r(2, 2) = (r1(1, 1)*r1(1, 1) - r1(1, -1)*r1(1, -1) - &
 
 3080            & r1(-1, 1)*r1(-1, 1) + r1(-1, -1)*r1(-1, -1)) / two
 
 3081        r(2, 1) = r1(1, 1)*r1(1, 0) - r1(-1, 1)*r1(-1, 0)
 
 3082        r(2, 0) = sqrt3 / two * (r1(1, 0)*r1(1, 0) - r1(-1, 0)*r1(-1, 0))
 
 3083        r(2, -1) = r1(1, -1)*r1(1, 0) - r1(-1, -1)*r1(-1, 0)
 
 3084        r(2, -2) = r1(1, 1)*r1(1, -1) - r1(-1, 1)*r1(-1, -1)
 
 3085        r(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
 
 3086        r(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
 
 3087        r(1, 0) = sqrt3 * r1(1, 0) * r1(0, 0)
 
 3088        r(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
 
 3089        r(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
 
 3090        r(0, 2) = sqrt3 / two * (r1(0, 1)*r1(0, 1) - r1(0, -1)*r1(0, -1))
 
 3091        r(0, 1) = sqrt3 * r1(0, 1) * r1(0, 0)
 
 3092        r(0, 0) = (three*r1(0, 0)*r1(0, 0)-one) / two
 
 3093        r(0, -1) = sqrt3 * r1(0, -1) * r1(0, 0)
 
 3094        r(0, -2) = sqrt3 * r1(0, 1) * r1(0, -1)
 
 3095        r(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
 
 3096        r(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
 
 3097        r(-1, 0) = sqrt3 * r1(-1, 0) * r1(0, 0)
 
 3098        r(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(0, -1)*r1(-1, 0)
 
 3099        r(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
 
 3100        r(-2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
 
 3101        r(-2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
 
 3102        r(-2, 0) = sqrt3 * r1(1, 0) * r1(-1, 0)
 
 3103        r(-2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
 
 3104        r(-2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
 
 3106            dst(7+m) = alpha * dot_product(r(-2:2, m), src(5:9))
 
 3115            scal_uvw_m(0) = dble(l)
 
 3117                u = sqrt(dble(l*l-m*m))
 
 3122            u = sqrt(dble(u*(u-one)))
 
 3125            scal_u_n(0) = dble(l)
 
 3126            scal_v_n(0) = -sqrt(dble(l*(l-1))) / sqrt2
 
 3129                u = sqrt(dble(l*l-n*n))
 
 3133                v = sqrt(v*(v-one)) / two
 
 3137                w = -sqrt(w*(w-one)) / two
 
 3141            u = sqrt(dble(2*l-1))
 
 3146            v = sqrt(dble((2*l-1)*(2*l-2))) / two
 
 3149            v = sqrt(dble(2*l*(2*l-1))) / two
 
 3152            scal_w_n(l-1) = zero
 
 3155            scal_w_n(1-l) = zero
 
 3158            v = r1(1, 1)*r_prev(l-1, l-1) - r1(1, -1)*r_prev(l-1, 1-l) - &
 
 3159                & r1(-1, 1)*r_prev(1-l, l-1) + r1(-1, -1)*r_prev(1-l, 1-l)
 
 3160            r(l, l) = v * scal_v_n(l) / scal_uvw_m(l)
 
 3163            v = r1(1, 1)*r_prev(1-l, l-1) - r1(1, -1)*r_prev(1-l, 1-l) + &
 
 3164                & r1(-1, 1)*r_prev(l-1, l-1) - r1(-1, -1)*r_prev(l-1, 1-l)
 
 3165            r(-l, l) = v * scal_v_n(l) / scal_uvw_m(l)
 
 3168            u = r1(0, 1)*r_prev(l-1, l-1) - r1(0, -1)*r_prev(l-1, 1-l)
 
 3169            v = r1(1, 1)*r_prev(l-2, l-1) - r1(1, -1)*r_prev(l-2, 1-l) - &
 
 3170                & r1(-1, 1)*r_prev(2-l, l-1) + r1(-1, -1)*r_prev(2-l, 1-l)
 
 3171            r(l-1, l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
 
 3174            u = r1(0, 1)*r_prev(1-l, l-1) - r1(0, -1)*r_prev(1-l, 1-l)
 
 3175            v = r1(1, 1)*r_prev(2-l, l-1) - r1(1, -1)*r_prev(2-l, 1-l) + &
 
 3176                & r1(-1, 1)*r_prev(l-2, l-1) - r1(-1, -1)*r_prev(l-2, 1-l)
 
 3177            r(1-l, l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
 
 3180            u = r1(0, 1)*r_prev(1, l-1) - r1(0, -1)*r_prev(1, 1-l)
 
 3181            v = r1(1, 1)*r_prev(0, l-1) - r1(1, -1)*r_prev(0, 1-l)
 
 3182            w = r1(1, 1)*r_prev(2, l-1) - r1(1, -1)*r_prev(2, 1-l) + &
 
 3183                & r1(-1, 1)*r_prev(-2, l-1) - r1(-1, -1)*r_prev(-2, 1-l)
 
 3184            r(1, l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3185            r(1, l) = r(1, l) / scal_uvw_m(l)
 
 3188            u = r1(0, 1)*r_prev(-1, l-1) - r1(0, -1)*r_prev(-1, 1-l)
 
 3189            v = r1(-1, 1)*r_prev(0, l-1) - r1(-1, -1)*r_prev(0, 1-l)
 
 3190            w = r1(1, 1)*r_prev(-2, l-1) - r1(1, -1)*r_prev(-2, 1-l) - &
 
 3191                & r1(-1, 1)*r_prev(2, l-1) + r1(-1, -1)*r_prev(2, 1-l)
 
 3192            r(-1, l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3193            r(-1, l) = r(-1, l) / scal_uvw_m(l)
 
 3196            u = r1(0, 1)*r_prev(0, l-1) - r1(0, -1)*r_prev(0, 1-l)
 
 3197            v = r1(1, 1)*r_prev(1, l-1) - r1(1, -1)*r_prev(1, 1-l) + &
 
 3198                & r1(-1, 1)*r_prev(-1, l-1) - r1(-1, -1)*r_prev(-1, 1-l)
 
 3199            r(0, l) = (u*scal_u_n(0) + v*scal_v_n(0)) / scal_uvw_m(l)
 
 3203                u = r1(0, 1)*r_prev(n, l-1) - r1(0, -1)*r_prev(n, 1-l)
 
 3204                v = r1(1, 1)*r_prev(n-1, l-1) - r1(1, -1)*r_prev(n-1, 1-l) - &
 
 3205                    & r1(-1, 1)*r_prev(1-n, l-1) + r1(-1, -1)*r_prev(1-n, 1-l)
 
 3206                w = r1(1, 1)*r_prev(n+1, l-1) - r1(1, -1)*r_prev(n+1, 1-l) + &
 
 3207                    & r1(-1, 1)*r_prev(-n-1, l-1) - r1(-1, -1)*r_prev(-n-1, 1-l)
 
 3208                r(n, l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3209                r(n, l) = r(n, l) / scal_uvw_m(l)
 
 3211                u = r1(0, 1)*r_prev(-n, l-1) - r1(0, -1)*r_prev(-n, 1-l)
 
 3212                v = r1(1, 1)*r_prev(1-n, l-1) - r1(1, -1)*r_prev(1-n, 1-l) + &
 
 3213                    & r1(-1, 1)*r_prev(n-1, l-1) - r1(-1, -1)*r_prev(n-1, 1-l)
 
 3214                w = r1(1, 1)*r_prev(-n-1, l-1) - r1(1, -1)*r_prev(-n-1, 1-l) - &
 
 3215                    & r1(-1, 1)*r_prev(n+1, l-1) + r1(-1, -1)*r_prev(n+1, 1-l)
 
 3216                r(-n, l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3217                r(-n, l) = r(-n, l) / scal_uvw_m(l)
 
 3220            dst(ind+l) = alpha * dot_product(src(ind-l:ind+l), r(-l:l, l))
 
 3222            v = r1(1, 1)*r_prev(l-1, 1-l) + r1(1, -1)*r_prev(l-1, l-1) - &
 
 3223                & r1(-1, 1)*r_prev(1-l, 1-l) - r1(-1, -1)*r_prev(1-l, l-1)
 
 3224            r(l, -l) = v * scal_v_n(l) / scal_uvw_m(l)
 
 3227            v = r1(1, 1)*r_prev(1-l, 1-l) + r1(1, -1)*r_prev(1-l, l-1) + &
 
 3228                & r1(-1, 1)*r_prev(l-1, 1-l) + r1(-1, -1)*r_prev(l-1, l-1)
 
 3229            r(-l, -l) = v * scal_v_n(l) / scal_uvw_m(l)
 
 3232            u = r1(0, 1)*r_prev(l-1, 1-l) + r1(0, -1)*r_prev(l-1, l-1)
 
 3233            v = r1(1, 1)*r_prev(l-2, 1-l) + r1(1, -1)*r_prev(l-2, l-1) - &
 
 3234                & r1(-1, 1)*r_prev(2-l, 1-l) - r1(-1, -1)*r_prev(2-l, l-1)
 
 3235            r(l-1, -l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
 
 3238            u = r1(0, 1)*r_prev(1-l, 1-l) + r1(0, -1)*r_prev(1-l, l-1)
 
 3239            v = r1(1, 1)*r_prev(2-l, 1-l) + r1(1, -1)*r_prev(2-l, l-1) + &
 
 3240                & r1(-1, 1)*r_prev(l-2, 1-l) + r1(-1, -1)*r_prev(l-2, l-1)
 
 3241            r(1-l, -l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
 
 3244            u = r1(0, 1)*r_prev(1, 1-l) + r1(0, -1)*r_prev(1, l-1)
 
 3245            v = r1(1, 1)*r_prev(0, 1-l) + r1(1, -1)*r_prev(0, l-1)
 
 3246            w = r1(1, 1)*r_prev(2, 1-l) + r1(1, -1)*r_prev(2, l-1) + &
 
 3247                & r1(-1, 1)*r_prev(-2, 1-l) + r1(-1, -1)*r_prev(-2, l-1)
 
 3248            r(1, -l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3249            r(1, -l) = r(1, -l) / scal_uvw_m(l)
 
 3252            u = r1(0, 1)*r_prev(-1, 1-l) + r1(0, -1)*r_prev(-1, l-1)
 
 3253            v = r1(-1, 1)*r_prev(0, 1-l) + r1(-1, -1)*r_prev(0, l-1)
 
 3254            w = r1(1, 1)*r_prev(-2, 1-l) + r1(1, -1)*r_prev(-2, l-1) - &
 
 3255                & r1(-1, 1)*r_prev(2, 1-l) - r1(-1, -1)*r_prev(2, l-1)
 
 3256            r(-1, -l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3257            r(-1, -l) = r(-1, -l) / scal_uvw_m(l)
 
 3260            u = r1(0, 1)*r_prev(0, 1-l) + r1(0, -1)*r_prev(0, l-1)
 
 3261            v = r1(1, 1)*r_prev(1, 1-l) + r1(1, -1)*r_prev(1, l-1) + &
 
 3262                & r1(-1, 1)*r_prev(-1, 1-l) + r1(-1, -1)*r_prev(-1, l-1)
 
 3263            r(0, -l) = (u*scal_u_n(0) + v*scal_v_n(0)) / scal_uvw_m(l)
 
 3267                u = r1(0, 1)*r_prev(n, 1-l) + r1(0, -1)*r_prev(n, l-1)
 
 3268                v = r1(1, 1)*r_prev(n-1, 1-l) + r1(1, -1)*r_prev(n-1, l-1) - &
 
 3269                    & r1(-1, 1)*r_prev(1-n, 1-l) - r1(-1, -1)*r_prev(1-n, l-1)
 
 3270                w = r1(1, 1)*r_prev(n+1, 1-l) + r1(1, -1)*r_prev(n+1, l-1) + &
 
 3271                    & r1(-1, 1)*r_prev(-n-1, 1-l) + r1(-1, -1)*r_prev(-n-1, l-1)
 
 3272                r(n, -l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3273                r(n, -l) = r(n, -l) / scal_uvw_m(l)
 
 3275                u = r1(0, 1)*r_prev(-n, 1-l) + r1(0, -1)*r_prev(-n, l-1)
 
 3276                v = r1(1, 1)*r_prev(1-n, 1-l) + r1(1, -1)*r_prev(1-n, l-1) + &
 
 3277                    & r1(-1, 1)*r_prev(n-1, 1-l) + r1(-1, -1)*r_prev(n-1, l-1)
 
 3278                w = r1(1, 1)*r_prev(-n-1, 1-l) + r1(1, -1)*r_prev(-n-1, l-1) - &
 
 3279                    & r1(-1, 1)*r_prev(n+1, 1-l) - r1(-1, -1)*r_prev(n+1, l-1)
 
 3280                r(-n, -l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3281                r(-n, -l) = r(-n, -l) / scal_uvw_m(l)
 
 3284            dst(ind-l) = alpha * dot_product(src(ind-l:ind+l), r(-l:l, -l))
 
 3288                v = r1(1, 0)*r_prev(l-1, m) - r1(-1, 0)*r_prev(1-l, m)
 
 3289                r(l, m) = v * scal_v_n(l) / scal_uvw_m(m)
 
 3292                v = r1(1, 0)*r_prev(1-l, m) + r1(-1, 0)*r_prev(l-1, m)
 
 3293                r(-l, m) = v * scal_v_n(l) / scal_uvw_m(m)
 
 3296                u = r1(0, 0) * r_prev(l-1, m)
 
 3297                v = r1(1, 0)*r_prev(l-2, m) - r1(-1, 0)*r_prev(2-l, m)
 
 3298                r(l-1, m) = u*scal_u_n(l-1) + v*scal_v_n(l-1)
 
 3299                r(l-1, m) = r(l-1, m) / scal_uvw_m(m)
 
 3302                u = r1(0, 0) * r_prev(1-l, m)
 
 3303                v = r1(1, 0)*r_prev(2-l, m) + r1(-1, 0)*r_prev(l-2, m)
 
 3304                r(1-l, m) = u*scal_u_n(l-1) + v*scal_v_n(l-1)
 
 3305                r(1-l, m) = r(1-l, m) / scal_uvw_m(m)
 
 3308                u = r1(0, 0) * r_prev(0, m)
 
 3309                v = r1(1, 0)*r_prev(1, m) + r1(-1, 0)*r_prev(-1, m)
 
 3310                r(0, m) = u*scal_u_n(0) + v*scal_v_n(0)
 
 3311                r(0, m) = r(0, m) / scal_uvw_m(m)
 
 3314                u = r1(0, 0) * r_prev(1, m)
 
 3315                v = r1(1, 0) * r_prev(0, m)
 
 3316                w = r1(1, 0)*r_prev(2, m) + r1(-1, 0)*r_prev(-2, m)
 
 3317                r(1, m) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3318                r(1, m) = r(1, m) / scal_uvw_m(m)
 
 3321                u = r1(0, 0) * r_prev(-1, m)
 
 3322                v = r1(-1, 0) * r_prev(0, m)
 
 3323                w = r1(1, 0)*r_prev(-2, m) - r1(-1, 0)*r_prev(2, m)
 
 3324                r(-1, m) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3325                r(-1, m) = r(-1, m) / scal_uvw_m(m)
 
 3329                    u = r1(0, 0) * r_prev(n, m)
 
 3330                    v = r1(1, 0)*r_prev(n-1, m) - r1(-1, 0)*r_prev(1-n, m)
 
 3331                    w = r1(1, 0)*r_prev(n+1, m) + r1(-1, 0)*r_prev(-1-n, m)
 
 3332                    r(n, m) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3333                    r(n, m) = r(n, m) / scal_uvw_m(m)
 
 3335                    u = r1(0, 0) * r_prev(-n, m)
 
 3336                    v = r1(1, 0)*r_prev(1-n, m) + r1(-1, 0)*r_prev(n-1, m)
 
 3337                    w = r1(1, 0)*r_prev(-n-1, m) - r1(-1, 0)*r_prev(n+1, m)
 
 3338                    r(-n, m) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3339                    r(-n, m) = r(-n, m) / scal_uvw_m(m)
 
 3342                dst(ind+m) = alpha * dot_product(src(ind-l:ind+l), r(-l:l, m))
 
 3348        dst(1) = beta*dst(1) + alpha*src(1)
 
 3353        dst(2) = beta*dst(2) + alpha*dot_product(r1(-1:1,-1), src(2:4))
 
 3354        dst(3) = beta*dst(3) + alpha*dot_product(r1(-1:1,0), src(2:4))
 
 3355        dst(4) = beta*dst(4) + alpha*dot_product(r1(-1:1,1), src(2:4))
 
 3362        r_prev(-p:p, -p:p) => work(1:l)
 
 3364        r(-p:p, -p:p) => work(l+1:m)
 
 3366        scal_uvw_m(-p:p) => work(m+1:l)
 
 3368        scal_u_n(-p:p) => work(l+1:m)
 
 3370        scal_v_n(-p:p) => work(m+1:l)
 
 3372        scal_w_n(-p:p) => work(l+1:m)
 
 3374        r(2, 2) = (r1(1, 1)*r1(1, 1) - r1(1, -1)*r1(1, -1) - &
 
 3375            & r1(-1, 1)*r1(-1, 1) + r1(-1, -1)*r1(-1, -1)) / two
 
 3376        r(2, 1) = r1(1, 1)*r1(1, 0) - r1(-1, 1)*r1(-1, 0)
 
 3377        r(2, 0) = sqrt3 / two * (r1(1, 0)*r1(1, 0) - r1(-1, 0)*r1(-1, 0))
 
 3378        r(2, -1) = r1(1, -1)*r1(1, 0) - r1(-1, -1)*r1(-1, 0)
 
 3379        r(2, -2) = r1(1, 1)*r1(1, -1) - r1(-1, 1)*r1(-1, -1)
 
 3380        r(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
 
 3381        r(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
 
 3382        r(1, 0) = sqrt3 * r1(1, 0) * r1(0, 0)
 
 3383        r(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
 
 3384        r(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
 
 3385        r(0, 2) = sqrt3 / two * (r1(0, 1)*r1(0, 1) - r1(0, -1)*r1(0, -1))
 
 3386        r(0, 1) = sqrt3 * r1(0, 1) * r1(0, 0)
 
 3387        r(0, 0) = (three*r1(0, 0)*r1(0, 0)-one) / two
 
 3388        r(0, -1) = sqrt3 * r1(0, -1) * r1(0, 0)
 
 3389        r(0, -2) = sqrt3 * r1(0, 1) * r1(0, -1)
 
 3390        r(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
 
 3391        r(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
 
 3392        r(-1, 0) = sqrt3 * r1(-1, 0) * r1(0, 0)
 
 3393        r(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(0, -1)*r1(-1, 0)
 
 3394        r(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
 
 3395        r(-2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
 
 3396        r(-2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
 
 3397        r(-2, 0) = sqrt3 * r1(1, 0) * r1(-1, 0)
 
 3398        r(-2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
 
 3399        r(-2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
 
 3401            dst(7+m) = beta*dst(7+m) + alpha*dot_product(r(-2:2, m), src(5:9))
 
 3410            scal_uvw_m(0) = dble(l)
 
 3412                u = sqrt(dble(l*l-m*m))
 
 3417            u = sqrt(dble(u*(u-one)))
 
 3420            scal_u_n(0) = dble(l)
 
 3421            scal_v_n(0) = -sqrt(dble(l*(l-1))) / sqrt2
 
 3424                u = sqrt(dble(l*l-n*n))
 
 3428                v = sqrt(v*(v-one)) / two
 
 3432                w = -sqrt(w*(w-one)) / two
 
 3436            u = sqrt(dble(2*l-1))
 
 3441            v = sqrt(dble((2*l-1)*(2*l-2))) / two
 
 3444            v = sqrt(dble(2*l*(2*l-1))) / two
 
 3447            scal_w_n(l-1) = zero
 
 3450            scal_w_n(1-l) = zero
 
 3453            v = r1(1, 1)*r_prev(l-1, l-1) - r1(1, -1)*r_prev(l-1, 1-l) - &
 
 3454                & r1(-1, 1)*r_prev(1-l, l-1) + r1(-1, -1)*r_prev(1-l, 1-l)
 
 3455            r(l, l) = v * scal_v_n(l) / scal_uvw_m(l)
 
 3458            v = r1(1, 1)*r_prev(1-l, l-1) - r1(1, -1)*r_prev(1-l, 1-l) + &
 
 3459                & r1(-1, 1)*r_prev(l-1, l-1) - r1(-1, -1)*r_prev(l-1, 1-l)
 
 3460            r(-l, l) = v * scal_v_n(l) / scal_uvw_m(l)
 
 3463            u = r1(0, 1)*r_prev(l-1, l-1) - r1(0, -1)*r_prev(l-1, 1-l)
 
 3464            v = r1(1, 1)*r_prev(l-2, l-1) - r1(1, -1)*r_prev(l-2, 1-l) - &
 
 3465                & r1(-1, 1)*r_prev(2-l, l-1) + r1(-1, -1)*r_prev(2-l, 1-l)
 
 3466            r(l-1, l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
 
 3469            u = r1(0, 1)*r_prev(1-l, l-1) - r1(0, -1)*r_prev(1-l, 1-l)
 
 3470            v = r1(1, 1)*r_prev(2-l, l-1) - r1(1, -1)*r_prev(2-l, 1-l) + &
 
 3471                & r1(-1, 1)*r_prev(l-2, l-1) - r1(-1, -1)*r_prev(l-2, 1-l)
 
 3472            r(1-l, l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
 
 3475            u = r1(0, 1)*r_prev(1, l-1) - r1(0, -1)*r_prev(1, 1-l)
 
 3476            v = r1(1, 1)*r_prev(0, l-1) - r1(1, -1)*r_prev(0, 1-l)
 
 3477            w = r1(1, 1)*r_prev(2, l-1) - r1(1, -1)*r_prev(2, 1-l) + &
 
 3478                & r1(-1, 1)*r_prev(-2, l-1) - r1(-1, -1)*r_prev(-2, 1-l)
 
 3479            r(1, l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3480            r(1, l) = r(1, l) / scal_uvw_m(l)
 
 3483            u = r1(0, 1)*r_prev(-1, l-1) - r1(0, -1)*r_prev(-1, 1-l)
 
 3484            v = r1(-1, 1)*r_prev(0, l-1) - r1(-1, -1)*r_prev(0, 1-l)
 
 3485            w = r1(1, 1)*r_prev(-2, l-1) - r1(1, -1)*r_prev(-2, 1-l) - &
 
 3486                & r1(-1, 1)*r_prev(2, l-1) + r1(-1, -1)*r_prev(2, 1-l)
 
 3487            r(-1, l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3488            r(-1, l) = r(-1, l) / scal_uvw_m(l)
 
 3491            u = r1(0, 1)*r_prev(0, l-1) - r1(0, -1)*r_prev(0, 1-l)
 
 3492            v = r1(1, 1)*r_prev(1, l-1) - r1(1, -1)*r_prev(1, 1-l) + &
 
 3493                & r1(-1, 1)*r_prev(-1, l-1) - r1(-1, -1)*r_prev(-1, 1-l)
 
 3494            r(0, l) = (u*scal_u_n(0) + v*scal_v_n(0)) / scal_uvw_m(l)
 
 3498                u = r1(0, 1)*r_prev(n, l-1) - r1(0, -1)*r_prev(n, 1-l)
 
 3499                v = r1(1, 1)*r_prev(n-1, l-1) - r1(1, -1)*r_prev(n-1, 1-l) - &
 
 3500                    & r1(-1, 1)*r_prev(1-n, l-1) + r1(-1, -1)*r_prev(1-n, 1-l)
 
 3501                w = r1(1, 1)*r_prev(n+1, l-1) - r1(1, -1)*r_prev(n+1, 1-l) + &
 
 3502                    & r1(-1, 1)*r_prev(-n-1, l-1) - r1(-1, -1)*r_prev(-n-1, 1-l)
 
 3503                r(n, l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3504                r(n, l) = r(n, l) / scal_uvw_m(l)
 
 3506                u = r1(0, 1)*r_prev(-n, l-1) - r1(0, -1)*r_prev(-n, 1-l)
 
 3507                v = r1(1, 1)*r_prev(1-n, l-1) - r1(1, -1)*r_prev(1-n, 1-l) + &
 
 3508                    & r1(-1, 1)*r_prev(n-1, l-1) - r1(-1, -1)*r_prev(n-1, 1-l)
 
 3509                w = r1(1, 1)*r_prev(-n-1, l-1) - r1(1, -1)*r_prev(-n-1, 1-l) - &
 
 3510                    & r1(-1, 1)*r_prev(n+1, l-1) + r1(-1, -1)*r_prev(n+1, 1-l)
 
 3511                r(-n, l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3512                r(-n, l) = r(-n, l) / scal_uvw_m(l)
 
 3515            dst(ind+l) = beta*dst(ind+l) + &
 
 3516                & alpha*dot_product(src(ind-l:ind+l), r(-l:l, l))
 
 3518            v = r1(1, 1)*r_prev(l-1, 1-l) + r1(1, -1)*r_prev(l-1, l-1) - &
 
 3519                & r1(-1, 1)*r_prev(1-l, 1-l) - r1(-1, -1)*r_prev(1-l, l-1)
 
 3520            r(l, -l) = v * scal_v_n(l) / scal_uvw_m(l)
 
 3523            v = r1(1, 1)*r_prev(1-l, 1-l) + r1(1, -1)*r_prev(1-l, l-1) + &
 
 3524                & r1(-1, 1)*r_prev(l-1, 1-l) + r1(-1, -1)*r_prev(l-1, l-1)
 
 3525            r(-l, -l) = v * scal_v_n(l) / scal_uvw_m(l)
 
 3528            u = r1(0, 1)*r_prev(l-1, 1-l) + r1(0, -1)*r_prev(l-1, l-1)
 
 3529            v = r1(1, 1)*r_prev(l-2, 1-l) + r1(1, -1)*r_prev(l-2, l-1) - &
 
 3530                & r1(-1, 1)*r_prev(2-l, 1-l) - r1(-1, -1)*r_prev(2-l, l-1)
 
 3531            r(l-1, -l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
 
 3534            u = r1(0, 1)*r_prev(1-l, 1-l) + r1(0, -1)*r_prev(1-l, l-1)
 
 3535            v = r1(1, 1)*r_prev(2-l, 1-l) + r1(1, -1)*r_prev(2-l, l-1) + &
 
 3536                & r1(-1, 1)*r_prev(l-2, 1-l) + r1(-1, -1)*r_prev(l-2, l-1)
 
 3537            r(1-l, -l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
 
 3540            u = r1(0, 1)*r_prev(1, 1-l) + r1(0, -1)*r_prev(1, l-1)
 
 3541            v = r1(1, 1)*r_prev(0, 1-l) + r1(1, -1)*r_prev(0, l-1)
 
 3542            w = r1(1, 1)*r_prev(2, 1-l) + r1(1, -1)*r_prev(2, l-1) + &
 
 3543                & r1(-1, 1)*r_prev(-2, 1-l) + r1(-1, -1)*r_prev(-2, l-1)
 
 3544            r(1, -l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3545            r(1, -l) = r(1, -l) / scal_uvw_m(l)
 
 3548            u = r1(0, 1)*r_prev(-1, 1-l) + r1(0, -1)*r_prev(-1, l-1)
 
 3549            v = r1(-1, 1)*r_prev(0, 1-l) + r1(-1, -1)*r_prev(0, l-1)
 
 3550            w = r1(1, 1)*r_prev(-2, 1-l) + r1(1, -1)*r_prev(-2, l-1) - &
 
 3551                & r1(-1, 1)*r_prev(2, 1-l) - r1(-1, -1)*r_prev(2, l-1)
 
 3552            r(-1, -l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3553            r(-1, -l) = r(-1, -l) / scal_uvw_m(l)
 
 3556            u = r1(0, 1)*r_prev(0, 1-l) + r1(0, -1)*r_prev(0, l-1)
 
 3557            v = r1(1, 1)*r_prev(1, 1-l) + r1(1, -1)*r_prev(1, l-1) + &
 
 3558                & r1(-1, 1)*r_prev(-1, 1-l) + r1(-1, -1)*r_prev(-1, l-1)
 
 3559            r(0, -l) = (u*scal_u_n(0) + v*scal_v_n(0)) / scal_uvw_m(l)
 
 3563                u = r1(0, 1)*r_prev(n, 1-l) + r1(0, -1)*r_prev(n, l-1)
 
 3564                v = r1(1, 1)*r_prev(n-1, 1-l) + r1(1, -1)*r_prev(n-1, l-1) - &
 
 3565                    & r1(-1, 1)*r_prev(1-n, 1-l) - r1(-1, -1)*r_prev(1-n, l-1)
 
 3566                w = r1(1, 1)*r_prev(n+1, 1-l) + r1(1, -1)*r_prev(n+1, l-1) + &
 
 3567                    & r1(-1, 1)*r_prev(-n-1, 1-l) + r1(-1, -1)*r_prev(-n-1, l-1)
 
 3568                r(n, -l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3569                r(n, -l) = r(n, -l) / scal_uvw_m(l)
 
 3571                u = r1(0, 1)*r_prev(-n, 1-l) + r1(0, -1)*r_prev(-n, l-1)
 
 3572                v = r1(1, 1)*r_prev(1-n, 1-l) + r1(1, -1)*r_prev(1-n, l-1) + &
 
 3573                    & r1(-1, 1)*r_prev(n-1, 1-l) + r1(-1, -1)*r_prev(n-1, l-1)
 
 3574                w = r1(1, 1)*r_prev(-n-1, 1-l) + r1(1, -1)*r_prev(-n-1, l-1) - &
 
 3575                    & r1(-1, 1)*r_prev(n+1, 1-l) - r1(-1, -1)*r_prev(n+1, l-1)
 
 3576                r(-n, -l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3577                r(-n, -l) = r(-n, -l) / scal_uvw_m(l)
 
 3580            dst(ind-l) = beta*dst(ind-l) + &
 
 3581                & alpha*dot_product(src(ind-l:ind+l), r(-l:l, -l))
 
 3585                v = r1(1, 0)*r_prev(l-1, m) - r1(-1, 0)*r_prev(1-l, m)
 
 3586                r(l, m) = v * scal_v_n(l) / scal_uvw_m(m)
 
 3589                v = r1(1, 0)*r_prev(1-l, m) + r1(-1, 0)*r_prev(l-1, m)
 
 3590                r(-l, m) = v * scal_v_n(l) / scal_uvw_m(m)
 
 3593                u = r1(0, 0) * r_prev(l-1, m)
 
 3594                v = r1(1, 0)*r_prev(l-2, m) - r1(-1, 0)*r_prev(2-l, m)
 
 3595                r(l-1, m) = u*scal_u_n(l-1) + v*scal_v_n(l-1)
 
 3596                r(l-1, m) = r(l-1, m) / scal_uvw_m(m)
 
 3599                u = r1(0, 0) * r_prev(1-l, m)
 
 3600                v = r1(1, 0)*r_prev(2-l, m) + r1(-1, 0)*r_prev(l-2, m)
 
 3601                r(1-l, m) = u*scal_u_n(l-1) + v*scal_v_n(l-1)
 
 3602                r(1-l, m) = r(1-l, m) / scal_uvw_m(m)
 
 3605                u = r1(0, 0) * r_prev(0, m)
 
 3606                v = r1(1, 0)*r_prev(1, m) + r1(-1, 0)*r_prev(-1, m)
 
 3607                r(0, m) = u*scal_u_n(0) + v*scal_v_n(0)
 
 3608                r(0, m) = r(0, m) / scal_uvw_m(m)
 
 3611                u = r1(0, 0) * r_prev(1, m)
 
 3612                v = r1(1, 0) * r_prev(0, m)
 
 3613                w = r1(1, 0)*r_prev(2, m) + r1(-1, 0)*r_prev(-2, m)
 
 3614                r(1, m) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3615                r(1, m) = r(1, m) / scal_uvw_m(m)
 
 3618                u = r1(0, 0) * r_prev(-1, m)
 
 3619                v = r1(-1, 0) * r_prev(0, m)
 
 3620                w = r1(1, 0)*r_prev(-2, m) - r1(-1, 0)*r_prev(2, m)
 
 3621                r(-1, m) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
 
 3622                r(-1, m) = r(-1, m) / scal_uvw_m(m)
 
 3626                    u = r1(0, 0) * r_prev(n, m)
 
 3627                    v = r1(1, 0)*r_prev(n-1, m) - r1(-1, 0)*r_prev(1-n, m)
 
 3628                    w = r1(1, 0)*r_prev(n+1, m) + r1(-1, 0)*r_prev(-1-n, m)
 
 3629                    r(n, m) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3630                    r(n, m) = r(n, m) / scal_uvw_m(m)
 
 3632                    u = r1(0, 0) * r_prev(-n, m)
 
 3633                    v = r1(1, 0)*r_prev(1-n, m) + r1(-1, 0)*r_prev(n-1, m)
 
 3634                    w = r1(1, 0)*r_prev(-n-1, m) - r1(-1, 0)*r_prev(n+1, m)
 
 3635                    r(-n, m) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
 
 3636                    r(-n, m) = r(-n, m) / scal_uvw_m(m)
 
 3639                dst(ind+m) = beta*dst(ind+m) + &
 
 3640                    & alpha*dot_product(src(ind-l:ind+l), r(-l:l, m))
 
 3669    integer, 
intent(in) :: p
 
 3670    real(dp), 
intent(in) :: vcos(p+1), vsin(p+1), alpha, src((p+1)*(p+1)), beta
 
 3672    real(dp), 
intent(inout) :: dst((p+1)*(p+1))
 
 3700    integer, 
intent(in) :: p
 
 3701    real(dp), 
intent(in) :: vcos(p+1), vsin(p+1), alpha, src((p+1)*(p+1)), beta
 
 3703    real(dp), 
intent(inout) :: dst((p+1)*(p+1))
 
 3705    integer :: l, m, ind
 
 3706    real(dp) :: v1, v2, v3, v4
 
 3708    if (alpha .eq. zero) 
then 
 3710        if (beta .eq. zero) 
then 
 3720    if (beta .eq. zero) 
then 
 3722        dst(1) = alpha*src(1)
 
 3728            dst(ind) = alpha*src(ind)
 
 3737                dst(ind+m) = alpha * (v1*v3-v2*v4)
 
 3739                dst(ind-m) = alpha * (v1*v4+v2*v3)
 
 3744        dst(1) = beta*dst(1) + alpha*src(1)
 
 3750            dst(ind) = beta*dst(ind) + alpha*src(ind)
 
 3759                dst(ind+m) = beta*dst(ind+m) + alpha*(v1*v3-v2*v4)
 
 3761                dst(ind-m) = beta*dst(ind-m) + alpha*(v1*v4+v2*v3)
 
 3790    integer, 
intent(in) :: p
 
 3791    real(dp), 
intent(in) :: vcos(p+1), vsin(p+1), alpha, src((p+1)*(p+1)), beta
 
 3793    real(dp), 
intent(inout) :: dst((p+1)*(p+1))
 
 3821    integer, 
intent(in) :: p
 
 3822    real(dp), 
intent(in) :: vcos(p+1), vsin(p+1), alpha, src((p+1)*(p+1)), beta
 
 3824    real(dp), 
intent(inout) :: dst((p+1)*(p+1))
 
 3826    integer :: l, m, ind
 
 3827    real(dp) :: v1, v2, v3, v4
 
 3829    if (alpha .eq. zero) 
then 
 3831        if (beta .eq. zero) 
then 
 3841    if (beta .eq. zero) 
then 
 3843        dst(1) = alpha*src(1)
 
 3848            dst(ind) = alpha*src(ind)
 
 3856                dst(ind+m) = alpha * (v1*v3+v2*v4)
 
 3858                dst(ind-m) = alpha * (v2*v3-v1*v4)
 
 3863        dst(1) = beta*dst(1) + alpha*src(1)
 
 3868            dst(ind) = beta*dst(ind) + alpha*src(ind)
 
 3876                dst(ind+m) = beta*dst(ind+m) + alpha*(v1*v3+v2*v4)
 
 3878                dst(ind-m) = beta*dst(ind-m) + alpha*(v2*v3-v1*v4)
 
 3918    integer, 
intent(in) :: p
 
 3919    real(dp), 
intent(in) :: ctheta, stheta, alpha, src((p+1)**2), beta
 
 3921    real(dp), 
intent(inout) :: dst((p+1)**2)
 
 3923    real(dp) :: work(2*(2*p+1)*(2*p+3)+p)
 
 3944    integer, 
intent(in) :: p
 
 3945    real(dp), 
intent(in) :: ctheta, stheta, alpha, src((p+1)**2), beta
 
 3947    real(dp), 
intent(out) :: dst((p+1)*(p+1))
 
 3949    real(dp), 
intent(out), 
target :: work(4*p*p+13*p+4)
 
 3951    real(dp) :: u, v, w, fl, fl2, tmp1, tmp2, vu(2), vv(2), vw(2), &
 
 3952        ctheta2, stheta2, cstheta
 
 3953    integer :: l, m, n, ind
 
 3955    real(dp), 
pointer :: r(:, :, :), r_prev(:, :, :), scal_uvw_m(:), &
 
 3956        & scal_u_n(:), scal_v_n(:), scal_w_n(:), r_swap(:, :, :), vsqr(:)
 
 3963    if (alpha .eq. zero) 
then 
 3965        if (beta .eq. zero) 
then 
 3975    if (beta .eq. zero) 
then 
 3978        dst(1) = alpha * src(1)
 
 3983        dst(2) = alpha * src(2)
 
 3984        dst(3) = alpha * (src(3)*ctheta - src(4)*stheta)
 
 3985        dst(4) = alpha * (src(3)*stheta + src(4)*ctheta)
 
 3990        l = 2 * (p+1) * (p+1)
 
 3991        r(1:2, 0:p, 0:p) => work(1:l)
 
 3993        r_prev(1:2, 0:p, 0:p) => work(l+1:m)
 
 3995        scal_uvw_m(0:p) => work(m+1:l)
 
 3997        scal_u_n(0:p-1) => work(l+1:m)
 
 3999        scal_v_n(0:p) => work(m+1:l)
 
 4001        scal_w_n(1:p-2) => work(l+1:m)
 
 4003        vsqr(1:p) => work(m+1:l)
 
 4005        ctheta2 = ctheta * ctheta
 
 4006        cstheta = ctheta * stheta
 
 4007        stheta2 = stheta * stheta
 
 4008        r(1, 2, 2) = (ctheta2 + one) / two
 
 4009        r(1, 1, 2) = cstheta
 
 4010        r(1, 0, 2) = sqrt3 / two * stheta2
 
 4011        dst(9) = alpha * (src(9)*r(1, 2, 2) + src(8)*r(1, 1, 2) + &
 
 4012            & src(7)*r(1, 0, 2))
 
 4013        r(1, 2, 1) = -cstheta
 
 4014        r(1, 1, 1) = ctheta2 - stheta2
 
 4015        r(1, 0, 1) = sqrt3 * cstheta
 
 4016        dst(8) = alpha * (src(9)*r(1, 2, 1) + src(8)*r(1, 1, 1) + &
 
 4017            & src(7)*r(1, 0, 1))
 
 4018        r(1, 2, 0) = sqrt3 / two * stheta2
 
 4019        r(1, 1, 0) = -sqrt3 * cstheta
 
 4020        r(1, 0, 0) = (three*ctheta2-one) / two
 
 4021        dst(7) = alpha * (src(9)*r(1, 2, 0) + src(8)*r(1, 1, 0) + &
 
 4022            & src(7)*r(1, 0, 0))
 
 4025        r(2, 2, 1) = -stheta
 
 4026        dst(6) = alpha * (src(6)*r(2, 1, 1) + src(5)*r(2, 2, 1))
 
 4029        dst(5) = alpha * (src(6)*r(2, 1, 2) + src(5)*r(2, 2, 2))
 
 4042            scal_uvw_m(0) = one / fl
 
 4045                u = sqrt(fl2 - vsqr(m))
 
 4046                scal_uvw_m(m) = one / u
 
 4049            u = sqrt(dble(u*(u-one)))
 
 4050            scal_uvw_m(l) = one / u
 
 4051            scal_u_n(0) = dble(l)
 
 4052            scal_v_n(0) = -sqrt(dble(l*(l-1))) / sqrt2
 
 4055                u = sqrt(fl2-vsqr(n))
 
 4061                v = sqrt(v*v-v) / two
 
 4064                w = -sqrt(w*w-w) / two
 
 4067            u = sqrt(dble(2*l-1))
 
 4069            v = sqrt(dble((2*l-1)*(2*l-2))) / two
 
 4071            v = sqrt(dble(2*l*(2*l-1))) / two
 
 4075            vv = ctheta*r_prev(:, l-1, l-1) + r_prev(2:1:-1, l-1, l-1)
 
 4076            r(:, l, l) = vv * scal_v_n(l) * scal_uvw_m(l)
 
 4077            tmp1 = src(ind+l) * r(1, l, l)
 
 4078            tmp2 = src(ind-l) * r(2, l, l)
 
 4080            vu = stheta * r_prev(:, l-1, l-1)
 
 4081            vv = ctheta*r_prev(:, l-2, l-1) + r_prev(2:1:-1, l-2, l-1)
 
 4082            r(:, l-1, l) = (vu*scal_u_n(l-1)+vv*scal_v_n(l-1)) * scal_uvw_m(l)
 
 4083            tmp1 = tmp1 + src(ind+l-1)*r(1, l-1, l)
 
 4084            tmp2 = tmp2 + src(ind-l+1)*r(2, l-1, l)
 
 4086            vu = stheta * r_prev(:, 1, l-1)
 
 4087            vv(1) = ctheta * r_prev(1, 0, l-1)
 
 4088            vv(2) = r_prev(1, 0, l-1)
 
 4089            vw = ctheta*r_prev(:, 2, l-1) - r_prev(2:1:-1, 2, l-1)
 
 4090            r(:, 1, l) = vu*scal_u_n(1) + vw*scal_w_n(1) + sqrt2*scal_v_n(1)*vv
 
 4091            r(:, 1, l) = r(:, 1, l) * scal_uvw_m(l)
 
 4092            tmp1 = tmp1 + src(ind+1)*r(1, 1, l)
 
 4093            tmp2 = tmp2 + src(ind-1)*r(2, 1, l)
 
 4095            u = stheta * r_prev(1, 0, l-1)
 
 4096            v = ctheta*r_prev(1, 1, l-1) - r_prev(2, 1, l-1)
 
 4097            r(1, 0, l) = (u*scal_u_n(0) + v*scal_v_n(0)) * scal_uvw_m(l)
 
 4098            tmp1 = tmp1 + src(ind)*r(1, 0, l)
 
 4102                vu = stheta * r_prev(:, n, l-1)
 
 4103                vv = ctheta*r_prev(:, n-1, l-1) + r_prev(2:1:-1, n-1, l-1)
 
 4104                vw = ctheta*r_prev(:, n+1, l-1) - r_prev(2:1:-1, n+1, l-1)
 
 4105                vu = vu*scal_u_n(n) + vv*scal_v_n(n) + vw*scal_w_n(n)
 
 4106                r(:, n, l) = vu * scal_uvw_m(l)
 
 4107                tmp1 = tmp1 + src(ind+n)*r(1, n, l)
 
 4108                tmp2 = tmp2 + src(ind-n)*r(2, n, l)
 
 4110            dst(ind+l) = alpha * tmp1
 
 4111            dst(ind-l) = alpha * tmp2
 
 4114            v = -stheta * r_prev(1, l-1, 0)
 
 4115            u = scal_v_n(l) * scal_uvw_m(0)
 
 4117            tmp1 = src(ind+l) * r(1, l, 0)
 
 4119            u = ctheta * r_prev(1, l-1, 0)
 
 4120            v = -stheta * r_prev(1, l-2, 0)
 
 4121            w = u*scal_u_n(l-1) + v*scal_v_n(l-1)
 
 4122            r(1, l-1, 0) = w * scal_uvw_m(0)
 
 4123            tmp1 = tmp1 + src(ind+l-1)*r(1, l-1, 0)
 
 4125            u = ctheta * r_prev(1, 0, 0)
 
 4126            v = -stheta * r_prev(1, 1, 0)
 
 4127            w = u*scal_u_n(0) + v*scal_v_n(0)
 
 4128            r(1, 0, 0) = w * scal_uvw_m(0)
 
 4129            tmp1 = tmp1 + src(ind)*r(1, 0, 0)
 
 4131            v = sqrt2*scal_v_n(1)*r_prev(1, 0, 0) + &
 
 4132                & scal_w_n(1)*r_prev(1, 2, 0)
 
 4133            u = ctheta * r_prev(1, 1, 0)
 
 4134            w = scal_u_n(1)*u - stheta*v
 
 4135            r(1, 1, 0) = w * scal_uvw_m(0)
 
 4136            tmp1 = tmp1 + src(ind+1)*r(1, 1, 0)
 
 4140                v = scal_v_n(n)*r_prev(1, n-1, 0) + &
 
 4141                    & scal_w_n(n)*r_prev(1, n+1, 0)
 
 4142                u = ctheta * r_prev(1, n, 0)
 
 4143                w = scal_u_n(n)*u - stheta*v
 
 4144                r(1, n, 0) = w * scal_uvw_m(0)
 
 4145                tmp1 = tmp1 + src(ind+n)*r(1, n, 0)
 
 4147            dst(ind) = alpha * tmp1
 
 4152                vv = -stheta * r_prev(:, l-1, m)
 
 4153                u = scal_v_n(l) * scal_uvw_m(m)
 
 4155                tmp1 = src(ind+l) * r(1, l, m)
 
 4156                tmp2 = src(ind-l) * r(2, l, m)
 
 4158                vu = ctheta * r_prev(:, l-1, m)
 
 4159                vv = -stheta * r_prev(:, l-2, m)
 
 4160                vw = vu*scal_u_n(l-1) + vv*scal_v_n(l-1)
 
 4161                r(:, l-1, m) = vw * scal_uvw_m(m)
 
 4162                tmp1 = tmp1 + src(ind+l-1)*r(1, l-1, m)
 
 4163                tmp2 = tmp2 + src(ind-l+1)*r(2, l-1, m)
 
 4165                u = ctheta * r_prev(1, 0, m)
 
 4166                v = -stheta * r_prev(1, 1, m)
 
 4167                w = u*scal_u_n(0) + v*scal_v_n(0)
 
 4168                r(1, 0, m) = w * scal_uvw_m(m)
 
 4169                tmp1 = tmp1 + src(ind)*r(1, 0, m)
 
 4171                v = sqrt2*scal_v_n(1)*r_prev(1, 0, m) + &
 
 4172                    & scal_w_n(1)*r_prev(1, 2, m)
 
 4173                u = ctheta * r_prev(1, 1, m)
 
 4174                w = scal_u_n(1)*u - stheta*v
 
 4175                r(1, 1, m) = w * scal_uvw_m(m)
 
 4176                tmp1 = tmp1 + src(ind+1)*r(1, 1, m)
 
 4178                u = ctheta * r_prev(2, 1, m)
 
 4179                w = -stheta * r_prev(2, 2, m)
 
 4180                v = u*scal_u_n(1) + w*scal_w_n(1)
 
 4181                r(2, 1, m) = v * scal_uvw_m(m)
 
 4182                tmp2 = tmp2 + src(ind-1)*r(2, 1, m)
 
 4186                    vv = scal_v_n(n)*r_prev(:, n-1, m) + &
 
 4187                        & scal_w_n(n)*r_prev(:, n+1, m)
 
 4188                    vu = ctheta * r_prev(:, n, m)
 
 4189                    vw = scal_u_n(n)*vu - stheta*vv
 
 4190                    r(:, n, m) = vw * scal_uvw_m(m)
 
 4191                    tmp1 = tmp1 + src(ind+n)*r(1, n, m)
 
 4192                    tmp2 = tmp2 + src(ind-n)*r(2, n, m)
 
 4194                dst(ind+m) = alpha * tmp1
 
 4195                dst(ind-m) = alpha * tmp2
 
 4199        stop 
"Not Implemented" 
 4226        & src_m, beta, dst_m)
 
 4228    integer, 
intent(in) :: p
 
 4229    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
 
 4230        & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
 
 4232    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 4234    real(dp) :: work(2*(p+1))
 
 4237        & src_m, beta, dst_m, work)
 
 4264        & src_m, beta, dst_m, work)
 
 4266    integer, 
intent(in) :: p
 
 4267    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
 
 4268        & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
 
 4270    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 4272    real(dp), 
intent(out), 
target :: work(2*(p+1))
 
 4274    real(dp) :: r1, r2, tmp1, tmp2, tmp3, res1, res2, pow_r1
 
 4275    integer :: j, k, n, indj, indjn, indjk1, indjk2
 
 4277    real(dp), 
pointer :: pow_r2(:)
 
 4279    if (alpha .eq. zero) 
then 
 4280        if (beta .eq. zero) 
then 
 4283            dst_m = beta * dst_m
 
 4291        pow_r2(1:p+1) => work(1:p+1)
 
 4300            pow_r2(j) = pow_r2(j-1) * r2
 
 4304        if (beta .eq. zero) 
then 
 4309                tmp1 = alpha * pow_r1 * vscales(indj)
 
 4310                pow_r1 = pow_r1 * r1
 
 4314                indjk1 = j*(j+1)/2 + 1
 
 4317                    indjn = (j-n)**2 + (j-n) + 1
 
 4318                    tmp3 = pow_r2(n+1) / &
 
 4319                        & vscales(indjn) * vcnk(indjk1+n)**2
 
 4320                    res1 = res1 + tmp3*src_m(indjn)
 
 4322                dst_m(indj) = tmp2 * res1
 
 4329                    indjk1 = (j-k)*(j-k+1)/2 + 1
 
 4330                    indjk2 = (j+k)*(j+k+1)/2 + 1
 
 4333                        indjn = (j-n)**2 + (j-n) + 1
 
 4334                        tmp3 = pow_r2(n+1) / &
 
 4335                            & vscales(indjn) * vcnk(indjk1+n) * &
 
 4337                        res1 = res1 + tmp3*src_m(indjn+k)
 
 4338                        res2 = res2 + tmp3*src_m(indjn-k)
 
 4340                    dst_m(indj+k) = tmp2 * res1
 
 4341                    dst_m(indj-k) = tmp2 * res2
 
 4350                tmp1 = alpha * pow_r1 * vscales(indj)
 
 4351                pow_r1 = pow_r1 * r1
 
 4355                indjk1 = j * (j+1) /2 + 1
 
 4358                    indjn = (j-n)**2 + (j-n) + 1
 
 4359                    tmp3 = pow_r2(n+1) / &
 
 4360                        & vscales(indjn) * vcnk(indjk1+n)**2
 
 4361                    res1 = res1 + tmp3*src_m(indjn)
 
 4363                dst_m(indj) = beta*dst_m(indj) + tmp2*res1
 
 4370                    indjk1 = (j-k)*(j-k+1)/2 + 1
 
 4371                    indjk2 = (j+k)*(j+k+1)/2 + 1
 
 4374                        indjn = (j-n)**2 + (j-n) + 1
 
 4375                        tmp3 = pow_r2(n+1) / &
 
 4376                            & vscales(indjn) * vcnk(indjk1+n) * &
 
 4378                        res1 = res1 + tmp3*src_m(indjn+k)
 
 4379                        res2 = res2 + tmp3*src_m(indjn-k)
 
 4381                    dst_m(indj+k) = beta*dst_m(indj+k) + tmp2*res1
 
 4382                    dst_m(indj-k) = beta*dst_m(indj-k) + tmp2*res2
 
 4389        if (beta .eq. zero) 
then 
 4394                do k = indj-j, indj+j
 
 4395                    dst_m(k) = tmp1 * src_m(k)
 
 4405                do k = indj-j, indj+j
 
 4406                    dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
 
 4436        & src_m, beta, dst_m)
 
 4438    integer, 
intent(in) :: p
 
 4439    real(dp), 
intent(in) :: z, src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
 
 4440        & alpha, src_m((p+1)*(p+1)), beta
 
 4442    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 4444    real(dp) :: work(2*p+1)
 
 4445    complex(dp) :: work_complex(2*p+1)
 
 4448        & alpha, src_m, beta, dst_m, work, work_complex)
 
 4475        & alpha, src_m, beta, dst_m, work, work_complex)
 
 4478    integer, 
intent(in) :: p
 
 4479    real(dp), 
intent(in) :: z, src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
 
 4480        & alpha, src_m((p+1)*(p+1)), beta
 
 4482    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 4484    real(dp), 
intent(out), 
target :: work(2*p+1)
 
 4485    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 4487    real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
 
 4488    integer :: j, k, n, indj, indn, l, NZ, ierr
 
 4490    real(dp) :: fact(2*p+1)
 
 4491    complex(dp) :: complex_argument
 
 4493    if (alpha .eq. zero) 
then 
 4494        if (beta .eq. zero) 
then 
 4497            dst_m = beta * dst_m
 
 4504        fact(j+1) = dble(j) * fact(j)
 
 4509        complex_argument = z
 
 4510        call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 4511        work(1) = sqrt(pi/(2*z)) * real(work_complex(1))
 
 4513            call cbesi(complex_argument, 1.5d0, 1, 2*p, &
 
 4514                & work_complex(2:2*p+1), nz, ierr)
 
 4515            work(2:2*p+1) = sqrt(pi/(2*z)) * real(work_complex(2:2*p+1))
 
 4519        if (beta .eq. zero) 
then 
 4530                    tmp1 = vscales(indn) * &
 
 4531                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 4532                        & src_sk(n+1) / vscales(indj) * dst_sk(j+1)
 
 4535                        tmp2 = (two**(-l)) / &
 
 4536                            & fact(l+1)*fact(2*l+1)/ &
 
 4537                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 4538                            & (work(j+n-l+1) / (z**l))
 
 4542                    res1 = res1 + tmp3*src_m(indn)
 
 4552                        tmp1 = vscales(indn+k) * &
 
 4553                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 4554                            & src_sk(n+1) / vscales(indj+k) * dst_sk(j+1)
 
 4557                            tmp2 = (two**(-l)) / &
 
 4558                                & fact(l+k+1)*fact(2*l+1)/ &
 
 4559                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 4560                                & (work(j+n-l+1) / (z**l))
 
 4564                        res1 = res1 + tmp3*src_m(indn+k)
 
 4565                        res2 = res2 + tmp3*src_m(indn-k)
 
 4567                    dst_m(indj+k) = res1
 
 4568                    dst_m(indj-k) = res2
 
 4582                    tmp1 = vscales(indn) * &
 
 4583                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 4584                        & src_sk(n+1) / vscales(indj) * dst_sk(j+1)
 
 4587                        tmp2 = (two**(-l)) / &
 
 4588                            & fact(l+1)*fact(2*l+1)/ &
 
 4589                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 4590                            & (work(j+n-l+1) / (z**l))
 
 4594                    res1 = res1 + tmp3*src_m(indn)
 
 4596                dst_m(indj) = beta*dst_m(indj) + res1
 
 4604                        tmp1 = vscales(indn+k) * &
 
 4605                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 4606                            & src_sk(n+1) / vscales(indj+k) * dst_sk(j+1)
 
 4609                            tmp2 = (two**(-l)) / &
 
 4610                                & fact(l+k+1)*fact(2*l+1)/ &
 
 4611                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 4612                                & (work(j+n-l+1) / (z**l))
 
 4616                        res1 = res1 + tmp3*src_m(indn+k)
 
 4617                        res2 = res2 + tmp3*src_m(indn-k)
 
 4619                    dst_m(indj+k) = beta*dst_m(indj+k) + res1
 
 4620                    dst_m(indj-k) = beta*dst_m(indj-k) + res2
 
 4628        if (beta .eq. zero) 
then 
 4633                do k = indj-j, indj+j
 
 4634                    dst_m(k) = tmp1 * src_m(k)
 
 4644                do k = indj-j, indj+j
 
 4645                    dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
 
 4675        & src_m, beta, dst_m)
 
 4677    integer, 
intent(in) :: p
 
 4678    real(dp), 
intent(in) :: z, src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
 
 4679        & alpha, src_m((p+1)*(p+1)), beta
 
 4681    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 4683    real(dp) :: work(2*p+1)
 
 4684    complex(dp) :: work_complex(2*p+1)
 
 4687        & alpha, src_m, beta, dst_m, work, work_complex)
 
 4714        & alpha, src_m, beta, dst_m, work, work_complex)
 
 4717    integer, 
intent(in) :: p
 
 4718    real(dp), 
intent(in) :: z, src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
 
 4719        & alpha, src_m((p+1)*(p+1)), beta
 
 4721    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 4723    real(dp), 
intent(out), 
target :: work(2*p+1)
 
 4724    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 4726    real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
 
 4727    integer :: j, k, n, indj, indn, l, NZ, ierr
 
 4729    real(dp) :: fact(2*p+1)
 
 4730    complex(dp) :: complex_argument
 
 4732    if (alpha .eq. zero) 
then 
 4733        if (beta .eq. zero) 
then 
 4736            dst_m = beta * dst_m
 
 4743        fact(j+1) = dble(j) * fact(j)
 
 4748        complex_argument = z
 
 4749        call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 4750        work(1) = sqrt(pi/(2*z)) * real(work_complex(1))
 
 4752            call cbesi(complex_argument, 1.5d0, 1, 2*p, &
 
 4753                & work_complex(2:2*p+1), nz, ierr)
 
 4754            work(2:2*p+1) = sqrt(pi/(2*z)) * real(work_complex(2:2*p+1))
 
 4758        if (beta .eq. zero) 
then 
 4769                    tmp1 = vscales(indn) * &
 
 4770                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 4771                        & src_sk(n+1) / vscales(indj) * dst_sk(j+1)
 
 4774                        tmp2 = (two**(-l)) / &
 
 4775                            & fact(l+1)*fact(2*l+1)/ &
 
 4776                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 4777                            & (work(j+n-l+1) / (z**l))
 
 4781                    res1 = res1 + tmp3*src_m(indj)
 
 4791                        tmp1 = vscales(indn+k) * &
 
 4792                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 4793                            & src_sk(n+1) / vscales(indj+k) * dst_sk(j+1)
 
 4796                            tmp2 = (two**(-l)) / &
 
 4797                                & fact(l+k+1)*fact(2*l+1)/ &
 
 4798                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 4799                                & (work(j+n-l+1) / (z**l))
 
 4803                        res1 = res1 + tmp3*src_m(indj+k)
 
 4804                        res2 = res2 + tmp3*src_m(indj-k)
 
 4806                    dst_m(indn+k) = res1
 
 4807                    dst_m(indn-k) = res2
 
 4821                    tmp1 = vscales(indn) * &
 
 4822                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 4823                        & src_sk(n+1) / vscales(indj) * dst_sk(j+1)
 
 4826                        tmp2 = (two**(-l)) / &
 
 4827                            & fact(l+1)*fact(2*l+1)/ &
 
 4828                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 4829                            & (work(j+n-l+1) / (z**l))
 
 4833                    res1 = res1 + tmp3*src_m(indj)
 
 4835                dst_m(indn) = beta*dst_m(indn) + res1
 
 4843                        tmp1 = vscales(indn+k) * &
 
 4844                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 4845                            & src_sk(n+1) / vscales(indj+k) * dst_sk(j+1)
 
 4848                            tmp2 = (two**(-l)) / &
 
 4849                                & fact(l+k+1)*fact(2*l+1)/ &
 
 4850                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 4851                                & (work(j+n-l+1) / (z**l))
 
 4855                        res1 = res1 + tmp3*src_m(indj+k)
 
 4856                        res2 = res2 + tmp3*src_m(indj-k)
 
 4858                    dst_m(indn+k) = beta*dst_m(indn+k) + res1
 
 4859                    dst_m(indn-k) = beta*dst_m(indn-k) + res2
 
 4867        if (beta .eq. zero) 
then 
 4872                do k = indj-j, indj+j
 
 4873                    dst_m(k) = tmp1 * src_m(k)
 
 4883                do k = indj-j, indj+j
 
 4884                    dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
 
 4916        & alpha, src_m, beta, dst_m)
 
 4919    integer, 
intent(in) :: p
 
 4920    real(dp), 
intent(in) :: src_sk(p+2), vscales((p+2)*(p+2)), &
 
 4921        & alpha, src_m((p+1)*(p+1)), beta
 
 4923    real(dp), 
intent(inout) :: dst_m((p+2)*(p+2))
 
 4926    real(dp) :: tmp1, tmp2, tmp3, res1, res2
 
 4927    integer :: j, k, n, indj, indn, l
 
 4929    real(dp) :: fact(2*p+1), fact2(p+2)
 
 4931    if (alpha .eq. zero) 
then 
 4932        if (beta .eq. zero) 
then 
 4935            dst_m = beta * dst_m
 
 4942        fact(j+1) = dble(j) * fact(j)
 
 4946        fact2(j+1) = dble(2*j+1) * fact2(j)
 
 4951    if (beta .eq. zero) 
then 
 4960        dst_m(1) = vscales(3) / src_sk(2) / vscales(1) * src_sk(1) / fact2(2) * &
 
 4973            tmp1 = vscales(indn) * &
 
 4974                & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 4975                & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
 
 4979            tmp2 = (two**(-l)) / &
 
 4980                & fact(l+1)*fact(2*l+1)/ &
 
 4981                & fact(l+1)/fact(l+1)/ &
 
 4985            res1 = res1 + tmp3*src_m(indn)
 
 4990            tmp1 = vscales(indn) * &
 
 4991                & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 4992                & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
 
 4996            tmp2 = (two**(-l)) / &
 
 4997                & fact(l+1)*fact(2*l+1)/ &
 
 4998                & fact(l+1)/fact(l+1)/ &
 
 5002            res1 = res1 + tmp3*src_m(indn)
 
 5012                tmp1 = vscales(indn+k) * &
 
 5013                    & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 5014                    & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5017                tmp2 = (two**(-l)) / &
 
 5018                    & fact(l+k+1)*fact(2*l+1)/ &
 
 5019                    & fact(l+1)/fact(l-k+1)/ &
 
 5023                res1 = res1 + tmp3*src_m(indn+k)
 
 5024                res2 = res2 + tmp3*src_m(indn-k)
 
 5029                tmp1 = vscales(indn+k) * &
 
 5030                    & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 5031                    & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5034                tmp2 = (two**(-l)) / &
 
 5035                    & fact(l+k+1)*fact(2*l+1)/ &
 
 5036                    & fact(l+1)/fact(l-k+1)/ &
 
 5040                res1 = res1 + tmp3*src_m(indn+k)
 
 5041                res2 = res2 + tmp3*src_m(indn-k)
 
 5042                dst_m(indj+k) = res1
 
 5043                dst_m(indj-k) = res2
 
 5053            tmp1 = vscales(indn+k) * &
 
 5054                & dble(2*j+1) * fact(n+k+1) / &
 
 5055                & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5058            tmp2 = (two**(-l)) / &
 
 5059                & fact(l+k+1)*fact(2*l+1)/ &
 
 5064            res1 = res1 + tmp3*src_m(indn+k)
 
 5065            res2 = res2 + tmp3*src_m(indn-k)
 
 5066            dst_m(indj+k) = res1
 
 5067            dst_m(indj-k) = res2
 
 5079        tmp1 = vscales(indn) * &
 
 5080            & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 5081            & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
 
 5084        tmp2 = (two**(-l)) / &
 
 5085            & fact(l+1)*fact(2*l+1)/ &
 
 5086            & fact(l+1)/fact(l+1)/ &
 
 5090        res1 = res1 + tmp3*src_m(indn)
 
 5100            tmp1 = vscales(indn+k) * &
 
 5101                & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 5102                & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5105            tmp2 = (two**(-l)) / &
 
 5106                & fact(l+k+1)*fact(2*l+1)/ &
 
 5107                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)/ &
 
 5111            res1 = res1 + tmp3*src_m(indn+k)
 
 5112            res2 = res2 + tmp3*src_m(indn-k)
 
 5113            dst_m(indj+k) = res1
 
 5114            dst_m(indj-k) = res2
 
 5127        tmp1 = vscales(indn) * &
 
 5128            & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 5129            & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
 
 5132        tmp2 = (two**(-l)) / &
 
 5133            & fact(l+1)*fact(2*l+1)/ &
 
 5134            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)/ &
 
 5138        res1 = res1 + tmp3*src_m(indn)
 
 5148            tmp1 = vscales(indn+k) * &
 
 5149                & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 5150                & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5153            tmp2 = (two**(-l)) / &
 
 5154                & fact(l+k+1)*fact(2*l+1)/ &
 
 5155                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)/ &
 
 5159            res1 = res1 + tmp3*src_m(indn+k)
 
 5160            res2 = res2 + tmp3*src_m(indn-k)
 
 5161            dst_m(indj+k) = res1
 
 5162            dst_m(indj-k) = res2
 
 5173        dst_m(1) = vscales(3) / src_sk(2) / vscales(1) * src_sk(1) / fact2(2) * &
 
 5174            & src_m(3) * alpha + dst_m(1)*beta
 
 5186            tmp1 = vscales(indn) * &
 
 5187                & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 5188                & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
 
 5192            tmp2 = (two**(-l)) / &
 
 5193                & fact(l+1)*fact(2*l+1)/ &
 
 5194                & fact(l+1)/fact(l+1)/ &
 
 5198            res1 = res1 + tmp3*src_m(indn)
 
 5203            tmp1 = vscales(indn) * &
 
 5204                & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 5205                & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
 
 5209            tmp2 = (two**(-l)) / &
 
 5210                & fact(l+1)*fact(2*l+1)/ &
 
 5211                & fact(l+1)/fact(l+1)/ &
 
 5215            res1 = res1 + tmp3*src_m(indn)
 
 5216            dst_m(indj) = beta*dst_m(indj) + res1
 
 5225                tmp1 = vscales(indn+k) * &
 
 5226                    & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 5227                    & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5230                tmp2 = (two**(-l)) / &
 
 5231                    & fact(l+k+1)*fact(2*l+1)/ &
 
 5232                    & fact(l+1)/fact(l-k+1)/ &
 
 5236                res1 = res1 + tmp3*src_m(indn+k)
 
 5237                res2 = res2 + tmp3*src_m(indn-k)
 
 5242                tmp1 = vscales(indn+k) * &
 
 5243                    & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 5244                    & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5247                tmp2 = (two**(-l)) / &
 
 5248                    & fact(l+k+1)*fact(2*l+1)/ &
 
 5249                    & fact(l+1)/fact(l-k+1)/ &
 
 5253                res1 = res1 + tmp3*src_m(indn+k)
 
 5254                res2 = res2 + tmp3*src_m(indn-k)
 
 5255                dst_m(indj+k) = beta*dst_m(indj+k) + res1
 
 5256                dst_m(indj-k) = beta*dst_m(indj-k) + res2
 
 5266            tmp1 = vscales(indn+k) * &
 
 5267                & dble(2*j+1) * fact(n+k+1) / &
 
 5268                & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5271            tmp2 = (two**(-l)) / &
 
 5272                & fact(l+k+1)*fact(2*l+1)/ &
 
 5277            res1 = res1 + tmp3*src_m(indn+k)
 
 5278            res2 = res2 + tmp3*src_m(indn-k)
 
 5279            dst_m(indj+k) = beta*dst_m(indj+k) + res1
 
 5280            dst_m(indj-k) = beta*dst_m(indj-k) + res2
 
 5292        tmp1 = vscales(indn) * &
 
 5293            & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 5294            & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
 
 5297        tmp2 = (two**(-l)) / &
 
 5298            & fact(l+1)*fact(2*l+1)/ &
 
 5299            & fact(l+1)/fact(l+1)/ &
 
 5303        res1 = res1 + tmp3*src_m(indn)
 
 5304        dst_m(indj) = beta*dst_m(indj) + res1
 
 5313            tmp1 = vscales(indn+k) * &
 
 5314                & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 5315                & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5318            tmp2 = (two**(-l)) / &
 
 5319                & fact(l+k+1)*fact(2*l+1)/ &
 
 5320                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)/ &
 
 5324            res1 = res1 + tmp3*src_m(indn+k)
 
 5325            res2 = res2 + tmp3*src_m(indn-k)
 
 5326            dst_m(indj+k) = beta*dst_m(indj+k) + res1
 
 5327            dst_m(indj-k) = beta*dst_m(indj-k) + res2
 
 5340        tmp1 = vscales(indn) * &
 
 5341            & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 5342            & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
 
 5345        tmp2 = (two**(-l)) / &
 
 5346            & fact(l+1)*fact(2*l+1)/ &
 
 5347            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)/ &
 
 5351        res1 = res1 + tmp3*src_m(indn)
 
 5352        dst_m(indj) = beta*dst_m(indj) + res1
 
 5361            tmp1 = vscales(indn+k) * &
 
 5362                & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 5363                & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
 
 5366            tmp2 = (two**(-l)) / &
 
 5367                & fact(l+k+1)*fact(2*l+1)/ &
 
 5368                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)/ &
 
 5372            res1 = res1 + tmp3*src_m(indn+k)
 
 5373            res2 = res2 + tmp3*src_m(indn-k)
 
 5374            dst_m(indj+k) = beta*dst_m(indj+k) + res1
 
 5375            dst_m(indj-k) = beta*dst_m(indj-k) + res2
 
 5403        & alpha, src_m, beta, dst_m)
 
 5405    integer, 
intent(in) :: p
 
 5406    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
 
 5407        & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
 
 5409    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5411    real(dp) :: work(2*(p+1))
 
 5414        & alpha, src_m, beta, dst_m, work)
 
 5441        & alpha, src_m, beta, dst_m, work)
 
 5443    integer, 
intent(in) :: p
 
 5444    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
 
 5445        & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
 
 5447    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5449    real(dp), 
intent(out), 
target :: work(2*(p+1))
 
 5451    real(dp) :: r1, r2, tmp1, tmp2, res1, res2
 
 5452    integer :: j, k, n, indj, indn, indjk1, indjk2
 
 5454    real(dp), 
pointer :: pow_r1(:), pow_r2(:)
 
 5456    if (alpha .eq. zero) 
then 
 5457        if (beta .eq. zero) 
then 
 5460            dst_m = beta * dst_m
 
 5467        pow_r1(1:p+1) => work(1:p+1)
 
 5468        pow_r2(1:p+1) => work(p+2:2*(p+1))
 
 5475            pow_r1(j) = pow_r1(j-1) * r1
 
 5476            pow_r2(j) = pow_r2(j-1) * r2
 
 5480        if (beta .eq. zero) 
then 
 5484                tmp1 = alpha * pow_r1(n+1) / vscales(indn)
 
 5491                    indjk1 = j*(j+1)/2 + 1
 
 5492                    tmp2 = tmp1 * vscales(indj) * pow_r2(j-n+1) * &
 
 5493                        & vcnk(indjk1+j-n)**2
 
 5494                    res1 = res1 + tmp2*src_m(indj)
 
 5505                        indjk1 = (j-k)*(j-k+1)/2 + 1
 
 5506                        indjk2 = (j+k)*(j+k+1)/2 + 1
 
 5507                        tmp2 = tmp1 * vscales(indj) * pow_r2(j-n+1) * &
 
 5508                            & vcnk(indjk1+j-n) * vcnk(indjk2+j-n)
 
 5509                        res1 = res1 + tmp2*src_m(indj+k)
 
 5510                        res2 = res2 + tmp2*src_m(indj-k)
 
 5512                    dst_m(indn+k) = res1
 
 5513                    dst_m(indn-k) = res2
 
 5521                tmp1 = alpha * pow_r1(n+1) / vscales(indn)
 
 5528                    indjk1 = j*(j+1)/2 + 1
 
 5529                    tmp2 = tmp1 * vscales(indj) * pow_r2(j-n+1) * &
 
 5530                        & vcnk(indjk1+j-n)**2
 
 5531                    res1 = res1 + tmp2*src_m(indj)
 
 5533                dst_m(indn) = beta*dst_m(indn) + res1
 
 5542                        indjk1 = (j-k)*(j-k+1)/2 + 1
 
 5543                        indjk2 = (j+k)*(j+k+1)/2 + 1
 
 5544                        tmp2 = tmp1 * vscales(indj) * pow_r2(j-n+1) * &
 
 5545                            & vcnk(indjk1+j-n) * vcnk(indjk2+j-n)
 
 5546                        res1 = res1 + tmp2*src_m(indj+k)
 
 5547                        res2 = res2 + tmp2*src_m(indj-k)
 
 5549                    dst_m(indn+k) = beta*dst_m(indn+k) + res1
 
 5550                    dst_m(indn-k) = beta*dst_m(indn-k) + res2
 
 5557        if (beta .eq. zero) 
then 
 5563                do k = indj-j, indj+j
 
 5564                    dst_m(k) = tmp1 * src_m(k)
 
 5575                do k = indj-j, indj+j
 
 5576                    dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
 
 5605    integer, 
intent(in) :: p
 
 5606    real(dp), 
intent(in) :: src_r, dst_r, alpha, src_m((p+1)*(p+1)), beta
 
 5608    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5610    real(dp) :: r1, tmp1
 
 5611    integer :: j, k, indj
 
 5613    if (beta .eq. zero) 
then 
 5616        dst_m = beta * dst_m
 
 5622        do k = indj-j, indj+j
 
 5623            dst_m(k) = dst_m(k) + src_m(k)*tmp1
 
 5651    integer, 
intent(in) :: p
 
 5652    real(dp), 
intent(in) :: src_r, dst_r, alpha, src_m((p+1)*(p+1)), beta
 
 5654    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5656    real(dp) :: r1, tmp1
 
 5657    integer :: j, k, indj
 
 5659    if (beta .eq. zero) 
then 
 5662        dst_m = beta * dst_m
 
 5668        do k = indj-j, indj+j
 
 5669            dst_m(k) = dst_m(k) + src_m(k)*tmp1
 
 5701        & src_m, beta, dst_m)
 
 5703    integer, 
intent(in) :: p
 
 5704    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 5705        & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
 
 5707    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5709    real(dp) :: work(6*p*p + 19*p + 8)
 
 5712        & src_m, beta, dst_m, work)
 
 5742        & src_m, beta, dst_m, work)
 
 5744    integer, 
intent(in) :: p
 
 5745    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 5746        & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
 
 5748    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5750    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 5752    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 5755    real(dp), 
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
 
 5757    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 5759    if (stheta .eq. zero) 
then 
 5762            & alpha, src_m, beta, dst_m, work)
 
 5768    tmp_m(1:m) => work(n+1:n+m) 
 
 5770    tmp_m2(1:m) => work(n+1:n+m) 
 
 5773    vcos => work(n+1:n+m) 
 
 5775    vsin => work(n+1:n+m) 
 
 5777    call trgev(cphi, sphi, p, vcos, vsin)
 
 5785        & tmp_m2, zero, tmp_m, work)
 
 5818        & src_m, beta, dst_m)
 
 5821    integer, 
intent(in) :: p
 
 5822    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 5823        & alpha, src_m((p+1)*(p+1)), beta, kappa
 
 5825    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5827    real(dp) :: work(6*p*p + 19*p + 8), src_sk(p+1), dst_sk(p+1), s1, s2, ck(3)
 
 5828    complex(dp) :: work_complex(2*p+1), z1, z2
 
 5833    s1 = sqrt(2/(pi*real(z1)))
 
 5834    s2 = sqrt(2/(pi*real(z2)))
 
 5835    call cbesk(z1, pt5, 1, 1, work_complex(1), nz, ierr)
 
 5836    src_sk(1) = s1 * real(work_complex(1))
 
 5837    call cbesk(z2, pt5, 1, 1, work_complex(1), nz, ierr)
 
 5838    dst_sk(1) = s2 * real(work_complex(1))
 
 5840        call cbesk(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 5841        src_sk(2:p+1) = s1 * real(work_complex(2:p+1))
 
 5842        call cbesk(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 5843        dst_sk(2:p+1) = s2 * real(work_complex(2:p+1))
 
 5848        & src_m, beta, dst_m, work, work_complex)
 
 5877        & src_m, beta, dst_m, work, work_complex)
 
 5879    integer, 
intent(in) :: p
 
 5880    real(dp), 
intent(in) :: c(3), src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
 
 5881        & alpha, src_m((p+1)*(p+1)), beta
 
 5883    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5885    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 5886    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 5888    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 5891    real(dp), 
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
 
 5893    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 5904    tmp_m(1:m) => work(n+1:n+m) 
 
 5906    tmp_m2(1:m) => work(n+1:n+m) 
 
 5909    vcos => work(n+1:n+m) 
 
 5911    vsin => work(n+1:n+m) 
 
 5913    call trgev(cphi, sphi, p, vcos, vsin)
 
 5921        & tmp_m2, zero, tmp_m, work, work_complex)
 
 5954    & alpha, src_m, beta, dst_m)
 
 5957    integer, 
intent(in) :: p
 
 5958    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 5959        & alpha, src_m((p+1)*(p+1)), beta, kappa
 
 5961    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 5963    real(dp) :: work(6*p*p + 19*p + 8), src_sk(p+1), dst_sk(p+1), s1, s2, ck(3)
 
 5964    complex(dp) :: work_complex(2*p+1), z1, z2
 
 5969    s1 = sqrt(2/(pi*real(z1)))
 
 5970    s2 = sqrt(2/(pi*real(z2)))
 
 5971    call cbesk(z1, pt5, 1, 1, work_complex(1), nz, ierr)
 
 5972    src_sk(1) = s1 * real(work_complex(1))
 
 5973    call cbesk(z2, pt5, 1, 1, work_complex(1), nz, ierr)
 
 5974    dst_sk(1) = s2 * real(work_complex(1))
 
 5976        call cbesk(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 5977        src_sk(2:p+1) = s1 * real(work_complex(2:p+1))
 
 5978        call cbesk(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 5979        dst_sk(2:p+1) = s2 * real(work_complex(2:p+1))
 
 5984        & alpha, src_m, beta, dst_m, work, work_complex)
 
 6013        & src_m, beta, dst_m, work, work_complex)
 
 6015    integer, 
intent(in) :: p
 
 6016    real(dp), 
intent(in) :: c(3), src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
 
 6017        & alpha, src_m((p+1)*(p+1)), beta
 
 6019    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 6021    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 6022    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 6024    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 6027    real(dp), 
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
 
 6029    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 6040    tmp_m(1:m) => work(n+1:n+m) 
 
 6042    tmp_m2(1:m) => work(n+1:n+m) 
 
 6045    vcos => work(n+1:n+m) 
 
 6047    vsin => work(n+1:n+m) 
 
 6049    call trgev(cphi, sphi, p, vcos, vsin)
 
 6057        & tmp_m2, zero, tmp_m, work, work_complex)
 
 6089        & src_m, beta, dst_m)
 
 6091    integer, 
intent(in) :: p
 
 6092    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 6093        & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
 
 6095    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 6097    real(dp) :: work(6*p*p + 19*p + 8)
 
 6100        & src_m, beta, dst_m, work)
 
 6128        & alpha, src_m, beta, dst_m, work)
 
 6130    integer, 
intent(in) :: p
 
 6131    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 6132        & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
 
 6134    real(dp), 
intent(inout) :: dst_m((p+1)*(p+1))
 
 6136    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 6138    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 6141    real(dp), 
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
 
 6143    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 6145    if (stheta .eq. 0) 
then 
 6148            & vcnk, alpha, src_m, beta, dst_m, work)
 
 6154    tmp_m(1:m) => work(n+1:n+m) 
 
 6156    tmp_m2(1:m) => work(n+1:n+m) 
 
 6159    vcos => work(n+1:n+m) 
 
 6161    vsin => work(n+1:n+m) 
 
 6163    call trgev(cphi, sphi, p, vcos, vsin)
 
 6171        & one, tmp_m2, zero, tmp_m, work)
 
 6179subroutine fmm_m2m_bessel_grad(p, sph_sk, vscales, sph_m, sph_m_grad)
 
 6181    integer, 
intent(in) :: p
 
 6182    real(dp), 
intent(in) :: sph_sk(p+2), vscales((p+2)**2), &
 
 6185    real(dp), 
intent(out) :: sph_m_grad((p+2)**2, 3)
 
 6187    real(dp), 
dimension(3, 3) :: zx_coord_transform, zy_coord_transform
 
 6188    real(dp) :: tmp((p+1)**2), tmp2((p+2)**2)
 
 6190    zx_coord_transform = zero
 
 6191    zx_coord_transform(3, 2) = one
 
 6192    zx_coord_transform(2, 3) = one
 
 6193    zx_coord_transform(1, 1) = one
 
 6194    zy_coord_transform = zero
 
 6195    zy_coord_transform(1, 2) = one
 
 6196    zy_coord_transform(2, 1) = one
 
 6197    zy_coord_transform(3, 3) = one
 
 6202        & one, tmp, zero, tmp2)
 
 6210        & one, tmp, zero, tmp2)
 
 6216        & one, sph_m, zero, sph_m_grad(:, 3))
 
 6217end subroutine fmm_m2m_bessel_grad
 
 6242        & src_l, beta, dst_l)
 
 6244    integer, 
intent(in) :: p
 
 6245    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
 
 6246        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
 
 6248    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 6250    real(dp) :: work(2*(p+1))
 
 6253        & src_l, beta, dst_l, work)
 
 6280        & src_l, beta, dst_l, work)
 
 6282    integer, 
intent(in) :: p
 
 6283    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
 
 6284        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
 
 6286    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 6288    real(dp), 
intent(out), 
target :: work(2*(p+1))
 
 6290    real(dp) :: r1, r2, tmp1, tmp2
 
 6291    integer :: j, k, n, indj, indn
 
 6293    real(dp), 
pointer :: pow_r1(:), pow_r2(:)
 
 6295    if (beta .eq. zero) 
then 
 6298        dst_l = beta * dst_l
 
 6303        pow_r1(1:p+1) => work(1:p+1)
 
 6304        pow_r2(1:p+1) => work(p+2:2*(p+1))
 
 6311            pow_r1(j) = pow_r1(j-1) * r1
 
 6312            pow_r2(j) = pow_r2(j-1) * r2
 
 6318                tmp1 = alpha * pow_r2(j+1) / vfact(j-k+1) / vfact(j+k+1) * &
 
 6322                    tmp2 = tmp1 * pow_r1(n+1) / vscales(indn) * &
 
 6323                        & vfact(n-k+1) * vfact(n+k+1) / vfact(n-j+1) / &
 
 6325                    if (mod(n+j, 2) .eq. 1) 
then 
 6329                        dst_l(indj) = dst_l(indj) + tmp2*src_l(indn)
 
 6331                        dst_l(indj+k) = dst_l(indj+k) + tmp2*src_l(indn+k)
 
 6332                        dst_l(indj-k) = dst_l(indj-k) + tmp2*src_l(indn-k)
 
 6343            do k = indj-j, indj+j
 
 6344                dst_l(k) = dst_l(k) + src_l(k)*tmp1
 
 6373        & src_l, beta, dst_l)
 
 6375    integer, 
intent(in) :: p
 
 6376    real(dp), 
intent(in) :: z, src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 6377        & alpha, src_l((p+1)*(p+1)), beta
 
 6379    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 6381    real(dp) :: work(2*p+1)
 
 6382    complex(dp) :: work_complex(2*p+1)
 
 6385        & alpha, src_l, beta, dst_l, work, work_complex)
 
 6412        & alpha, src_l, beta, dst_l, work, work_complex)
 
 6415    integer, 
intent(in) :: p
 
 6416    real(dp), 
intent(in) :: z, src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 6417        & alpha, src_l((p+1)*(p+1)), beta
 
 6419    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 6421    real(dp), 
intent(out), 
target :: work(2*p+1)
 
 6422    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 6424    real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
 
 6425    integer :: j, k, n, indj, indn, l, NZ, ierr
 
 6427    real(dp) :: fact(2*p+1)
 
 6428    complex(dp) :: complex_argument
 
 6430    if (alpha .eq. zero) 
then 
 6431        if (beta .eq. zero) 
then 
 6434            dst_l = beta * dst_l
 
 6441        fact(j+1) = dble(j) * fact(j)
 
 6446        complex_argument = z
 
 6447        call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 6448        work(1) = sqrt(pi/(2*z)) * real(work_complex(1))
 
 6450            call cbesi(complex_argument, 1.5d0, 1, 2*p, &
 
 6451                & work_complex(2:2*p+1), nz, ierr)
 
 6452            work(2:2*p+1) = sqrt(pi/(2*z)) * real(work_complex(2:2*p+1))
 
 6456        if (beta .eq. zero) 
then 
 6467                    tmp1 = vscales(indn) * ((-one)**(j+n)) * &
 
 6468                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 6469                        & src_si(n+1) / vscales(indj) * dst_si(j+1) * alpha
 
 6472                        tmp2 = (two**(-l)) / &
 
 6473                            & fact(l+1)*fact(2*l+1)/ &
 
 6474                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 6475                            & (work(j+n-l+1) / (z**l))
 
 6479                    res1 = res1 + tmp3*src_l(indn)
 
 6489                        tmp1 = vscales(indn+k) * ((-one)**(j+n)) * &
 
 6490                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 6491                            & src_si(n+1) / vscales(indj+k) * dst_si(j+1) * alpha
 
 6494                            tmp2 = (two**(-l)) / &
 
 6495                                & fact(l+k+1)*fact(2*l+1)/ &
 
 6496                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 6497                                & (work(j+n-l+1) / (z**l))
 
 6501                        res1 = res1 + tmp3*src_l(indn+k)
 
 6502                        res2 = res2 + tmp3*src_l(indn-k)
 
 6504                    dst_l(indj+k) = res1
 
 6505                    dst_l(indj-k) = res2
 
 6519                    tmp1 = vscales(indn) * ((-one)**(j+n)) * &
 
 6520                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 6521                        & src_si(n+1) / vscales(indj) * dst_si(j+1) * alpha
 
 6524                        tmp2 = (two**(-l)) / &
 
 6525                            & fact(l+1)*fact(2*l+1)/ &
 
 6526                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 6527                            & (work(j+n-l+1) / (z**l))
 
 6531                    res1 = res1 + tmp3*src_l(indn)
 
 6533                dst_l(indj) = beta*dst_l(indj) + res1
 
 6541                        tmp1 = vscales(indn+k) * ((-one)**(j+n)) * &
 
 6542                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 6543                            & src_si(n+1) / vscales(indj+k) * dst_si(j+1) * alpha
 
 6546                            tmp2 = (two**(-l)) / &
 
 6547                                & fact(l+k+1)*fact(2*l+1)/ &
 
 6548                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 6549                                & (work(j+n-l+1) / (z**l))
 
 6553                        res1 = res1 + tmp3*src_l(indn+k)
 
 6554                        res2 = res2 + tmp3*src_l(indn-k)
 
 6556                    dst_l(indj+k) = beta*dst_l(indj+k) + res1
 
 6557                    dst_l(indj-k) = beta*dst_l(indj-k) + res2
 
 6565        if (beta .eq. zero) 
then 
 6570                do k = indj-j, indj+j
 
 6571                    dst_l(k) = tmp1 * src_l(k)
 
 6581                do k = indj-j, indj+j
 
 6582                    dst_l(k) = beta*dst_l(k) + tmp1*src_l(k)
 
 6612        & src_l, beta, dst_l)
 
 6614    integer, 
intent(in) :: p
 
 6615    real(dp), 
intent(in) :: z, src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 6616        & alpha, src_l((p+1)*(p+1)), beta
 
 6618    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 6620    real(dp) :: work(2*p+1)
 
 6621    complex(dp) :: work_complex(2*p+1)
 
 6624        & alpha, src_l, beta, dst_l, work, work_complex)
 
 6651        & alpha, src_l, beta, dst_l, work, work_complex)
 
 6654    integer, 
intent(in) :: p
 
 6655    real(dp), 
intent(in) :: z, src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 6656        & alpha, src_l((p+1)*(p+1)), beta
 
 6658    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 6660    real(dp), 
intent(out), 
target :: work(2*p+1)
 
 6661    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 6663    real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
 
 6664    integer :: j, k, n, indj, indn, l, NZ, ierr
 
 6666    real(dp) :: fact(2*p+1)
 
 6667    complex(dp) :: complex_argument
 
 6669    if (alpha .eq. zero) 
then 
 6670        if (beta .eq. zero) 
then 
 6673            dst_l = beta * dst_l
 
 6680        fact(j+1) = dble(j) * fact(j)
 
 6685        complex_argument = z
 
 6686        call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 6687        work(1) = sqrt(pi/(2*z)) * real(work_complex(1))
 
 6689            call cbesi(complex_argument, 1.5d0, 1, 2*p, &
 
 6690                & work_complex(2:2*p+1), nz, ierr)
 
 6691            work(2:2*p+1) = sqrt(pi/(2*z)) * real(work_complex(2:2*p+1))
 
 6695        if (beta .eq. zero) 
then 
 6706                    tmp1 = vscales(indn) * ((-one)**(j+n)) * &
 
 6707                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 6708                        & src_si(n+1) / vscales(indj) * dst_si(j+1) * alpha
 
 6711                        tmp2 = (two**(-l)) / &
 
 6712                            & fact(l+1)*fact(2*l+1)/ &
 
 6713                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 6714                            & (work(j+n-l+1) / (z**l))
 
 6718                    res1 = res1 + tmp3*src_l(indj)
 
 6728                        tmp1 = vscales(indn+k) * ((-one)**(j+n)) * &
 
 6729                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 6730                            & src_si(n+1) / vscales(indj+k) * dst_si(j+1) * alpha
 
 6733                            tmp2 = (two**(-l)) / &
 
 6734                                & fact(l+k+1)*fact(2*l+1)/ &
 
 6735                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 6736                                & (work(j+n-l+1) / (z**l))
 
 6740                        res1 = res1 + tmp3*src_l(indj+k)
 
 6741                        res2 = res2 + tmp3*src_l(indj-k)
 
 6743                    dst_l(indn+k) = res1
 
 6744                    dst_l(indn-k) = res2
 
 6758                    tmp1 = vscales(indn) * ((-one)**(j+n)) * &
 
 6759                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 6761                        & src_si(n+1) / vscales(indj) * dst_si(j+1) * alpha
 
 6764                        tmp2 = (two**(-l)) / &
 
 6765                            & fact(l+1)*fact(2*l+1)/ &
 
 6766                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 6767                            & (work(j+n-l+1) / (z**l))
 
 6771                    res1 = res1 + tmp3*src_l(indj)
 
 6773                dst_l(indn) = beta*dst_l(indn) + res1
 
 6781                        tmp1 = vscales(indn+k) * ((-one)**(j+n)) * &
 
 6782                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 6783                            & src_si(n+1) / vscales(indj+k) * dst_si(j+1) * alpha
 
 6786                            tmp2 = (two**(-l)) / &
 
 6787                                & fact(l+k+1)*fact(2*l+1)/ &
 
 6788                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 6789                                & (work(j+n-l+1) / (z**l))
 
 6793                        res1 = res1 + tmp3*src_l(indj+k)
 
 6794                        res2 = res2 + tmp3*src_l(indj-k)
 
 6796                    dst_l(indn+k) = beta*dst_l(indn+k) + res1
 
 6797                    dst_l(indn-k) = beta*dst_l(indn-k) + res2
 
 6805        if (beta .eq. zero) 
then 
 6810                do k = indj-j, indj+j
 
 6811                    dst_l(k) = tmp1 * src_l(k)
 
 6821                do k = indj-j, indj+j
 
 6822                    dst_l(k) = beta*dst_l(k) + tmp1*src_l(k)
 
 6854        & alpha, src_l, beta, dst_l)
 
 6857    integer, 
intent(in) :: p
 
 6858    real(dp), 
intent(in) :: src_si(p+2), vscales((p+2)*(p+2)), &
 
 6859        & alpha, src_l((p+1)*(p+1)), beta
 
 6861    real(dp), 
intent(inout) :: dst_l((p+2)*(p+2))
 
 6864        & -alpha, src_l, beta, dst_l)
 
 6890        & alpha, src_l, beta, dst_l)
 
 6892    integer, 
intent(in) :: p
 
 6893    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
 
 6894        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
 
 6896    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 6898    real(dp) :: work(2*(p+1))
 
 6901        & alpha, src_l, beta, dst_l, work)
 
 6928        & alpha, src_l, beta, dst_l, work)
 
 6930    integer, 
intent(in) :: p
 
 6931    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
 
 6932        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
 
 6934    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 6936    real(dp), 
intent(out), 
target :: work(2*(p+1))
 
 6938    real(dp) :: r1, r2, tmp1, tmp2
 
 6939    integer :: j, k, n, indj, indn
 
 6941    real(dp), 
pointer :: pow_r1(:), pow_r2(:)
 
 6943    if (beta .eq. zero) 
then 
 6946        dst_l = beta * dst_l
 
 6951        pow_r1(1:p+1) => work(1:p+1)
 
 6952        pow_r2(1:p+1) => work(p+2:2*(p+1))
 
 6959            pow_r1(j) = pow_r1(j-1) * r1
 
 6960            pow_r2(j) = pow_r2(j-1) * r2
 
 6965                tmp1 = alpha * pow_r2(j+1) / vfact(j-k+1) / vfact(j+k+1) * &
 
 6969                    tmp2 = tmp1 * pow_r1(n+1) / vscales(indn) * &
 
 6970                        & vfact(n-k+1) * vfact(n+k+1) / vfact(n-j+1) / &
 
 6972                    if (mod(n+j, 2) .eq. 1) 
then 
 6976                        dst_l(indn) = dst_l(indn) + tmp2*src_l(indj)
 
 6978                        dst_l(indn+k) = dst_l(indn+k) + tmp2*src_l(indj+k)
 
 6979                        dst_l(indn-k) = dst_l(indn-k) + tmp2*src_l(indj-k)
 
 6990            do k = indj-j, indj+j
 
 6991                dst_l(k) = dst_l(k) + src_l(k)*tmp1
 
 7019    integer, 
intent(in) :: p
 
 7020    real(dp), 
intent(in) :: src_r, dst_r, alpha, src_l((p+1)*(p+1)), beta
 
 7022    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7024    real(dp) :: r1, tmp1
 
 7025    integer :: j, k, indj
 
 7027    if (beta .eq. zero) 
then 
 7030        dst_l = beta * dst_l
 
 7036        do k = indj-j, indj+j
 
 7037            dst_l(k) = dst_l(k) + src_l(k)*tmp1
 
 7064    integer, 
intent(in) :: p
 
 7065    real(dp), 
intent(in) :: src_r, dst_r, alpha, src_l((p+1)*(p+1)), beta
 
 7067    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7069    real(dp) :: r1, tmp1
 
 7070    integer :: j, k, indj
 
 7072    if (beta .eq. zero) 
then 
 7075        dst_l = beta * dst_l
 
 7081        do k = indj-j, indj+j
 
 7082            dst_l(k) = dst_l(k) + src_l(k)*tmp1
 
 7114        & src_l, beta, dst_l)
 
 7116    integer, 
intent(in) :: p
 
 7117    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 7118        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
 
 7120    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7122    real(dp) :: work(6*p*p+19*p+8)
 
 7125        & src_l, beta, dst_l, work)
 
 7155        & src_l, beta, dst_l, work)
 
 7157    integer, 
intent(in) :: p
 
 7158    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 7159        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
 
 7161    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7163    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 7165    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 7168    real(dp), 
pointer :: tmp_l(:), tmp_l2(:), vcos(:), vsin(:)
 
 7170    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 7172    if (stheta .eq. zero) 
then 
 7175            & alpha, src_l, beta, dst_l, work)
 
 7181    tmp_l(1:m) => work(n+1:n+m) 
 
 7183    tmp_l2(1:m) => work(n+1:n+m) 
 
 7186    vcos => work(n+1:n+m) 
 
 7188    vsin => work(n+1:n+m) 
 
 7190    call trgev(cphi, sphi, p, vcos, vsin)
 
 7198        & tmp_l2, zero, tmp_l, work)
 
 7231        & src_l, beta, dst_l)
 
 7234    integer, 
intent(in) :: p
 
 7235    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 7236        & alpha, src_l((p+1)*(p+1)), beta, kappa
 
 7238    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7240    real(dp) :: work(6*p*p + 19*p + 8), src_si(p+1), dst_si(p+1), s1, s2
 
 7241    complex(dp) :: work_complex(2*p+1), z1, z2
 
 7246    s1 = sqrt(pi/(2*real(z1)))
 
 7247    s2 = sqrt(pi/(2*real(z2)))
 
 7248    call cbesi(z1, pt5, 1, 1, work_complex(1), nz, ierr)
 
 7249    src_si(1) = s1 * real(work_complex(1))
 
 7250    call cbesi(z2, pt5, 1, 1, work_complex(1), nz, ierr)
 
 7251    dst_si(1) = s2 * real(work_complex(1))
 
 7253        call cbesi(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 7254        src_si(2:p+1) = s1 * real(work_complex(2:p+1))
 
 7255        call cbesi(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 7256        dst_si(2:p+1) = s2 * real(work_complex(2:p+1))
 
 7260        & src_l, beta, dst_l, work, work_complex)
 
 7289        & src_l, beta, dst_l, work, work_complex)
 
 7291    integer, 
intent(in) :: p
 
 7292    real(dp), 
intent(in) :: c(3), src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 7293        & alpha, src_l((p+1)*(p+1)), beta
 
 7295    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7297    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 7298    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 7300    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 7303    real(dp), 
pointer :: tmp_l(:), tmp_l2(:), vcos(:), vsin(:)
 
 7305    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 7315    tmp_l(1:m) => work(n+1:n+m) 
 
 7317    tmp_l2(1:m) => work(n+1:n+m) 
 
 7320    vcos => work(n+1:n+m) 
 
 7322    vsin => work(n+1:n+m) 
 
 7324    call trgev(cphi, sphi, p, vcos, vsin)
 
 7332        & one, tmp_l2, zero, tmp_l, work, work_complex)
 
 7365        & src_l, beta, dst_l)
 
 7368    integer, 
intent(in) :: p
 
 7369    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 7370        & alpha, src_l((p+1)*(p+1)), beta, kappa
 
 7372    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7374    real(dp) :: work(6*p*p + 19*p + 8), src_si(p+1), dst_si(p+1), s1, s2
 
 7375    complex(dp) :: work_complex(2*p+1), z1, z2
 
 7380    s1 = sqrt(pi/(2*real(z1)))
 
 7381    s2 = sqrt(pi/(2*real(z2)))
 
 7382    call cbesi(z1, pt5, 1, 1, work_complex(1), nz, ierr)
 
 7383    src_si(1) = s1 * real(work_complex(1))
 
 7384    call cbesi(z2, pt5, 1, 1, work_complex(1), nz, ierr)
 
 7385    dst_si(1) = s2 * real(work_complex(1))
 
 7387        call cbesi(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 7388        src_si(2:p+1) = s1 * real(work_complex(2:p+1))
 
 7389        call cbesi(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 7390        dst_si(2:p+1) = s2 * real(work_complex(2:p+1))
 
 7394        & src_l, beta, dst_l, work, work_complex)
 
 7423        & alpha, src_l, beta, dst_l, work, work_complex)
 
 7425    integer, 
intent(in) :: p
 
 7426    real(dp), 
intent(in) :: c(3), src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 7427        & alpha, src_l((p+1)*(p+1)), beta
 
 7429    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7431    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 7432    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 7434    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 7437    real(dp), 
pointer :: tmp_l(:), tmp_l2(:), vcos(:), vsin(:)
 
 7439    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 7450    tmp_l(1:m) => work(n+1:n+m) 
 
 7452    tmp_l2(1:m) => work(n+1:n+m) 
 
 7455    vcos => work(n+1:n+m) 
 
 7457    vsin => work(n+1:n+m) 
 
 7459    call trgev(cphi, sphi, p, vcos, vsin)
 
 7467        & one, tmp_l2, zero, tmp_l, work, work_complex)
 
 7475subroutine fmm_l2l_bessel_grad(p, sph_si, vscales, sph_l, sph_l_grad)
 
 7477    integer, 
intent(in) :: p
 
 7478    real(dp), 
intent(in) :: sph_si(p+2), vscales((p+2)**2), sph_l((p+1)**2)
 
 7480    real(dp), 
intent(out) :: sph_l_grad((p+2)**2, 3)
 
 7482    real(dp), 
dimension(3, 3) :: zx_coord_transform, zy_coord_transform
 
 7483    real(dp) :: tmp((p+1)**2), tmp2((p+2)**2)
 
 7485    zx_coord_transform = zero
 
 7486    zx_coord_transform(3, 2) = one
 
 7487    zx_coord_transform(2, 3) = one
 
 7488    zx_coord_transform(1, 1) = one
 
 7489    zy_coord_transform = zero
 
 7490    zy_coord_transform(1, 2) = one
 
 7491    zy_coord_transform(2, 1) = one
 
 7492    zy_coord_transform(3, 3) = one
 
 7497        & one, tmp, zero, tmp2)
 
 7505        & one, tmp, zero, tmp2)
 
 7511        & one, sph_l, zero, sph_l_grad(:, 3))
 
 7512end subroutine fmm_l2l_bessel_grad
 
 7540        & alpha, src_l, beta, dst_l)
 
 7542    integer, 
intent(in) :: p
 
 7543    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 7544        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
 
 7546    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7548    real(dp) :: work(6*p*p + 19*p + 8)
 
 7551        & src_l, beta, dst_l, work)
 
 7581        & alpha, src_l, beta, dst_l, work)
 
 7583    integer, 
intent(in) :: p
 
 7584    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 7585        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
 
 7587    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 7589    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 7591    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 7594    real(dp), 
pointer :: tmp_l(:), tmp_l2(:), vcos(:), vsin(:)
 
 7596    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 7598    if (stheta .eq. zero) 
then 
 7601            & vfact, alpha, src_l, beta, dst_l, work)
 
 7607    tmp_l(1:m) => work(n+1:n+m) 
 
 7609    tmp_l2(1:m) => work(n+1:n+m) 
 
 7612    vcos => work(n+1:n+m) 
 
 7614    vsin => work(n+1:n+m) 
 
 7616    call trgev(cphi, sphi, p, vcos, vsin)
 
 7624        & one, tmp_l2, zero, tmp_l, work)
 
 7656        & m2l_ztranslate_coef, alpha, src_m, beta, dst_l)
 
 7658    integer, 
intent(in) :: pm, pl
 
 7659    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((pm+pl+1)*(pm+pl+1)), &
 
 7660        & m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
 
 7663    real(dp), 
intent(inout) :: dst_l((pl+1)*(pl+1))
 
 7665    real(dp) :: work((pm+2)*(pm+1))
 
 7668        & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
 
 7696        & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
 
 7698    integer, 
intent(in) :: pm, pl
 
 7699    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((pm+pl+1)*(pm+pl+1)), &
 
 7700        & m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
 
 7703    real(dp), 
intent(inout) :: dst_l((pl+1)*(pl+1))
 
 7705    real(dp), 
intent(out), 
target :: work((pm+2)*(pm+1))
 
 7707    real(dp) :: tmp1, r1, r2, res1, res2, pow_r2
 
 7708    integer :: j, k, n, indj, indk1, indk2
 
 7710    real(dp), 
pointer :: src_m2(:), pow_r1(:)
 
 7712    if (alpha .eq. zero) 
then 
 7713        if (beta .eq. zero) 
then 
 7716            dst_l = beta * dst_l
 
 7723    if (z .eq. zero) 
then 
 7728    src_m2(1:n) => work(1:n)
 
 7729    pow_r1(1:pm+1) => work(n+1:n+pm+1)
 
 7738        pow_r1(j+1) = tmp1 / vscales(indj)
 
 7746        src_m2(j+1) = pow_r1(j+1) * src_m(indj)
 
 7756            src_m2(indk1+j-k) = pow_r1(j+1) * src_m(indj+k)
 
 7757            src_m2(indk2+j-k) = pow_r1(j+1) * src_m(indj-k)
 
 7763    if (beta .eq. zero) 
then 
 7768            tmp1 = alpha * vscales(indj) * pow_r2
 
 7769            pow_r2 = pow_r2 * r2
 
 7773                res1 = res1 + m2l_ztranslate_coef(n+1, 1, j+1)*src_m2(n+1)
 
 7775            dst_l(indj) = tmp1 * res1
 
 7780                indk1 = pm + 2 + (2*pm-k+2)*(k-1)
 
 7781                indk2 = indk1 + pm - k + 1
 
 7787                        & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)* &
 
 7790                        & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)* &
 
 7793                dst_l(indj+k) = tmp1 * res1
 
 7794                dst_l(indj-k) = tmp1 * res2
 
 7802            tmp1 = alpha * vscales(indj) * pow_r2
 
 7803            pow_r2 = pow_r2 * r2
 
 7806                res1 = res1 + m2l_ztranslate_coef(n+1, 1, j+1)*src_m2(n+1)
 
 7808            dst_l(indj) = beta*dst_l(indj) + tmp1*res1
 
 7812                indk1 = pm + 2 + (2*pm-k+2)*(k-1)
 
 7813                indk2 = indk1 + pm - k + 1
 
 7818                        & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)* &
 
 7821                        & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)* &
 
 7824                dst_l(indj+k) = beta*dst_l(indj+k) + tmp1*res1
 
 7825                dst_l(indj-k) = beta*dst_l(indj-k) + tmp1*res2
 
 7855        & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
 
 7857    integer, 
intent(in) :: pl, pm
 
 7858    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((pm+pl+1)*(pm+pl+1)), &
 
 7859        & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1), alpha, &
 
 7860        & src_l((pl+1)*(pl+1)), beta
 
 7862    real(dp), 
intent(inout) :: dst_m((pm+1)*(pm+1))
 
 7864    real(dp) :: work((pl+2)*(pl+1))
 
 7867        & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
 
 7894        & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
 
 7896    integer, 
intent(in) :: pl, pm
 
 7897    real(dp), 
intent(in) :: z, src_r, dst_r, vscales((pm+pl+1)*(pm+pl+1)), &
 
 7898        & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1), alpha, &
 
 7899        & src_l((pl+1)*(pl+1)), beta
 
 7901    real(dp), 
intent(inout) :: dst_m((pm+1)*(pm+1))
 
 7903    real(dp), 
intent(out), 
target :: work((pl+2)*(pl+1))
 
 7905    real(dp) :: tmp1, r1, r2, pow_r1, res1, res2
 
 7906    integer :: j, k, n, indj, indn, indk1, indk2
 
 7908    real(dp), 
pointer :: pow_r2(:), src_l2(:)
 
 7910    if (alpha .eq. zero) 
then 
 7911        if (beta .eq. zero) 
then 
 7914            dst_m = beta * dst_m
 
 7921    if (z .eq. zero) 
then 
 7926    src_l2(1:n) => work(1:n)
 
 7927    pow_r2(1:pl+1) => work(n+1:n+pl+1)
 
 7936        pow_r2(j+1) = tmp1 * vscales(indj)
 
 7944        src_l2(j+1) = pow_r2(j+1) * src_l(indj)
 
 7953            src_l2(indk1+j-k) = pow_r2(j+1) * src_l(indj+k)
 
 7954            src_l2(indk2+j-k) = pow_r2(j+1) * src_l(indj-k)
 
 7960    if (beta .eq. zero) 
then 
 7964            tmp1 = alpha * pow_r1 / vscales(indn)
 
 7965            pow_r1 = pow_r1 * r1
 
 7969                res1 = res1 + m2l_ztranslate_adj_coef(j+1, 1, n+1)*src_l2(j+1)
 
 7971            dst_m(indn) = tmp1 * res1
 
 7975                indk1 = pl + 2 + (2*pl-k+2)*(k-1)
 
 7976                indk2 = indk1 + pl - k + 1
 
 7982                        & m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1)* &
 
 7985                        & m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1)* &
 
 7988                dst_m(indn+k) = tmp1 * res1
 
 7989                dst_m(indn-k) = tmp1 * res2
 
 7997            tmp1 = alpha * pow_r1 / vscales(indn)
 
 7998            pow_r1 = pow_r1 * r1
 
 8002                res1 = res1 + m2l_ztranslate_adj_coef(j+1, 1, n+1)*src_l2(j+1)
 
 8004            dst_m(indn) = beta*dst_m(indn) + tmp1*res1
 
 8008                indk1 = pl + 2 + (2*pl-k+2)*(k-1)
 
 8009                indk2 = indk1 + pl - k + 1
 
 8015                        & m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1)* &
 
 8018                        & m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1)* &
 
 8021                dst_m(indn+k) = beta*dst_m(indn+k) + tmp1*res1
 
 8022                dst_m(indn-k) = beta*dst_m(indn-k) + tmp1*res2
 
 8055        & m2l_ztranslate_coef, alpha, src_m, beta, dst_l)
 
 8057    integer, 
intent(in) :: pm, pl
 
 8058    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((pm+pl+1)**2), &
 
 8059        & m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
 
 8062    real(dp), 
intent(inout) :: dst_l((pl+1)*(pl+1))
 
 8064    real(dp) :: work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
 
 8067        & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
 
 8099        & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
 
 8101    integer, 
intent(in) :: pm, pl
 
 8102    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((pm+pl+1)**2), &
 
 8103        & m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
 
 8106    real(dp), 
intent(inout) :: dst_l((pl+1)*(pl+1))
 
 8108    real(dp), 
intent(out), 
target :: &
 
 8109        & work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
 
 8111    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 8114    real(dp), 
pointer :: tmp_ml(:), tmp_ml2(:), vcos(:), vsin(:)
 
 8116    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 8118    if (stheta .eq. zero) 
then 
 8121            & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
 
 8128    tmp_ml(1:m) => work(n+1:n+m) 
 
 8130    tmp_ml2(1:m) => work(n+1:n+m) 
 
 8133    vcos => work(n+1:n+m) 
 
 8135    vsin => work(n+1:n+m) 
 
 8137    call trgev(cphi, sphi, p, vcos, vsin)
 
 8145        & m2l_ztranslate_coef, one, tmp_ml2, zero, tmp_ml, work)
 
 8178        & src_m, beta, dst_l)
 
 8181    integer, 
intent(in) :: p
 
 8182    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 8183        & alpha, src_m((p+1)*(p+1)), beta, kappa
 
 8185    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 8187    real(dp) :: work(6*p*p + 19*p + 8), src_sk(p+1), dst_si(p+1), s1, s2
 
 8188    complex(dp) :: work_complex(2*p+1), z1, z2
 
 8193    s1 = sqrt(2/(pi*real(z1)))
 
 8194    s2 = sqrt(pi/(2*real(z2)))
 
 8195    call cbesk(z1, pt5, 1, 1, work_complex(1), nz, ierr)
 
 8196    src_sk(1) = s1 * real(work_complex(1))
 
 8197    call cbesi(z2, pt5, 1, 1, work_complex(1), nz, ierr)
 
 8198    dst_si(1) = s2 * real(work_complex(1))
 
 8200        call cbesk(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 8201        src_sk(2:p+1) = s1 * real(work_complex(2:p+1))
 
 8202        call cbesi(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 8203        dst_si(2:p+1) = s2 * real(work_complex(2:p+1))
 
 8207        & src_m, beta, dst_l, work, work_complex)
 
 8236        & src_m, beta, dst_l, work, work_complex)
 
 8238    integer, 
intent(in) :: p
 
 8239    real(dp), 
intent(in) :: c(3), src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 8240        & alpha, src_m((p+1)*(p+1)), beta
 
 8242    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 8244    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 8245    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 8247    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 8250    real(dp), 
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
 
 8252    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 8263    tmp_m(1:m) => work(n+1:n+m) 
 
 8265    tmp_m2(1:m) => work(n+1:n+m) 
 
 8268    vcos => work(n+1:n+m) 
 
 8270    vsin => work(n+1:n+m) 
 
 8272    call trgev(cphi, sphi, p, vcos, vsin)
 
 8280        & tmp_m2, zero, tmp_m, work, work_complex)
 
 8313        & src_m, beta, dst_l)
 
 8316    integer, 
intent(in) :: p
 
 8317    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
 
 8318        & alpha, src_m((p+1)*(p+1)), beta, kappa
 
 8320    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 8322    real(dp) :: work(6*p*p + 19*p + 8), src_sk(p+1), dst_si(p+1), s1, s2, ck(3)
 
 8323    complex(dp) :: work_complex(2*p+1), z1, z2
 
 8328    s1 = sqrt(2/(pi*real(z1)))
 
 8329    s2 = sqrt(pi/(2*real(z2)))
 
 8330    call cbesk(z1, pt5, 1, 1, work_complex(1), nz, ierr)
 
 8331    src_sk(1) = s1 * real(work_complex(1))
 
 8332    call cbesi(z2, pt5, 1, 1, work_complex(1), nz, ierr)
 
 8333    dst_si(1) = s2 * real(work_complex(1))
 
 8335        call cbesk(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 8336        src_sk(2:p+1) = s1 * real(work_complex(2:p+1))
 
 8337        call cbesi(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
 
 8338        dst_si(2:p+1) = s2 * real(work_complex(2:p+1))
 
 8343        & alpha, src_m, beta, dst_l, work, work_complex)
 
 8372        & src_m, beta, dst_l, work, work_complex)
 
 8374    integer, 
intent(in) :: p
 
 8375    real(dp), 
intent(in) :: c(3), src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 8376        & alpha, src_m((p+1)*(p+1)), beta
 
 8378    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 8380    real(dp), 
intent(out), 
target :: work(6*p*p + 19*p + 8)
 
 8381    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 8383    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 8386    real(dp), 
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
 
 8388    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 8399    tmp_m(1:m) => work(n+1:n+m) 
 
 8401    tmp_m2(1:m) => work(n+1:n+m) 
 
 8404    vcos => work(n+1:n+m) 
 
 8406    vsin => work(n+1:n+m) 
 
 8408    call trgev(cphi, sphi, p, vcos, vsin)
 
 8416        & tmp_m2, zero, tmp_m, work, work_complex)
 
 8446        & src_m, beta, dst_l)
 
 8448    integer, 
intent(in) :: p
 
 8449    real(dp), 
intent(in) :: z, src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 8450        & alpha, src_m((p+1)*(p+1)), beta
 
 8452    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 8454    real(dp) :: work(2*p+1)
 
 8455    complex(dp) :: work_complex(2*p+1)
 
 8458        & alpha, src_m, beta, dst_l, work, work_complex)
 
 8485        & alpha, src_m, beta, dst_l, work, work_complex)
 
 8488    integer, 
intent(in) :: p
 
 8489    real(dp), 
intent(in) :: z, src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 8490        & alpha, src_m((p+1)*(p+1)), beta
 
 8492    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 8494    real(dp), 
intent(out), 
target :: work(2*p+1)
 
 8495    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 8497    real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
 
 8498    integer :: j, k, n, indj, indn, l, NZ, ierr
 
 8500    real(dp) :: fact(2*p+1)
 
 8501    complex(dp) :: complex_argument
 
 8503    if (alpha .eq. zero) 
then 
 8504        if (beta .eq. zero) 
then 
 8507            dst_l = beta * dst_l
 
 8514        fact(j+1) = dble(j) * fact(j)
 
 8519        complex_argument = z
 
 8520        call cbesk(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 8521        work(1) = sqrt(2/(pi*z)) * real(work_complex(1))
 
 8523            call cbesk(complex_argument, 1.5d0, 1, 2*p, &
 
 8524                & work_complex(2:2*p+1), nz, ierr)
 
 8525            work(2:2*p+1) = sqrt(2/(pi*z)) * real(work_complex(2:2*p+1))
 
 8529        if (beta .eq. zero) 
then 
 8540                    tmp1 = vscales(indn) * ((-one)**(n)) * &
 
 8541                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 8542                        & src_sk(n+1) / vscales(indj) * dst_si(j+1)
 
 8545                        tmp2 = ((-two)**(-l)) / &
 
 8546                            & fact(l+1)*fact(2*l+1)/ &
 
 8547                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 8548                            & (work(j+n-l+1) / (z**l))
 
 8552                    res1 = res1 + tmp3*src_m(indn)
 
 8562                        tmp1 = vscales(indn+k) * ((-one)**(n)) * &
 
 8563                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 8564                            & src_sk(n+1) / vscales(indj+k) * dst_si(j+1)
 
 8567                            tmp2 = ((-two)**(-l)) / &
 
 8568                                & fact(l+k+1)*fact(2*l+1)/ &
 
 8569                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 8570                                & (work(j+n-l+1) / (z**l))
 
 8574                        res1 = res1 + tmp3*src_m(indn+k)
 
 8575                        res2 = res2 + tmp3*src_m(indn-k)
 
 8577                    dst_l(indj+k) = res1
 
 8578                    dst_l(indj-k) = res2
 
 8592                    tmp1 = vscales(indn) * ((-one)**(n)) * &
 
 8593                        & dble(2*j+1) * fact(j+1) * fact(n+1) / &
 
 8594                        & src_sk(n+1) / vscales(indj) * dst_si(j+1)
 
 8597                        tmp2 = ((-two)**(-l)) / &
 
 8598                            & fact(l+1)*fact(2*l+1)/ &
 
 8599                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 8600                            & (work(j+n-l+1) / (z**l))
 
 8604                    res1 = res1 + tmp3*src_m(indn)
 
 8606                dst_l(indj) = beta*dst_l(indj) + res1
 
 8614                        tmp1 = vscales(indn+k) * ((-one)**(n)) * &
 
 8615                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
 
 8616                            & src_sk(n+1) / vscales(indj+k) * dst_si(j+1)
 
 8619                            tmp2 = ((-two)**(-l)) / &
 
 8620                                & fact(l+k+1)*fact(2*l+1)/ &
 
 8621                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 8622                                & (work(j+n-l+1) / (z**l))
 
 8626                        res1 = res1 + tmp3*src_m(indn+k)
 
 8627                        res2 = res2 + tmp3*src_m(indn-k)
 
 8629                    dst_l(indj+k) = beta*dst_l(indj+k) + res1
 
 8630                    dst_l(indj-k) = beta*dst_l(indj-k) + res2
 
 8638        if (beta .eq. zero) 
then 
 8643                do k = indj-j, indj+j
 
 8644                    dst_l(k) = tmp1 * src_m(k)
 
 8654                do k = indj-j, indj+j
 
 8655                    dst_l(k) = beta*dst_l(k) + tmp1*src_m(k)
 
 8685        & src_m, beta, dst_l)
 
 8687    integer, 
intent(in) :: p
 
 8688    real(dp), 
intent(in) :: z, src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 8689        & alpha, src_m((p+1)*(p+1)), beta
 
 8691    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 8693    real(dp) :: work(2*p+1)
 
 8694    complex(dp) :: work_complex(2*p+1)
 
 8697        & alpha, src_m, beta, dst_l, work, work_complex)
 
 8724        & alpha, src_m, beta, dst_l, work, work_complex)
 
 8727    integer, 
intent(in) :: p
 
 8728    real(dp), 
intent(in) :: z, src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
 
 8729        & alpha, src_m((p+1)*(p+1)), beta
 
 8731    real(dp), 
intent(inout) :: dst_l((p+1)*(p+1))
 
 8733    real(dp), 
intent(out), 
target :: work(2*p+1)
 
 8734    complex(dp), 
intent(out) :: work_complex(2*p+1)
 
 8736    real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
 
 8737    integer :: j, k, n, indj, indn, l, NZ, ierr
 
 8739    real(dp) :: fact(2*p+1)
 
 8740    complex(dp) :: complex_argument
 
 8742    if (alpha .eq. zero) 
then 
 8743        if (beta .eq. zero) 
then 
 8746            dst_l = beta * dst_l
 
 8753        fact(j+1) = dble(j) * fact(j)
 
 8758        complex_argument = z
 
 8759        call cbesk(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
 
 8760        work(1) = sqrt(2/(pi*z)) * real(work_complex(1))
 
 8762            call cbesk(complex_argument, 1.5d0, 1, 2*p, &
 
 8763                & work_complex(2:2*p+1), nz, ierr)
 
 8764            work(2:2*p+1) = sqrt(2/(pi*z)) * real(work_complex(2:2*p+1))
 
 8768        if (beta .eq. zero) 
then 
 8779                    tmp1 = vscales(indn) * ((-one)**(n)) * &
 
 8780                        & dble(2*j+1) * fact(j+1) * fact(n+1) * &
 
 8781                        & src_sk(j+1) / vscales(indj) / dst_si(n+1)
 
 8784                        tmp2 = ((-two)**(-l)) / &
 
 8785                            & fact(l+1)*fact(2*l+1)/ &
 
 8786                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 8787                            & (work(j+n-l+1) / (z**l))
 
 8791                    res1 = res1 + tmp3*src_m(indj)
 
 8801                        tmp1 = vscales(indn+k) * ((-one)**(n)) * &
 
 8802                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) * &
 
 8803                            & src_sk(j+1) / vscales(indj+k) / dst_si(n+1)
 
 8806                            tmp2 = ((-two)**(-l)) / &
 
 8807                                & fact(l+k+1)*fact(2*l+1)/ &
 
 8808                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 8809                                & (work(j+n-l+1) / (z**l))
 
 8813                        res1 = res1 + tmp3*src_m(indj+k)
 
 8814                        res2 = res2 + tmp3*src_m(indj-k)
 
 8816                    dst_l(indn+k) = res1
 
 8817                    dst_l(indn-k) = res2
 
 8831                    tmp1 = vscales(indn) * ((-one)**(n)) * &
 
 8832                        & dble(2*j+1) * fact(j+1) * fact(n+1) * &
 
 8833                        & src_sk(j+1) / vscales(indj) / dst_si(n+1)
 
 8836                        tmp2 = ((-two)**(-l)) / &
 
 8837                            & fact(l+1)*fact(2*l+1)/ &
 
 8838                            & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
 
 8839                            & (work(j+n-l+1) / (z**l))
 
 8843                    res1 = res1 + tmp3*src_m(indj)
 
 8845                dst_l(indn) = beta*dst_l(indn) + res1
 
 8853                        tmp1 = vscales(indn+k) * ((-one)**(n)) * &
 
 8854                            & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) * &
 
 8855                            & src_sk(j+1) / vscales(indj+k) / dst_si(n+1)
 
 8858                            tmp2 = ((-two)**(-l)) / &
 
 8859                                & fact(l+k+1)*fact(2*l+1)/ &
 
 8860                                & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
 
 8861                                & (work(j+n-l+1) / (z**l))
 
 8865                        res1 = res1 + tmp3*src_m(indj+k)
 
 8866                        res2 = res2 + tmp3*src_m(indj-k)
 
 8868                    dst_l(indn+k) = beta*dst_l(indn+k) + res1
 
 8869                    dst_l(indn-k) = beta*dst_l(indn-k) + res2
 
 8877        if (beta .eq. zero) 
then 
 8882                do k = indj-j, indj+j
 
 8883                    dst_l(k) = tmp1 * src_m(k)
 
 8893                do k = indj-j, indj+j
 
 8894                    dst_l(k) = beta*dst_l(k) + tmp1*src_m(k)
 
 8929        & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
 
 8931    integer, 
intent(in) :: pl, pm
 
 8932    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((pm+pl+1)**2), &
 
 8933        & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1), alpha, &
 
 8934        & src_l((pl+1)*(pl+1)), beta
 
 8936    real(dp), 
intent(inout) :: dst_m((pm+1)*(pm+1))
 
 8938    real(dp) :: work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
 
 8941        & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
 
 8973        & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
 
 8975    integer, 
intent(in) :: pl, pm
 
 8976    real(dp), 
intent(in) :: c(3), src_r, dst_r, vscales((pm+pl+1)**2), &
 
 8977        & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1), alpha, &
 
 8978        & src_l((pl+1)*(pl+1)), beta
 
 8980    real(dp), 
intent(inout) :: dst_m((pm+1)*(pm+1))
 
 8982    real(dp), 
intent(out), 
target :: &
 
 8983        & work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
 
 8985    real(dp) :: rho, ctheta, stheta, cphi, sphi
 
 8988    real(dp), 
pointer :: tmp_ml(:), tmp_ml2(:), vcos(:), vsin(:)
 
 8990    call carttosph(c, rho, ctheta, stheta, cphi, sphi)
 
 8992    if (stheta .eq. 0) 
then 
 8995            & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
 
 9002    tmp_ml(1:m) => work(n+1:n+m) 
 
 9004    tmp_ml2(1:m) => work(n+1:n+m) 
 
 9007    vcos => work(n+1:n+m) 
 
 9009    vsin => work(n+1:n+m) 
 
 9011    call trgev(cphi, sphi, p, vcos, vsin)
 
 9019        & m2l_ztranslate_adj_coef, one, tmp_ml2, zero, tmp_ml, work)
 
Compile-time constants and definitions.
 
Harmonics-related core routines.
 
subroutine ylmscale(p, vscales, v4pi2lp1, vscales_rel)
Compute scaling factors of real normalized spherical harmonics.
 
subroutine fmm_m2p_adj(c, src_q, dst_r, p, vscales_rel, beta, dst_m)
Adjoint M2P operation.
 
subroutine fmm_l2p_adj_work(c, src_q, dst_r, p, vscales_rel, beta, dst_l, work)
Adjoint L2P.
 
subroutine fmm_m2m_ztranslate(z, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m)
Direct M2M translation over OZ axis.
 
subroutine fmm_l2l_scale_adj(src_r, dst_r, p, alpha, src_l, beta, dst_l)
Adjoint scale L2L, when spherical harmonics are centered in the same point.
 
subroutine fmm_m2m_bessel_rotation_adj_work(c, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_m2l_bessel_rotation(c, src_r, dst_r, kappa, p, vscales, alpha, src_m, beta, dst_l)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_m2l_rotation(c, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l)
Direct M2L translation by 4 rotations and 1 translation.
 
subroutine fmm_l2l_bessel_derivative_ztranslate_work(src_si, p, vscales, alpha, src_l, beta, dst_l)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2m_ztranslate_work(z, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Direct M2M translation over OZ axis.
 
subroutine fmm_sph_rotate_oxz_work(p, ctheta, stheta, alpha, src, beta, dst, work)
Transform spherical harmonics in the OXZ plane.
 
subroutine fmm_m2m_ztranslate_adj(z, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m)
Adjoint M2M translation over OZ axis.
 
subroutine fmm_m2m_bessel_derivative_ztranslate_work(src_sk, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation over OZ axis.
 
subroutine fmm_l2l_rotation_adj(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l)
Adjoint L2L translation by 4 rotations and 1 translation.
 
subroutine fmm_m2l_ztranslate(z, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l)
Direct M2L translation over OZ axis.
 
subroutine fmm_m2p_bessel_adj(c, src_q, dst_r, kappa, p, vscales, beta, dst_m)
Adjoint M2P operation.
 
subroutine fmm_m2m_rotation(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_l2l_scale(src_r, dst_r, p, alpha, src_l, beta, dst_l)
Scale L2L, when spherical harmonics are centered in the same point.
 
subroutine fmm_sph_rotate_oxz(p, ctheta, stheta, alpha, src, beta, dst)
Transform spherical harmonics in the OXZ plane.
 
subroutine fmm_m2l_ztranslate_adj_work(z, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
Adjoint M2L translation over OZ axis.
 
subroutine fmm_m2m_bessel_ztranslate(z, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation over OZ axis.
 
subroutine fmm_l2p(c, src_r, p, vscales, alpha, src_l, beta, dst_v)
Accumulate potential, induced by local spherical harmonics.
 
subroutine polleg_work(ctheta, stheta, p, vplm, work)
Compute associated Legendre polynomials.
 
subroutine fmm_l2l_ztranslate_work(z, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Direct L2L translation over OZ axis.
 
subroutine fmm_l2p_adj(c, src_q, dst_r, p, vscales_rel, beta, dst_l)
Adjoint of L2P.
 
subroutine fmm_m2l_bessel_ztranslate(z, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l)
Direct M2M translation over OZ axis.
 
subroutine fmm_p2m(c, src_q, dst_r, p, vscales, beta, dst_m)
Accumulate a multipole expansion induced by a particle of a given charge.
 
subroutine fmm_l2l_bessel_ztranslate_adj(z, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2m_scale(src_r, dst_r, p, alpha, src_m, beta, dst_m)
Scale M2M, when spherical harmonics are centered in the same point.
 
subroutine fmm_sph_rotate_oz_adj(p, vcos, vsin, alpha, src, beta, dst)
Rotate spherical harmonics around OZ axis in an opposite direction.
 
subroutine fmm_l2l_ztranslate_adj(z, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l)
Adjoint L2L translation over OZ axis.
 
subroutine fmm_l2l_bessel_ztranslate_work(z, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct M2M translation over OZ axis.
 
subroutine fmm_sph_rotate_oz_adj_work(p, vcos, vsin, alpha, src, beta, dst)
Rotate spherical harmonics around OZ axis in an opposite direction.
 
subroutine trgev(cphi, sphi, p, vcos, vsin)
Compute arrays of  and .
 
subroutine fmm_m2p_adj_work(c, src_q, dst_r, p, vscales_rel, beta, dst_m, work)
Adjoint M2P operation.
 
subroutine carttosph(x, rho, ctheta, stheta, cphi, sphi)
Convert input cartesian coordinate into spherical coordinate.
 
subroutine fmm_l2l_bessel_rotation_adj_work(c, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct L2L translation by 4 rotations and 1 translation.
 
subroutine fmm_m2p_bessel_adj_work(c, src_q, dst_sk, p, vscales, beta, dst_m, work_complex, work)
Adjoint M2P operation.
 
subroutine fmm_m2m_bessel_rotation_work(c, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine ylmbas(x, rho, ctheta, stheta, cphi, sphi, p, vscales, vylm, vplm, vcos, vsin)
Compute all spherical harmonics up to a given degree at a given point.
 
subroutine fmm_l2l_rotation_work(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Direct L2L translation by 4 rotations and 1 translation.
 
subroutine fmm_m2l_ztranslate_work(z, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
Direct M2L translation over OZ axis.
 
subroutine polleg(ctheta, stheta, p, vplm)
Compute associated Legendre polynomials.
 
subroutine fmm_l2l_bessel_ztranslate_adj_work(z, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2p_work(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v, work)
Accumulate potential, induced by multipole spherical harmonics.
 
subroutine fmm_l2l_rotation(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l)
Direct L2L translation by 4 rotations and 1 translation.
 
subroutine fmm_l2l_bessel_rotation(c, src_r, dst_r, kappa, p, vscales, alpha, src_l, beta, dst_l)
Direct L2L translation by 4 rotations and 1 translation.
 
subroutine fmm_m2l_rotation_adj(c, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
Adjoint M2L translation by 4 rotations and 1 translation.
 
subroutine fmm_sph_rotate_oz(p, vcos, vsin, alpha, src, beta, dst)
Rotate spherical harmonics around OZ axis.
 
subroutine fmm_m2m_rotation_work(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_sph_transform(p, r1, alpha, src, beta, dst)
Transform coefficients of spherical harmonics to a new cartesion system.
 
subroutine fmm_l2p_bessel(c, src_r, p, vscales, alpha, src_l, beta, dst_v)
Accumulate potential, induced by local spherical harmonics.
 
subroutine fmm_m2p_mat(c, r, p, vscales, mat)
Accumulate potentials, induced by each multipole spherical harmonic.
 
subroutine fmm_m2l_rotation_work(c, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
Direct M2L translation by 4 rotations and 1 translation.
 
subroutine fmm_l2l_ztranslate_adj_work(z, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Adjoint L2L translation over OZ axis.
 
subroutine fmm_l2l_bessel_rotation_adj(c, src_r, dst_r, kappa, p, vscales, alpha, src_l, beta, dst_l)
Direct L2L translation by 4 rotations and 1 translation.
 
subroutine fmm_l2l_rotation_adj_work(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Adjoint L2L translation by 4 rotations and 1 translation.
 
subroutine fmm_m2m_bessel_ztranslate_adj(z, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation over OZ axis.
 
subroutine fmm_l2l_bessel_ztranslate(z, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2l_bessel_ztranslate_adj_work(z, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2m_bessel_ztranslate_adj_work(z, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2m_ztranslate_adj_work(z, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Adjoint M2M translation over OZ axis.
 
subroutine fmm_l2l_ztranslate(z, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l)
Direct L2L translation over OZ axis.
 
subroutine fmm_sph_transform_work(p, r1, alpha, src, beta, dst, work)
Transform spherical harmonics to a new cartesion system of coordinates.
 
subroutine fmm_m2l_bessel_rotation_work(c, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_m2l_bessel_rotation_adj(c, src_r, dst_r, kappa, p, vscales, alpha, src_m, beta, dst_l)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_l2l_bessel_rotation_work(c, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct L2L translation by 4 rotations and 1 translation.
 
subroutine fmm_m2l_rotation_adj_work(c, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
Adjoint M2L translation by 4 rotations and 1 translation.
 
subroutine fmm_m2p(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v)
Accumulate potential, induced by multipole spherical harmonics.
 
subroutine fmm_m2l_bessel_rotation_adj_work(c, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_m2m_scale_adj(src_r, dst_r, p, alpha, src_m, beta, dst_m)
Adjoint scale M2M, when spherical harmonics are centered in the same point.
 
subroutine fmm_m2m_bessel_ztranslate_work(z, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2l_bessel_ztranslate_work(z, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2l_bessel_ztranslate_adj(z, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l)
Direct M2M translation over OZ axis.
 
subroutine fmm_m2p_bessel(c, src_r, p, vscales, alpha, src_m, beta, dst_v)
Accumulate potential, induced by multipole spherical harmonics.
 
subroutine fmm_constants(dmax, pm, pl, vcnk, m2l_ztranslate_coef, m2l_ztranslate_adj_coef)
Compute FMM-related constants.
 
subroutine fmm_p2m_work(c, src_q, dst_r, p, vscales, beta, dst_m, work)
Accumulate a multipole expansion induced by a particle of a given charge.
 
subroutine fmm_sph_rotate_oz_work(p, vcos, vsin, alpha, src, beta, dst)
Rotate spherical harmonics around OZ axis.
 
subroutine fmm_l2p_bessel_work(c, p, vscales, SI_ri, alpha, src_l, beta, dst_v, work_complex, work)
Accumulate Bessel potential, induced by local spherical harmonics.
 
subroutine fmm_m2m_bessel_rotation(c, src_r, dst_r, kappa, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_m2l_ztranslate_adj(z, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
Adjoint M2L translation over OZ axis.
 
subroutine fmm_m2m_rotation_adj_work(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Adjoint M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_m2p_bessel_work(c, p, vscales, SK_ri, alpha, src_m, beta, dst_v, work_complex, work)
Accumulate Bessel potential, induced by multipole spherical harmonics.
 
subroutine fmm_m2m_bessel_rotation_adj(c, src_r, dst_r, kappa, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation by 4 rotations and 1 translation.
 
subroutine fmm_l2p_work(c, src_r, p, vscales_rel, alpha, src_l, beta, dst_v, work)
Accumulate potential, induced by local spherical harmonics.
 
subroutine fmm_m2m_rotation_adj(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m)
Adjoint M2M translation by 4 rotations and 1 translation.