subroutine rfdpar1(be,rate,m)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rfdpar
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract:  This routine finds the quadratic factors and (if present) the
!            single linear factor of the generic quasi-Gaussian recursive
!            filter of degree m. 
!            
!
! program history log:
!   1998-01-01  purser
!   2004-06-22  treadon - update documentation
!   2004-06-23  purser  - update documentation
!   2011-07-03  todling - use general interface to FORTRAN intrisic math
!
!   input argument list:
!     m    - the degree of the recursive filter
!
!   output argument list:
!     be   - Recursive filter "beta" coefficients for characteristic modes
!     rate - Decay rate factors for the characteristic modes of the filter
!
!   remarks:
!            This is one of several possible variants of 
!            the construction procedure described in Purser et al, 2003, MWR
!            131, 1524--1535. An earlier version, with additional material,
!            is Purser et al. 2001: NCEP Office Note 431. The style of
!            filtering implemented here treats the (in general, complex)
!            characteristic eigenmodes of the idealized degree-m recursive
!            filter additively. This style was motivated by evidence that
!            this approach led to better control of amplitude in the case
!            of spatial variations of filtering or grid scale. However, the
!            exact construction and application of the filter parameters in
!            this style differs from the published descriptions (referenced
!            above) which do not attempt to analyze and separate the spatial
!            eigenmodes of the filter as is done here.
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use constants, only: zero,half,one
  implicit none

  integer(i_kind)          ,intent(in   ) :: m
  real(r_kind),dimension(*),intent(  out) :: be,rate

  integer(i_kind),parameter:: nn=12
  real(r_kind),parameter:: qcrit=0.001_r_kind


  logical polish
  integer(i_kind) j2,jreal,ipow,jimag,kmod2,i
  real(r_kind) qa,q,r
  real(r_kind),dimension(nn,nn):: van
  complex(r_kind) ca,cb,cc
  complex(r_kind),dimension(nn):: croot
  complex(r_kind),dimension(0:nn):: cof

  data polish/.true./

  polish=.true.
  kmod2=mod(m,2)
! Set up the coefficients of the polynomial in z of degree m that approximates
! the function, exp(z/2):
  cof=zero
  cof(0)=one
  do i=1,m
     cof(i)=half*cof(i-1)/float(i)
  enddo
! Locate the m roots of this polynomial:
  call zroots(cof,m,croot,polish)

! If m is even, all roots are complex; if odd, there is one real root,
! which is the first in the list, croot, returned by subr. zroots:
  if(kmod2==1)then    ! treat the single real root:
     r=-real(croot(1),r_kind)
     q=-aimag(croot(1))
     qa=abs(q)
     if(qa>qcrit) then
        write(6,*)'RFDPAR1:  imaginary root qa=',qa,' > qcrit=',qcrit
        call stop2(69)
     endif
     r=sqrt(r)
     rate(1)=r
     van(1,1)=one
     van(2,1)=r
     do i=3,m
        van(i,1)=van(i-1,1)*r*r
     enddo
  endif
! The complex roots occur in conjugate pairs and emerge from zroots
! ordered by their real parts. It is therefore sufficient to use stride-2
! when looping through the part of the list, croot, that contains the
! complex roots.
  do j2=2,m,2         ! loop over remaining independent complex roots
     jreal=kmod2+j2-1
     jimag=kmod2+j2
     ca=-croot(j2)
     cb=sqrt(ca)
     r=real(cb,r_kind)
     q=aimag(cb)
     if(r<zero)then
        cb=-cb
        r=-r
        q=-q
     endif
     rate(jreal)=r
     rate(jimag)=q
     van(1,jreal)=one
     van(1,jimag)=zero
     do i=2,m
        ipow=i*2-3
        cc=cb**ipow
        van(i,jreal)=real(cc,r_kind)
        van(i,jimag)=-aimag(cc)
     enddo
  enddo
  be(1)=one
  do i=2,m
     be(i)=zero
  enddo
  call linmm(van,be,m,1,nn,m)
  return
end subroutine rfdpar1


subroutine rfdpar2(be,rate,turn,samp,m)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rfdpar2
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract:  This routine sets the "turning conditions" for a recursive
!            that intercepts a domain boundary.
!
! program history log:
!   1998-01-01  purser
!   2004-06-23  purser - added documentation
!   2011-07-03  todling - use mathkinds for math-intrisic
!
!   input argument list:
!     be   - Recursive filter coefficients listed by characteristic modes
!     rate - Decay rate factors for the characteristic modes of the filter
!     m    - degree of the recursive filter
!
!   output argument list:
!     turn - square matrix of "turning coefficients" for bounded domain
!     samp - amplitude adjustment factor related to normalization.
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: zero,half,two
  implicit none

  integer(i_kind)            ,intent(in   ) :: m
  real(r_kind),dimension(m)  ,intent(in   ) :: be,rate
  real(r_kind),dimension(m,m),intent(  out) :: turn
  real(r_kind)               ,intent(  out) :: samp
  
  integer(i_kind) kr,li,lr,ki,kmod2
  real(r_kind) beki,bekr,be1,r1,s
  complex(r_kind) cbe,crl,cl1,c1l,crk,crj,cbeh,clk,clj
  complex(r_kind) clkr,clki

  kmod2=mod(m,2)
  s=zero
  if(kmod2==1)then     ! the first root is real:
     r1=rate(1)
     be1=be(1)
     turn(1,1)=be1/(two*r1)
     s=s+turn(1,1)*be1
     do lr=kmod2+1,m,2
        li=lr+1
        cbe=cmplx(be(lr),be(li),r_kind)
        crl=cmplx(rate(lr),rate(li),r_kind)
        cl1=cbe/(r1+crl)
        turn(lr,1)=real(cl1,r_kind)
        turn(li,1)=aimag(cl1)
        c1l=be1/(r1+crl)
        turn(1,lr)=real(c1l,r_kind)
        turn(1,li)=-aimag(c1l)
        s=s+turn(lr,1)*be1+turn(1,lr)*be(lr)+turn(1,li)*be(li)
     enddo
  endif
  do kr=kmod2+1,m,2
     ki=kr+1
     crk=cmplx(rate(kr),rate(ki),r_kind)
     crj=conjg(crk)
     bekr=be(kr)
     beki=be(ki)
     do lr=kmod2+1,m,2
        li=lr+1
        cbeh=half*cmplx(be(lr),be(li),r_kind)
        crl=cmplx(rate(lr),rate(li),r_kind)
        clk=cbeh/(crl+crk)
        clj=cbeh/(crl+crj)
        clkr=clk+clj
        clki=clk-clj
        turn(lr,kr)= real(clkr,r_kind)
        turn(li,kr)= aimag(clkr)
        turn(lr,ki)=-aimag(clki)
        turn(li,ki)= real(clki,r_kind)
        s=s+turn(lr,kr)*bekr+turn(lr,ki)*beki
     enddo
  enddo
  samp=half/s
  return
end subroutine rfdpar2


subroutine rfdparv(dsh,rate,al,n,m)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rfdparv
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract: Convert filter mode decay rates and effective grid steps into
!           actual recursive filter "alpha" coefficients for each 
!           characteristic mode of the filter.
!            
! program history log:
!   1998-01-01  purser
!   2004-06-23  purser - added documentation
!   2004-08-18  derber - remove complex variables, add intent
!
!   input argument list:
!     dsh  - Effective grid step in covariance-scale units
!     rate - Decay rate factors for the characteristic modes of the filter
!     n    - size of the grid in the active direction
!     m    - degree of the recursive filter
!
!   output argument list:
!     al   - "alpha" coefficients for each of the characteristic modes 
!
!$$$
  use kinds, only: r_kind,i_kind
  implicit none

  integer(i_kind)            ,intent(in   ) :: n,m
  real(r_kind),dimension(n)  ,intent(in   ) :: dsh
  real(r_kind),dimension(m)  ,intent(in   ) :: rate
  real(r_kind),dimension(n,m),intent(  out) :: al

  integer(i_kind) i,kr,ki,kmod2
  real(r_kind):: c0i,c0r,c1,c2
  
  kmod2=mod(m,2)
  if(kmod2 == 1)then
     c0r=-rate(1)
     do i=1,n
        al(i,1)=exp(c0r*dsh(i))
     enddo
  endif
  do kr=kmod2+1,m,2
     ki=kr+1
     c0i=-rate(ki)
     c0r=-rate(kr)
     do i=1,n
        c1=c0i*dsh(i)
        c2=exp(c0r*dsh(i))
        al(i,kr)=c2*cos(c1)
        al(i,ki)=c2*sin(c1)
     end do
  enddo
  return
end subroutine rfdparv


subroutine linmm(a,b,m,mm,na,nb)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    linmm
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract:  This routine carries out the linear inversion for a system of
!            m equations with mm right-hand side vectors.
!
! program history log:
!   1998-01-01  purser
!   2004-06-23  purser  - added documentation
!   2013-01-26  parrish - WCOSS debug compile error -- change input a from intent(in) to intent(inout)
!
!   input argument list:
!     a    - the system matrix (not preserved on output)
!     b    - the right-hand side vectors (output as solution vectors)
!     m    - the order of the system matrix
!     mm   - the number of right-hand side vectors.
!     na   - first fortran dimension of array, a
!     nb   - forst fortran dimension of array, b
!
!   output argument list:
!     b    - the mm solution vectors
!
!$$$
  use kinds, only: r_kind,i_kind
  implicit none

  integer(i_kind)             ,intent(in   ) :: m,mm,na,nb
  real(r_kind),dimension(na,*),intent(inout) :: a
  real(r_kind),dimension(nb,*),intent(inout) :: b

  integer(i_kind),dimension(m):: ipiv    ! <- numerical pivot sequence
  real(r_kind) d

  call ldum(a,ipiv,d,m,na) ! <- perform an L-D-U decomposition with pivoting
  call udlmm(a,b,ipiv,m,mm,na,nb) ! <- perform back-substitutions

  return
end subroutine linmm


subroutine ldum(a,ipiv,d,m,na)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ldum
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract:  This routine performs an L-D-U decomposition (in place) with
!            a partial pivoting sequence
!
! program history log:
!   1998-01-01  purser
!   2004-06-23  purser  - added documentation
!
!   input argument list:
!     a    - matrix of a linear system 
!     m    - order of the system matrix
!     na   - first fortran dimension of matrix, a
!
!   output argument list:
!     a    - pivoted L-D-U decomposed system matrix
!     ipiv - pivoting sequence code (list of successive row-swap targets)
!     d    - (real) determinant sign-change flag, +1 or -1.
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: zero,one
  implicit none

  integer(i_kind)                ,intent(in   ) :: m,na
  integer(i_kind),dimension(m)   ,intent(  out) :: ipiv
  real(r_kind)   ,dimension(na,*),intent(inout) :: a
  real(r_kind)                   ,intent(  out) :: d

  integer(i_kind),parameter:: nn=500

  integer(i_kind) k,ibig,jm,i,jp,j
  real(r_kind) aa,aam,t,abig,ajji,aij,ajj
  real(r_kind),dimension(nn):: s

  if(m>nn) then
     write(6,*)'LDUM:  matrix order m=',m,' greater than assumed limit nn=',nn
     call stop2(67)
  endif
  do i=1,m
     aam=zero
     do j=1,m
        aa=abs(a(i,j))
        if(aa>aam) then
           aam=aa
        end if
     enddo
     if (aam==zero) then
        write(6,*)'LDUM:  row ',i,' of matrix in ldum vanishes'
        call stop2(65)
     endif
     s(i)=one/aam
  enddo
  d=one
  ipiv(m)=m ! <- set default to "no swap"
  do j=1,m-1
     jp=j+1
     abig=s(j)*abs(a(j,j))
     ibig=j
     do i=jp,m
        aa=s(i)*abs(a(i,j))
        if(aa>abig)then
           ibig=i
           abig=aa
        endif
     enddo
!  swap rows, recording changed sign of determinant
     ipiv(j)=ibig
     if(ibig/=j)then
        d=-d ! <- change sign of determinant each time a row-swap occurs
        do k=1,m
           t=a(j,k)
           a(j,k)=a(ibig,k)
           a(ibig,k)=t
        enddo
        s(ibig)=s(j)
     endif
     ajj=a(j,j)
     if(ajj==zero)then
        jm=j-1
        write(6,*)'LDUM:  ***failure*** matrix singular, rank=',jm
        call stop2(66)
     endif
     ajji=one/ajj
     do i=jp,m
        aij=ajji*a(i,j)
        a(i,j)=aij
        do k=jp,m
           a(i,k)=a(i,k)-aij*a(j,k)
        enddo
     enddo
  enddo
  return
end subroutine ldum


subroutine udlmm(a,b,ipiv,m,mm,na,nb)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    udlmm
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract:  This routine performs the L, the D and the U back-substitution
!            steps on each of the mm right-hand side vectors, taking due 
!            account of the encoded pivoting sequence. Note: the calls to
!            subr. dsbvr could be replaced by the f90 intrinsic, "dot_product",
!            (or "sum"), not available when this routine was originally coded. 
!
! program history log:
!   1998-01-01  purser
!   2004-06-23  purser  - added documentation
!
!   input argument list:
!     a    - L-D-U pivoted factors of matrix, stored in place
!     b    - mm right-hand sides
!     ipiv - pivoting sequence of transposition image indexes
!     m    - order of matrix
!     mm   - number of right-hand sides
!     na   - first fortran dimension of a
!     nb   - first fortran dimension of b
!
!   output argument list:
!     b    - mm solution vectors
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: one
  implicit none

  integer(i_kind)                ,intent(in   ) :: na,nb,m,mm
  integer(i_kind),dimension(m)   ,intent(in   ) :: ipiv
  real(r_kind)   ,dimension(na,*),intent(in   ) :: a
  real(r_kind)   ,dimension(nb,*),intent(inout) :: b

  integer(i_kind) k,i,l
  real(r_kind) s,aiii
  do k=1,mm !loop over columns of b
     do i=1,m
        l=ipiv(i)
        s=b(l,k)
        b(l,k)=b(i,k)
        call dsbvr(b(1,k),a(i,1),s,i-1,na)
        b(i,k)=s
     enddo
     b(m,k)=b(m,k)/a(m,m)
     do i=m-1,1,-1
        aiii=one/a(i,i)
        call dsbvr(b(i+1,k),a(i,i+1),b(i,k),m-i,na)
        b(i,k)=b(i,k)*aiii
     enddo
  enddo
  return
end subroutine udlmm


subroutine dsbvr(d,a,s,m,na)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    dsbvr
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract:  Subtract the dot product of vector d and matrix row a(1,:)
!            from scalar s. (Note: this is a carry-over from the pre-f90 era;
!            a more succinct version of the same task can now be accomplished
!            by using the dot_product intrinsic in the routine that calls this
!            one.)
!
! program history log:
!   1998-01-01  purser
!   2009-08-10  lueken updated subprogram doc block
!
!   input argument list:
!    na,m
!    s
!    d
!    a
!
!   output argument list:
!    s
!
! attributes:
!   language: f90
!   machine:
!
!$$$
  use kinds, only: r_kind,i_kind
  implicit none

  integer(i_kind)             ,intent(in   ) :: na,m
  real(r_kind)                ,intent(inout) :: s
  real(r_kind),dimension(m)   ,intent(in   ) :: d
  real(r_kind),dimension(na,*),intent(in   ) :: a

  integer(i_kind) i

  do i=1,m
     s=s-d(i)*a(1,i)
  enddo
  return
end subroutine dsbvr


subroutine zroots(a,m,roots,polish)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    zroots
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract:  This routine is the root-finding procedure adapted from
!            Press et al. 1992: Numerical Recipes for fortran 77, Cambridge.
!
! program history log:
!   1998-01-01  purser
!   2004-06-23  purser  - added documentation
!   2004-10-28  treadon - replace "tiny" with "tiny_r_kind"
!   2005-03-29  treadon - define small = sqrt(tiny_r_kind)
!   2011-07-03  todling - use mathkinds for math-intrisic
!
!   input argument list:
!     a    - complex array of m+1 polynomial coefficients
!     m    - degree of polynomial
!     polish - logical flag to determine whether to "polish" the root values
!
!   output argument list:
!     roots - array of m complex roots listed in increasing order of real prts
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: zero,two,tiny_r_kind
  implicit none

  logical                     ,intent(in   ) :: polish
  integer(i_kind)             ,intent(in   ) :: m
  complex(r_kind),dimension(*),intent(in   ) :: a
  complex(r_kind),dimension(m),intent(  out) :: roots

  integer(i_kind),parameter:: maxm=101

  integer(i_kind) j,i,jj
  real(r_kind)::  small,twosmall2
  complex(r_kind) x,b,c
  complex(r_kind),dimension(maxm):: ad


  small=sqrt(tiny_r_kind)
  twosmall2=two*small*small

  do j=1,m+1
     ad(j)=a(j)
  end do
  do j=m,1,-1
     x=cmplx(zero,zero,r_kind)
     call laguer(ad,j,x,small,.false.)
     if(abs(aimag(x)) <= twosmall2*abs(real(x,r_kind))) then
        x=cmplx(real(x,r_kind),zero,r_kind)
     end if
     roots(j)=x
     b=ad(j+1)
     do jj=j,1,-1
        c=ad(jj)
        ad(jj)=b
        b=x*b+c
     end do
  end do
  if (polish) then
     do j=1,m
        call laguer(a,m,roots(j),small,.true.)
     end do
  endif
  do j=2,m
     x=roots(j)
     do i=j-1,1,-1
        if(real(roots(i),r_kind)<=real(x,r_kind))then
           roots(i+1)=x
           cycle
        end if
        roots(i+1)=roots(i)
     end do
     roots(1)=x
  end do
  return
end subroutine zroots


subroutine laguer(a,m,x,small,polish)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    laguer
!     prgmmr:    purser      org:  np20              date:  1998-01-01
!
! abstract:  This routine is an adaptation of the "Numerical Recipes"
!            root-refining "Laguerre's Method" subroutine. See Numerical
!            Recipes in Fortran 77, 2nd Ed., by Press, Teukolsky, Vetterling
!            and Flannery, Cambridge.
!
! program history log:
!   1998-01-01  purser
!   2004-06-23  purser  - update documentation
!   2011-07-02  todling - adjust smalls for single precision
!
!   input argument list:
!     a    - complex array of m+1 polynomial coefficients
!     m    - degree of the polynomial
!     small - a small real +ve number comparable with machine precision
!     polish - logical flag to determine whether to "polish" result to
!              high precision
!
!   output argument list:
!     x    - one of the complex roots of the polynomial
! 
!$$$
  use kinds, only: r_kind,i_kind,r_single,r_double
  use constants, only: zero,two
  implicit none

  logical                     ,intent(in   ) :: polish
  complex(r_kind),dimension(*),intent(in   ) :: a(*)
  complex(r_kind)             ,intent(inout) :: x
  integer(i_kind)             ,intent(in   ) :: m
  real(r_kind)                ,intent(in   ) :: small

  integer(i_kind),parameter:: maxit=100

  integer(i_kind) iter,j
  real(r_kind) abx,cdx,err,dxold
  real(r_kind) smalls
  complex(r_kind) dx,x1,b,d,f,g,h,sq,gp,gm,g2

  if(r_kind==r_single) then
     smalls=6.8e-4_r_kind
  else if(r_kind==r_double) then
     smalls=6.8e-8_r_kind
  else
     write(6,*)'LAGUER:  unsupported precision r_kind=', r_kind
     call stop2(99)
  endif

  dxold=abs(x)
  do iter=1,maxit
     b=a(m+1)
     err=abs(b)
     d=zero
     f=zero
     abx=abs(x)
     do j=m,1,-1
        f=x*f+d
        d=x*d+b
        b=x*b+a(j)
        err=abs(b)+abx*err
     end do
     err=smalls*err
     if(abs(b)<=err) then
        return
     else
        g=d/b
        g2=g*g
        h=g2-two*f/b
        sq=sqrt((m-1)*(m*h-g2))
        gp=g+sq
        gm=g-sq
        if(abs(gp)<abs(gm)) then
           gp=gm
        end if
        dx=m/gp
     endif
     x1=x-dx
     if(x==x1)return
     x=x1
     cdx=abs(dx)
     dxold=cdx
     if(.not.polish)then
        if(cdx<=small*abs(x))return
     endif
  end do
  write(6,*)'LAGUER:  ***warning*** too many iterations'
  return
end subroutine laguer