subroutine smoothrf(work,nlevs)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    smoothrf    perform horizontal part of background error
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: smoothrf perform horizontal part of background error
!
! program history log:
!   2000-03-15  wu
!   2004-05-06  derber - combine regional, add multiple layers
!   2004-08-27  kleist - new berror variable
!   2004-10-26  wu - give smallest RF half weight for regional wind variables
!   2004-11-03  treadon - pass horizontal scale weighting factors through berror
!   2004-11-22  derber - add openMP
!   2005-03-09  wgu/kleist - square hzscl in totwgt calculation
!   2005-05-27  kleist/parrish - add option to use new patch interpolation
!                  if (norsp==0) will default to polar cascade
!   2005-11-16  wgu - set nmix=nr+1+(ny-nlat)/2 to make sure
!                  nmix+nrmxb=nr no matter what number nlat is.   
!   2010-05-05  derber create diag2tr - diag2nh -diag2sh routines to simplify smoothrf routines
!   2010-05-22  todling - remove implicit ordering requirement in nvar_id
!   2010-10-08  derber - optimize and clean up
!
!   input argument list:
!     work     - horizontal fields to be smoothed
!     nlevs    - number of vertical levels for smoothing
!
!   output argument list:
!     work     - smoothed horizontal field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon,regional
  use constants, only: zero,half
  use berror, only: ii,jj,ii1,jj1,ii2,jj2,slw,slw1,slw2, &
       nx,ny,nf,hzscl,hswgt,nfg,nhscrf
  use control_vectors, only:  nrf_var
  use mpimod, only:  nvar_id
  use smooth_polcarf, only: smooth_polcas,smooth_polcasa
  implicit none

! Declare passed variables
  integer(i_kind)                        ,intent(in   ) :: nlevs
  real(r_kind),dimension(nlat,nlon,nlevs),intent(inout) :: work

! Declare local variables
  integer(i_kind) j,i
  integer(i_kind) k,kk,kkk

  real(r_kind),dimension(nhscrf):: totwgt
  real(r_kind),allocatable,dimension(:,:) :: pall
  real(r_kind),dimension(nlat,nlon,3*nlevs) :: workout


! Regional case
  if(regional)then
!$omp parallel do  schedule(dynamic,1) private(k,j,totwgt), &
!$omp  shared(nlevs,hswgt,hzscl,nrf_var,nvar_id,work,nx,ny,ii,jj,slw,nhscrf),default(none)
     do k=1,nlevs

!       apply horizontal recursive filters
        do j=1,nhscrf
           totwgt(j)=hswgt(j)*hzscl(j)*hzscl(j)
        end do
        
        if(nrf_var(nvar_id(k))=='sf'.or.nrf_var(nvar_id(k))=='vp')then
           totwgt(3)=half*totwgt(3)
        end if
        
        call rfxyyx(work(1,1,k),ny,nx,ii(1,1,1,k),&
             jj(1,1,1,k),slw(1,k),totwgt)
        
     end do

! Global case
  else

     do j=1,nhscrf
        totwgt(j)=hswgt(j)*hzscl(j)*hzscl(j)
     end do
     
!$omp parallel do  schedule(dynamic,1) private(kk) &
!$omp private(i,j,k,kkk,pall), &
!$omp shared(nlevs,workout,nx,ny,ii,jj,slw,nhscrf,nf,nfg,totwgt, &
!$omp        ii2,slw1,slw2,jj2,jj1,nlat,nlon,ii1,work),default(none)
     do kk=1,3*nlevs

        k=(kk-1)/3+1
        kkk=mod(kk-1,3)+1

        do j=1,nlon
          do i=1,nlat
            workout(i,j,kk)=zero
          end do
        end do

!       Recursive filter applications

        if(kkk == 1)then

!         equatorial/mid-latitude band
          allocate(pall(ny,nx))
          call grid2tr(work(1,1,k),pall)
          call rfxyyx(pall,ny,nx,ii(1,1,1,k),jj(1,1,1,k),slw(1,k),totwgt)
          call grid2tr_ad(workout(1,1,kk),pall)
          deallocate(pall)

        else if(kkk == 2)then

!         North pole patch --interpolate - recursive filter - adjoint interpolate
          allocate(pall(-nf:nf,-nf:nf))
          call grid2nh(work(1,1,k),pall)
          call rfxyyx(pall,nfg,nfg,ii1(1,1,1,k),jj1(1,1,1,k),slw1(1,k),totwgt)
          call grid2nh_ad(workout(1,1,kk),pall)
          deallocate(pall)
 
        else if (kkk == 3)then

!         South pole patch --interpolate - recursive filter - adjoint interpolate
          allocate(pall(-nf:nf,-nf:nf))
          call grid2sh(work(1,1,k),pall)
          call rfxyyx(pall,nfg,nfg,ii2(1,1,1,k),jj2(1,1,1,k),slw2(1,k),totwgt)
          call grid2sh_ad(workout(1,1,kk),pall)
          deallocate(pall)

        end if
        
!    End of kk loop over 3*nlevs
     end do

!    Sum up three different patches  for each level
     do kk=1,nlevs
       kkk=(kk-1)*3
       do j=1,nlon
         do i=1,nlat
           work(i,j,kk)=workout(i,j,kkk+1) + workout(i,j,kkk+2) + &
                        workout(i,j,kkk+3)
         end do
       end do
     end do

! End of global block
  end if

  return
end subroutine smoothrf

subroutine grid2tr(work,p1all)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grid2tr    perform transformation from grid to tropics patch
!   prgmmr: derber           org: np2                date: 2010-04-29
!
! abstract: grid2tr perform transformation from grid to tropics patch
!
! program history log:
!   2010-04-29  derber
!   input argument list:
!     work     - horizontal field to be transformed
!
!   output argument list:
!     p1all    - output tropics field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon
  use constants, only: zero
  use berror, only: bl,bl2,nx,ny,ndy,ndx,nmix,ndx2,nymx
  implicit none

! Declare passed variables
  real(r_kind),dimension(nlat,nlon),intent(in)  :: work
  real(r_kind),dimension(ny,nx),intent(out)     :: p1all

! Declare local variables
  integer(i_kind) j,i,i2,i1,j1

! -----------------------------------------------------------------------------
  do j=1,nx
     do i=1,ny
        p1all(i,j)=zero
     end do
  end do
! Extract central patch (band) from full grid (work --> p1)
! Blending zones
  do i=1,ndx
     i1=i-ndx+nlon
     i2=nx-ndx+i
     do j=1,ny
        j1=j+ndy
        p1all(j,i) =work(j1,i1)      ! left (west) blending zone
        p1all(j,i2)=work(j1,i)       ! right (east) blending zone
     enddo
  enddo

! Middle zone (no blending)
  do i=ndx+1,nx-ndx
     i1=i-ndx
     do j=1,ny
        p1all(j,i)=work(j+ndy,i1)
     enddo
  enddo

! Apply blending coefficients to central patch
  do i=1,ndx2
     i1=ndx2+1-i
     i2=nx-ndx2+i
     do j=1,ny
        p1all(j,i) =p1all(j,i) *bl(i1)  ! left (west) blending zone
        p1all(j,i2)=p1all(j,i2)*bl(i)   ! right (east) blending zone
     enddo
  enddo

! bl2 of p1
  do i=1,nx
     do j=1,nmix
        p1all(j,i)=p1all(j,i)*bl2(nmix+1-j)
     enddo
     do j=nymx+1,ny
        p1all(j,i)=p1all(j,i)*bl2(j-nymx)
     enddo
  enddo


  return
end subroutine grid2tr

subroutine grid2tr_ad(work,p1all)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grid2tr    perform adjoint of transformation from grid to tropics patch
!   prgmmr: derber           org: np2                date: 2010-04-29
!
! abstract: grid2tr perform adjoint of transformation from grid to tropics patch
!
! program history log:
!   2010-04-29  derber
!   input argument list:
!     p1all    - input nh field
!
!   output argument list:
!     work     - horizontal field to be transformed
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon
  use constants, only: zero
  use berror, only: bl,bl2,nx,ny,ndx,ndy,ndx2,nmix,nymx
  implicit none

! Declare passed variables
  real(r_kind),dimension(nlat,nlon),intent(inout)  :: work
  real(r_kind),dimension(ny,nx),intent(inout)      :: p1all

! Declare local variables
  integer(i_kind) j,i,i2,i1,j1

! -----------------------------------------------------------------------------

! Equatorial patch
! Adjoint of central patch blending on left/right sides of patch
  do i=1,ndx2
     i1=ndx2+1-i
     i2=nx-ndx2+i
     do j=1,ny
        p1all(j,i) =p1all(j,i) *bl(i1)   ! left (west) blending zone
        p1all(j,i2)=p1all(j,i2)*bl(i)    ! right (east) blending zone
     enddo
  enddo

! bl2 of p1
  do i=1,nx
     do j=1,nmix
        p1all(j,i)=p1all(j,i)*bl2(nmix+1-j)
     enddo
     do j=nymx+1,ny
        p1all(j,i)=p1all(j,i)*bl2(j-nymx)
     enddo
  enddo

!   Adjoint of transfer between central band and full grid (p1 --> work)
  do i=1,ndx
     i1=i-ndx+nlon
     i2=nx-ndx+i
     do j=1,ny
        j1=j+ndy
        work(j1,i1)=work(j1,i1)+p1all(j,i)  ! left (west) blending zone
        work(j1,i) =work(j1,i) +p1all(j,i2) ! right (east) blending zone
     enddo
  enddo

! Middle zone (no blending)
  do i=ndx+1,nx-ndx
     i1=i-ndx
     do j=1,ny
        j1=j+ndy
        work(j1,i1)=work(j1,i1)+p1all(j,i)
     enddo
  enddo

  return
end subroutine grid2tr_ad
subroutine grid2nh(work,pall)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grid2nh    perform transformation from grid to nh patch
!   prgmmr: derber           org: np2                date: 2010-04-29
!
! abstract: grid2nh perform transformation from grid to nh patch
!
! program history log:
!   2010-04-29  derber
!   input argument list:
!     work     - horizontal field to be transformed
!
!   output argument list:
!     pall     - output nh field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon
  use constants, only: zero
  use berror, only: wtaxs,wtxrs,inaxs,inxrs,bl2,mr,nr,nf,ndy,norm,nxem
  use smooth_polcarf, only: norsp,smooth_polcas,smooth_polcasa
  implicit none

! Declare passed variables
  real(r_kind),dimension(nlat,nlon),intent(in)  :: work
  real(r_kind),dimension(-nf:nf,-nf:nf),intent(inout)    :: pall

! Declare local variables
  real(r_kind),dimension(nlon+1,mr:nr)    :: p2all
  integer(i_kind) j,i,j1

! -----------------------------------------------------------------------------
  do j=mr,nr
     do i=1,nlon+1
        p2all(i,j)=zero
     end do
  end do
! North pole patch(p2) -- blending and transfer to grid

  do i=1,nlon
!    Load field into patches
     do j=mr,nr
        p2all(i,j)=work(nlat-j,i)
     enddo
  enddo
! Apply blending coefficients
  do j=ndy,nr
     j1=j-ndy+1
     do i=1,nlon
        p2all(i,j)=p2all(i,j)*bl2(j1)
     enddo
  enddo
! Interpolation to polar grid
  if(norsp>0) then
     call smooth_polcasa(pall,p2all)
  else
     call polcasa(pall,p2all,nxem,norm,nlon,wtaxs,wtxrs,inaxs,inxrs,nf,mr,nr)
  end if

  return
end subroutine grid2nh

subroutine grid2nh_ad(work,pall)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grid2nh    perform adjoint of transformation from grid to nh patch
!   prgmmr: derber           org: np2                date: 2010-04-29
!
! abstract: grid2nh perform adjoint of transformation from grid to nh patch
!
! program history log:
!   2010-04-29  derber
!   input argument list:
!     pall     - input nh field
!
!   output argument list:
!     work     - horizontal field to be transformed
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon
  use constants, only: zero
  use berror, only: wtaxs,wtxrs,inaxs,inxrs,bl2,mr,nr,nf,ndy,norm,nxem
  use smooth_polcarf, only: norsp,smooth_polcas,smooth_polcasa
  implicit none

! Declare passed variables
  real(r_kind),dimension(nlat,nlon),intent(inout)  :: work
  real(r_kind),dimension(-nf:nf,-nf:nf),intent(inout)    :: pall

! Declare local variables
  real(r_kind),dimension(nlon+1,mr:nr)     :: p2all
  integer(i_kind) j,i,j1

! -----------------------------------------------------------------------------

! Adjoint of interpolation to polar grid
  if(norsp>0) then
     call smooth_polcas(pall,p2all)
  else
     call polcas(pall,p2all,nxem,norm,nlon,wtaxs,wtxrs,inaxs,inxrs,nf,mr,nr)
  end if
! North pole patch(p2) -- adjoint of blending and transfer to grid

! Apply blending coefficients
  do j=ndy,nr
     j1=j-ndy+1
     do i=1,nlon
        p2all(i,j)=p2all(i,j)*bl2(j1)
     enddo
  enddo

  do i=1,nlon
!    Load field into patches
     do j=mr,nr
        work(nlat-j,i)=work(nlat-j,i)+p2all(i,j)
     enddo
  enddo

  return
end subroutine grid2nh_ad

subroutine grid2sh(work,pall)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grid2sh    perform transformation from grid to sh patch
!   prgmmr: derber           org: np2                date: 2010-04-29
!
! abstract: grid2sh perform transformation from grid to sh patch
!
! program history log:
!   2010-04-29  derber
!   input argument list:
!     work     - horizontal field to be transformed
!
!   output argument list:
!     pall     - output sh field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon
  use constants, only: zero
  use berror, only: wtaxs,wtxrs,inaxs,inxrs,bl2,mr,nr,nf,ndy,norm,nxem
  use smooth_polcarf, only: norsp,smooth_polcasa
  implicit none

! Declare passed variables
  real(r_kind),dimension(nlat,nlon),intent(in)  :: work
  real(r_kind),dimension(-nf:nf,-nf:nf),intent(inout)    :: pall

! Declare local variables
  real(r_kind),dimension(nlon+1,mr:nr)    :: p3all
  integer(i_kind) j,i,j1

! -----------------------------------------------------------------------------

  do j=mr,nr
     do i=1,nlon+1
        p3all(i,j)=zero
     end do
  end do
! south pole patch(p3) -- blending and transfer to grid

  do i=1,nlon
!    Load field into patches
     do j=mr,nr
        p3all(i,j)=work(j+1,i)
     enddo
  enddo
! Apply blending coefficients
  do j=ndy,nr
     j1=j-ndy+1
     do i=1,nlon
        p3all(i,j)=p3all(i,j)*bl2(j1)
     enddo
  enddo
! Interpolate to polar grid
  if(norsp>0) then
     call smooth_polcasa(pall,p3all)
  else
     call polcasa(pall,p3all,nxem,norm,nlon,wtaxs,wtxrs,inaxs,inxrs,nf,mr,nr)
  end if

  return
end subroutine grid2sh

subroutine grid2sh_ad(work,pall)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grid2sh    perform adjoint of transformation from grid to sh patch
!   prgmmr: derber           org: np2                date: 2010-04-29
!
! abstract: grid2sh perform adjoint of transformation from grid to sh patch
!
! program history log:
!   2010-04-29  derber
!   input argument list:
!     pall     - input sh field
!
!   output argument list:
!     work     - horizontal field to be transformed
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon
  use constants, only: zero
  use berror, only: wtaxs,wtxrs,inaxs,inxrs,bl2,mr,nr,nf,ndy,norm,nxem
  use smooth_polcarf, only: norsp,smooth_polcas
  implicit none

! Declare passed variables
  real(r_kind),dimension(nlat,nlon),intent(inout)        :: work
  real(r_kind),dimension(-nf:nf,-nf:nf),intent(inout)    :: pall

! Declare local variables
  real(r_kind),dimension(nlon+1,mr:nr)     :: p3all
  integer(i_kind) j,i,j1


! Interpolate to polar grid
  if(norsp>0) then
     call smooth_polcas(pall,p3all)
  else
     call polcas(pall,p3all,nxem,norm,nlon,wtaxs,wtxrs,inaxs,inxrs,nf,mr,nr)
  end if

! South pole patch(p2) -- adjoint of blending and transfer to grid

! Apply blending coefficients
  do j=ndy,nr
     j1=j-ndy+1
     do i=1,nlon
        p3all(i,j)=p3all(i,j)*bl2(j1)
     enddo
  enddo

  do i=1,nlon
!    Load field into patches
     do j=mr,nr
        work(j+1,i)=work(j+1,i)+p3all(i,j)
     enddo
  enddo

  return
end subroutine grid2sh_ad


subroutine rfxyyx(p1,nx,ny,iix,jjx,dssx,totwgt)
  
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rfxyyx      perform horizontal smoothing
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: smoothrf perform self-adjoint horizontal smoothing. nsloop
!           smoothing fields.
!
! program history log:
!   2000-03-15  wu
!   2004-08-24  derber - change indexing add rfhyt to speed things up
!
!   input argument list:
!     p1       - horizontal field to be smoothed
!     nx       - first dimension of p1
!     ny       - second dimension of p1
!     iix      - array of pointers for smoothing table (first dimension)
!     jjx      - array of pointers for smoothing table (second dimension)
!     dssx     - renormalization constants including variance
!     wgt      - weight (empirical*expected)
!
!   output argument list:
!                 all after horizontal smoothing
!     p1       - horizontal field which has been smoothed
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only:  zero
  use berror, only: be,table,ndeg,nhscrf
  implicit none

! Declare passed variables
  integer(i_kind)                        ,intent(in   ) :: nx,ny
  integer(i_kind),dimension(nx,ny,nhscrf),intent(in   ) :: iix,jjx
  real(r_kind),dimension(nx,ny)          ,intent(inout) :: p1
  real(r_kind),dimension(nx,ny)          ,intent(in   ) :: dssx
  real(r_kind),dimension(nhscrf)         ,intent(in   ) :: totwgt

! Declare local variables
  integer(i_kind) ix,iy,i,j,im,n

  real(r_kind),dimension(nx,ny):: p2,p1t,p1in
  real(r_kind),dimension(nx,ny,ndeg):: alx,aly

! Zero local arrays
  do iy=1,ny
     do ix=1,nx
        p1in(ix,iy)=p1(ix,iy)
        p1(ix,iy)=zero
     enddo
  enddo

! Loop over number of scales
 
  do n=1,nhscrf

     do iy=1,ny
        do ix=1,nx
           p2(ix,iy)=zero
        enddo
     enddo
     do im=1,ndeg
        do j=1,ny
           do i=1,nx
              alx(i,j,im)=table(iix(i,j,n),im)
              aly(i,j,im)=table(jjx(i,j,n),im)
           enddo
        enddo
     enddo

     call rfhx0(p1in,p2,nx,ny,ndeg,alx,be)

     call rfhyt(p2,p1t,nx,ny,ndeg,aly,be)

     do iy=1,ny
        do ix=1,nx
           p1t(ix,iy)=p1t(ix,iy)*dssx(ix,iy)*totwgt(n)
        enddo
     enddo

     call rfhy(p1t,p2,nx,ny,ndeg,aly,be)

     call rfhx0(p2,p1,nx,ny,ndeg,alx,be)

! end loop over number of horizontal scales
  end do

  return
end subroutine rfxyyx
! -----------------------------------------------------------------------------
subroutine rfhx0(p1,p2,nx,ny,ndeg,alx,be)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:   rfhx0        performs x component of recursive filter
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: performs x component of recursive filter
!
! program history log:
!   2000-03-15  wu
!   2004-05-06  derber  combine regional, add multiple layers
!   2004-08-24  derber change indexing to 1-nx,1-ny
!
!   input argument list:
!     p1       - field to be smoothed
!     nx       - first dimension of p1
!     ny       - second dimension of p1
!     ndeg     - degree of smoothing   
!     alx      - smoothing coefficients
!     be       - smoothing coefficients
!
!   output argument list:
!     p2       - field after smoothing
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only:  zero
  implicit none

  integer(i_kind)                   ,intent(in   ) :: nx,ny,ndeg
  real(r_kind),dimension(nx,ny)     ,intent(in   ) :: p1
  real(r_kind),dimension(ndeg)      ,intent(in   ) :: be
  real(r_kind),dimension(nx,ny,ndeg),intent(in   ) :: alx
  real(r_kind),dimension(nx,ny)     ,intent(inout) :: p2

  integer(i_kind) ix,iy,kr,ki

  real(r_kind) gakr,gaki,dekr,deki,bekr,beki
  real(r_kind),dimension(ndeg,ny)   :: gap,dep

  do iy=1,ny
    do ix=1,ndeg
      gap(ix,iy)=zero
      dep(ix,iy)=zero
    end do
  end do

  if (mod(ndeg,2) == 1) then  

!    Advancing filter:
     do ix=1,nx
        do iy=1,ny
           gap(1,iy)=alx(ix,iy,1)*gap(1,iy)+be(1)*p1(ix,iy)
           p2(ix,iy)=p2(ix,iy)+gap(1,iy)
        enddo

                           ! treat remaining complex roots:
        do kr=2,ndeg,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do iy=1,ny
              gakr=gap(kr,iy)
              gaki=gap(ki,iy)
              gap(kr,iy)=alx(ix,iy,kr)*gakr-alx(ix,iy,ki)*gaki+bekr*p1(ix,iy)
              gap(ki,iy)=alx(ix,iy,ki)*gakr+alx(ix,iy,kr)*gaki+beki*p1(ix,iy)
              p2(ix,iy)=p2(ix,iy)+gap(kr,iy)
           enddo
        enddo
     enddo

! Backing filter:
     do ix=nx,1,-1
!       treat real roots
        do iy=1,ny
           p2(ix,iy)=p2(ix,iy)+dep(1,iy)
           dep(1,iy)=alx(ix,iy,1)*(dep(1,iy)+be(1)*p1(ix,iy))
        enddo
                           ! treat remaining complex roots:
        do kr=2,ndeg,2   ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           do iy=1,ny
              p2(ix,iy)=p2(ix,iy)+dep(kr,iy)
              dekr=dep(kr,iy)+bekr*p1(ix,iy)
              deki=dep(ki,iy)+beki*p1(ix,iy)
              dep(kr,iy)=alx(ix,iy,kr)*dekr-alx(ix,iy,ki)*deki
              dep(ki,iy)=alx(ix,iy,ki)*dekr+alx(ix,iy,kr)*deki
           enddo
        enddo
     enddo

  else
     do iy=1,ny

        !       Advancing filter
        ! treat remaining complex roots:
        do kr=1,ndeg,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              gakr=gap(kr,iy)
              gaki=gap(ki,iy)
              gap(kr,iy)=alx(ix,iy,kr)*gakr-alx(ix,iy,ki)*gaki+bekr*p1(ix,iy)
              gap(ki,iy)=alx(ix,iy,ki)*gakr+alx(ix,iy,kr)*gaki+beki*p1(ix,iy)
              p2(ix,iy)=p2(ix,iy)+gap(kr,iy)
              
           end do
           
        !       Backing filter:
        ! treat remaining complex roots:
           do ix=nx,1,-1
              p2(ix,iy)=p2(ix,iy)+dep(kr,iy)
              dekr=dep(kr,iy)+bekr*p1(ix,iy)
              deki=dep(ki,iy)+beki*p1(ix,iy)
              dep(kr,iy)=alx(ix,iy,kr)*dekr-alx(ix,iy,ki)*deki
              dep(ki,iy)=alx(ix,iy,ki)*dekr+alx(ix,iy,kr)*deki
              
           enddo
        end do
     end do
  endif
  return
end subroutine rfhx0
! -----------------------------------------------------------------------------
subroutine rfhyt(p1,p2,nx,ny,ndegy,aly,be)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:   rfhyt        performs x component of recursive filter
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: performs x component of recursive filter
!
! program history log:
!   2000-03-15  wu
!   2004-05-06  derber  combine regional, add multiple layers
!   2004-08-24  derber create rfhyt from rfhy - remove unnecessary computations
!                      remove unused parameters - change indexing
!
!   input argument list:
!     p1       - field to be smoothed
!     nx       - first dimension of p1
!     ny       - second dimension of p1
!     ndegy    - degree of smoothing y direction
!     aly      - smoothing coefficients y direction
!     be       - smoothing coefficients
!
!   output argument list:
!     p2       - field after smoothing
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$

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

  integer(i_kind)                    ,intent(in   ) :: nx,ny,ndegy
  real(r_kind),dimension(nx,ny)      ,intent(in   ) :: p1
  real(r_kind),dimension(nx,ny,ndegy),intent(in   ) :: aly
  real(r_kind),dimension(ndegy)      ,intent(in   ) :: be
  real(r_kind),dimension(nx,ny)      ,intent(  out) :: p2

  integer(i_kind) ix,iy,kr,ki,ly

  real(r_kind),dimension(nx,ndegy):: gap,dep
  real(r_kind) gakr,gaki,dekr,deki
  real(r_kind) beki,bekr


  do ly=1,ndegy
     do ix=1,nx
        gap(ix,ly)=zero
        dep(ix,ly)=zero
     enddo
  enddo

  if (mod(ndegy,2) == 1) then

! Advancing filter:
     do iy=1,ny
!       treat the real root:
        do ix=1,nx
           gap(ix,1)=aly(ix,iy,1)*gap(ix,1)+be(1)*p1(ix,iy)
           p2(ix,iy)=gap(ix,1)
        enddo
                           ! treat remaining complex roots:
        do kr=2,ndegy,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              gakr=gap(ix,kr)
              gaki=gap(ix,ki)
              gap(ix,kr)=aly(ix,iy,kr)*gakr-aly(ix,iy,ki)*gaki+bekr*p1(ix,iy)
              gap(ix,ki)=aly(ix,iy,ki)*gakr+aly(ix,iy,kr)*gaki+beki*p1(ix,iy)
              p2(ix,iy)=p2(ix,iy)+gap(ix,kr)
           enddo
        enddo
     enddo

! Backing filter:
     do iy=ny,1,-1
!       treat the real root:
        do ix=1,nx
           p2(ix,iy)=p2(ix,iy)+dep(ix,1)
           dep(ix,1)=aly(ix,iy,1)*(dep(ix,1)+be(1)*p1(ix,iy))
        enddo
                           ! treat remaining complex roots:
        do kr=2,ndegy,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              p2(ix,iy)=p2(ix,iy)+dep(ix,kr)
              dekr=dep(ix,kr)+bekr*p1(ix,iy)
              deki=dep(ix,ki)+beki*p1(ix,iy)
              dep(ix,kr)=aly(ix,iy,kr)*dekr-aly(ix,iy,ki)*deki
              dep(ix,ki)=aly(ix,iy,ki)*dekr+aly(ix,iy,kr)*deki
           enddo
        enddo
     enddo

  else  

!    Advancing filter:
     do iy=1,ny
        do ix=1,nx
           p2(ix,iy)=zero
        enddo
        ! treat remaining complex roots:
        do kr=1,ndegy,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              gakr=gap(ix,kr)
              gaki=gap(ix,ki)
              gap(ix,kr)=aly(ix,iy,kr)*gakr-aly(ix,iy,ki)*gaki+bekr*p1(ix,iy)
              gap(ix,ki)=aly(ix,iy,ki)*gakr+aly(ix,iy,kr)*gaki+beki*p1(ix,iy)
              p2(ix,iy)=p2(ix,iy)+gap(ix,kr)
           enddo
        enddo
     enddo
     
!    Backing filter:
     do iy=ny,1,-1
        ! treat remaining complex roots:
        do kr=1,ndegy,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              p2(ix,iy)=p2(ix,iy)+dep(ix,kr)
              dekr=dep(ix,kr)+bekr*p1(ix,iy)
              deki=dep(ix,ki)+beki*p1(ix,iy)
              dep(ix,kr)=aly(ix,iy,kr)*dekr-aly(ix,iy,ki)*deki
              dep(ix,ki)=aly(ix,iy,ki)*dekr+aly(ix,iy,kr)*deki
           enddo
        enddo
     enddo
  endif
  return
end subroutine rfhyt
! -----------------------------------------------------------------------------
subroutine rfhy(p1,p2,nx,ny,ndegy,aly,be)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:   rfhy         performs x component of recursive filter
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: performs x component of recursive filter
!
! program history log:
!   2000-03-15  wu
!   2004-05-06  derber  combine regional, add multiple layers
!   2004-08-24  derber  remove unused parameters and unnecessary computation
!                       change indexing
!
!   input argument list:
!     p1       - field to be smoothed
!     nx       - first dimension of p1
!     ny       - second dimension of p1
!     ndegy    - degree of smoothing y direction
!     aly      - smoothing coefficients y direction
!     be       - smoothing coefficients
!
!   output argument list:
!     p2       - field after smoothing
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$

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

  integer(i_kind)                    ,intent(in   ) :: nx,ny,ndegy
  real(r_kind),dimension(nx,ny)      ,intent(in   ) :: p1
  real(r_kind),dimension(nx,ny,ndegy),intent(in   ) :: aly
  real(r_kind),dimension(ndegy)      ,intent(in   ) :: be
  real(r_kind),dimension(nx,ny)      ,intent(  out) :: p2

  integer(i_kind) ix,iy,kr,ki,ly

  real(r_kind) gakr,gaki,dekr,deki
  real(r_kind) beki,bekr

  real(r_kind),dimension(nx,ndegy):: gap,dep


  do ly=1,ndegy
     do ix=1,nx
        gap(ix,ly)=zero
        dep(ix,ly)=zero
     enddo
  end do

  if (mod(ndegy,2) == 1) then

! Advancing filter:
     do iy=1,ny
!       treat the real root:
        do ix=1,nx
           gap(ix,1)=aly(ix,iy,1)*gap(ix,1)+be(1)*p1(ix,iy)
           p2(ix,iy)=gap(ix,1)
        enddo
                           ! treat remaining complex roots:
        do kr=2,ndegy,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              gakr=gap(ix,kr)
              gaki=gap(ix,ki)
              gap(ix,kr)=aly(ix,iy,kr)*gakr-aly(ix,iy,ki)*gaki+bekr*p1(ix,iy)
              gap(ix,ki)=aly(ix,iy,ki)*gakr+aly(ix,iy,kr)*gaki+beki*p1(ix,iy)
              p2(ix,iy)=p2(ix,iy)+gap(ix,kr)
           enddo
        enddo
     enddo

! Backing filter:
     do iy=ny,1,-1
!       treat the real root:
        do ix=1,nx
           p2(ix,iy)=p2(ix,iy)+dep(ix,1)
           dep(ix,1)=aly(ix,iy,1)*(dep(ix,1)+be(1)*p1(ix,iy))
        enddo
                           ! treat remaining complex roots:
        do kr=2,ndegy,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              p2(ix,iy)=p2(ix,iy)+dep(ix,kr)
              dekr=dep(ix,kr)+bekr*p1(ix,iy)
              deki=dep(ix,ki)+beki*p1(ix,iy)
              dep(ix,kr)=aly(ix,iy,kr)*dekr-aly(ix,iy,ki)*deki
              dep(ix,ki)=aly(ix,iy,ki)*dekr+aly(ix,iy,kr)*deki
           enddo
        enddo
     enddo

  else  

!    Advancing filter:
     do iy=1,ny
        do ix=1,nx
           p2(ix,iy)=zero
        enddo
        ! treat remaining complex roots:
        do kr=1,ndegy,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              gakr=gap(ix,kr)
              gaki=gap(ix,ki)
              gap(ix,kr)=aly(ix,iy,kr)*gakr-aly(ix,iy,ki)*gaki+bekr*p1(ix,iy)
              gap(ix,ki)=aly(ix,iy,ki)*gakr+aly(ix,iy,kr)*gaki+beki*p1(ix,iy)
              p2(ix,iy)=p2(ix,iy)+gap(ix,kr)
           enddo
        enddo
     enddo
     
!    Backing filter:
     do iy=ny,1,-1
        ! treat remaining complex roots:
        do kr=1,ndegy,2  ! <-- index of "real" components
           ki=kr+1      ! <-- index of "imag" components
           bekr=be(kr)
           beki=be(ki)
           do ix=1,nx
              p2(ix,iy)=p2(ix,iy)+dep(ix,kr)
              dekr=dep(ix,kr)+bekr*p1(ix,iy)
              deki=dep(ix,ki)+beki*p1(ix,iy)
              dep(ix,kr)=aly(ix,iy,kr)*dekr-aly(ix,iy,ki)*deki
              dep(ix,ki)=aly(ix,iy,ki)*dekr+aly(ix,iy,kr)*deki
           enddo
        enddo
     enddo
  endif
  return
end subroutine rfhy
! ------------------------------------------------------------------------------
subroutine sqrt_smoothrf(z,work,nlevs)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    sqrt_smoothrf    perform sqrt horizontal part of background error
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: smoothrf perform horizontal part of background error
!
! program history log:
!   2000-03-15  wu
!   2004-05-06  derber - combine regional, add multiple layers
!   2004-08-27  kleist - new berror variable
!   2004-10-26  wu - give smallest RF half weight for regional wind variables
!   2004-11-03  treadon - pass horizontal scale weighting factors through berror
!   2004-11-22  derber - add openMP
!   2005-03-09  wgu/kleist - square hzscl in totwgt calculation
!   2005-05-27  kleist/parrish - add option to use new patch interpolation
!                  if (norsp==0) will default to polar cascade
!   2005-11-16  wgu - set nmix=nr+1+(ny-nlat)/2 to make sure
!                  nmix+nrmxb=nr no matter what number nlat is.   
!   2010-05-22  todling - remove implicit ordering requirement in nvar_id
!
!   input argument list:
!     work     - horizontal fields to be smoothed
!     nlevs    - number of vertical levels for smoothing
!
!   output argument list:
!     work     - smoothed horizontal field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon,regional,nnnn1o
  use jfunc,only: nval_lenz
  use constants, only: zero,half
  use berror, only: ii,jj,ii1,jj1,nhscrf,&
       ii2,jj2,slw,slw1,slw2,nx,ny,nf,hzscl,hswgt,nfg,nfnf
  use control_vectors, only:  nrf_var
  use mpimod, only:  nvar_id
  use smooth_polcarf, only: smooth_polcas
  implicit none

! Declare passed variables
  integer(i_kind)                        ,intent(in   ) :: nlevs
  real(r_kind),dimension(nval_lenz)      ,intent(in   ) :: z
  real(r_kind),dimension(nlat,nlon,nlevs),intent(inout) :: work

! Declare local variables
  integer(i_kind) j,i
  integer(i_kind) k,iz,kk,kkk

  real(r_kind),dimension(nhscrf):: totwgt
  real(r_kind),allocatable,dimension(:,:)::pall,zloc
  real(r_kind),dimension(nlat,nlon,3*nlevs) :: workout

! Regional case
  if(regional)then
!$omp parallel do  schedule(dynamic,1) private(k,i,j,iz,totwgt,zloc), &
!$omp  shared(nlevs,nhscrf,hswgt,hzscl,nrf_var,nvar_id,nlat,nlon,nnnn1o,slw,nx,ii,jj,z,work,ny),default(none)
     do k=1,nlevs

        allocate(zloc(nlat*nlon,nhscrf))
!       apply horizontal recursive filters
        do j=1,nhscrf
           totwgt(j)=sqrt(hswgt(j)*hzscl(j)*hzscl(j))
        end do
        
        if(nrf_var(nvar_id(k))=='sf'.or.nrf_var(nvar_id(k))=='vp')then
           totwgt(3)=sqrt(half)*totwgt(3)
        end if

        do j=1,nhscrf
           iz=nlat*nlon*(k-1)+nlat*nlon*nnnn1o*(j-1)
           do i=1,nlat*nlon
              zloc(i,j)=z(i+iz)
           end do
        end do
        
        call sqrt_rfxyyx(zloc,work(1,1,k),ny,nx,ii(1,1,1,k),&
             jj(1,1,1,k),slw(1,k),totwgt)
        
        deallocate(zloc)
     end do

! Global case
  else

     do j=1,nhscrf
        totwgt(j)=sqrt(hswgt(j))*hzscl(j)
     end do
     
!       zero output array
     do k=1,nlevs
     end do
!$omp parallel do  schedule(dynamic,1) private(kk) &
!$omp private(i,j,k,iz,kkk,pall,zloc), &
!$omp   shared(nlevs,nlon,nlat,nhscrf,nfnf,nx,ny,nnnn1o,nfg,nf,z,ii,jj,ii1,jj1,ii2,jj2,slw2, &
!$omp          slw,slw1,workout,totwgt),default(none)
     do kk=1,3*nlevs

        k=(kk-1)/3+1
        kkk=mod(kk-1,3)+1

        do i=1,nlon
           do j=1,nlat
              workout(j,i,kk)=zero
           end do
        end do
!       Recursive filter applications

        if(kkk == 1)then

!          First do equatorial/mid-latitude band

           allocate(pall(ny,nx),zloc(ny*nx,nhscrf))
           do j=1,nhscrf
              iz=(ny*nx+2*nfnf)*(k-1+nnnn1o*(j-1))
              do i=1,ny*nx
                 zloc(i,j)=z(i+iz)
              end do
           end do
           call sqrt_rfxyyx(zloc,pall,ny,nx,ii(1,1,1,k),jj(1,1,1,k),slw(1,k),totwgt)
           call grid2tr_ad(workout(1,1,kk),pall)
           deallocate(pall,zloc)

        else if(kkk == 2)then

!          North pole patch --interpolate - recursive filter - adjoint interpolate

           allocate(pall(-nf:nf,-nf:nf),zloc(nfnf,nhscrf))
           do j=1,nhscrf
              iz=(ny*nx+2*nfnf)*(k-1+nnnn1o*(j-1))+ny*nx
              do i=1,nfnf
                 zloc(i,j)=z(i+iz)
              end do
           end do
           call sqrt_rfxyyx(zloc,pall,nfg,nfg,ii1(1,1,1,k),jj1(1,1,1,k),slw1(1,k),totwgt)
           call grid2nh_ad(workout(1,1,kk),pall)
           deallocate(pall,zloc)
        else if(kkk == 3)then

!          South pole patch --interpolate - recursive filter - adjoint interpolate

           allocate(pall(-nf:nf,-nf:nf),zloc(nfnf,nhscrf))
           do j=1,nhscrf
              iz=(ny*nx+2*nfnf)*(k-1+nnnn1o*(j-1))+ny*nx+nfnf
              do i=1,nfnf
                 zloc(i,j)=z(i+iz)
              end do
           end do
           call sqrt_rfxyyx(zloc,pall,nfg,nfg,ii2(1,1,1,k),jj2(1,1,1,k),slw2(1,k),totwgt)
           call grid2sh_ad(workout(1,1,kk),pall)
           deallocate(pall,zloc)
        end if
        
!    End of k loop over nlevs
     end do

!    Sum up three different patches  for each level
     do kk=1,nlevs
       kkk=(kk-1)*3
       do j=1,nlon
         do i=1,nlat
           work(i,j,kk)=workout(i,j,kkk+1) + workout(i,j,kkk+2) + &
                        workout(i,j,kkk+3)
         end do
       end do
     end do

! End of global block
  end if

  return
end subroutine sqrt_smoothrf
! ------------------------------------------------------------------------------
subroutine sqrt_smoothrf_ad(z,work,nlevs)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    smoothrf    perform horizontal part of background error
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: smoothrf perform horizontal part of background error
!
! program history log:
!   2000-03-15  wu
!   2004-05-06  derber - combine regional, add multiple layers
!   2004-08-27  kleist - new berror variable
!   2004-10-26  wu - give smallest RF half weight for regional wind variables
!   2004-11-03  treadon - pass horizontal scale weighting factors through berror
!   2004-11-22  derber - add openMP
!   2005-03-09  wgu/kleist - square hzscl in totwgt calculation
!   2005-05-27  kleist/parrish - add option to use new patch interpolation
!                   if (norsp==0) will default to polar cascade
!   2005-11-16  wgu - set nmix=nr+1+(ny-nlat)/2 to make sure
!                   nmix+nrmxb=nr no matter what number nlat is.   
!   2010-05-22  todling - remove implicit ordering requirement in nvar_id
!   2013-01-26  parrish - change work from intent(in) to intent(inout)
!
!   input argument list:
!     work     - horizontal fields to be smoothed
!     nlevs    - number of vertical levels for smoothing
!
!   output argument list:
!     z        - smoothed horizontal field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: nlat,nlon,nnnn1o,regional
  use jfunc,only: nval_lenz
  use constants, only: zero,half
  use berror, only: ii,jj,ii1,jj1,nhscrf,&
       ii2,jj2,slw,slw1,slw2,nx,ny,nf,hzscl,hswgt,nfg,nfnf
  use control_vectors, only:  nrf_var
  use mpimod, only:  nvar_id
  implicit none

! Declare passed variables
  integer(i_kind)                        ,intent(in   ) :: nlevs
  real(r_kind),dimension(nval_lenz)      ,intent(inout) :: z
  real(r_kind),dimension(nlat,nlon,nlevs),intent(inout) :: work

! Declare local variables
  integer(i_kind) j,i
  integer(i_kind) k,iz,kk,kkk

  real(r_kind),dimension(nhscrf):: totwgt
  real(r_kind),allocatable,dimension(:,:):: zloc,pall


! Regional case
  if(regional)then
!$omp parallel do  schedule(dynamic,1) private(k,i,j,iz,totwgt,zloc), &
!$omp  shared(nlevs,nhscrf,hswgt,hzscl,nrf_var,nvar_id,nlat,nlon,nnnn1o,slw,ii,jj,z,nx,ny,work),default(none)
     do k=1,nlevs

        allocate(zloc(nlat*nlon,nhscrf))
!       apply horizontal recursive filters
        do j=1,nhscrf
           totwgt(j)=sqrt(hswgt(j)*hzscl(j)*hzscl(j))
        end do

        if(nrf_var(nvar_id(k))=='sf'.or.nrf_var(nvar_id(k))=='vp')then
           totwgt(3)=sqrt(half)*totwgt(3)
        end if

        call sqrt_rfxyyx_ad(zloc,work(1,1,k),ny,nx,ii(1,1,1,k),&
             jj(1,1,1,k),slw(1,k),totwgt)

        do j=1,nhscrf
           iz=nlat*nlon*(k-1)+nlat*nlon*nnnn1o*(j-1)
           do i=1,nlat*nlon
              z(i+iz)=zloc(i,j)
           end do
        end do
        deallocate(zloc)
        
     end do

! Global case
  else

     do j=1,nhscrf
        totwgt(j)=sqrt(hswgt(j))*hzscl(j)
     end do
     

!$omp parallel do  schedule(dynamic,1) private(kk) &
!$omp private(i,j,k,iz,kkk,pall,zloc), &
!$omp   shared(nlevs,nlon,nlat,nhscrf,nfnf,nx,ny,nnnn1o,z,ii,jj,ii1,jj1,slw,slw1, &
!$omp          ii2,jj2,slw2,nfg,nf,work,totwgt),default(none)
     do kk=1,3*nlevs

        k=(kk-1)/3+1
        kkk=mod(kk-1,3)+1

!       Recursive filter applications

        if(kkk == 1)then

!          First do equatorial/mid-latitude band
           allocate(pall(ny,nx),zloc(ny*nx,nhscrf))
           call grid2tr(work(1,1,k),pall)
           call sqrt_rfxyyx_ad(zloc,pall,ny,nx,ii(1,1,1,k),jj(1,1,1,k),slw(1,k),totwgt)
           do j=1,nhscrf
              iz=(ny*nx+2*nfnf)*(k-1+nnnn1o*(j-1))
              do i=1,ny*nx
                 z(i+iz)=z(i+iz)+zloc(i,j)
              end do
           end do
           deallocate(pall,zloc)

        else if(kkk == 2)then

!          North pole patch --interpolate - recursive filter - adjoint interpolate

           allocate(pall(-nf:nf,-nf:nf),zloc(nfnf,nhscrf))
           call grid2nh(work(1,1,k),pall)
           call sqrt_rfxyyx_ad(zloc,pall,nfg,nfg,ii1(1,1,1,k),jj1(1,1,1,k),slw1(1,k),totwgt)
           do j=1,nhscrf
              iz=(ny*nx+2*nfnf)*(k-1+nnnn1o*(j-1))+ny*nx
              do i=1,nfnf
                 z(i+iz)=z(i+iz)+zloc(i,j)
              end do
           end do
           deallocate(pall,zloc)

        else if(kkk == 3)then

!          South pole patch --interpolate - recursive filter - adjoint interpolate

           allocate(pall(-nf:nf,-nf:nf),zloc(nfnf,nhscrf))
           call grid2sh(work(1,1,k),pall)
           call sqrt_rfxyyx_ad(zloc,pall,nfg,nfg,ii2(1,1,1,k),jj2(1,1,1,k),slw2(1,k),totwgt)

           do j=1,nhscrf
              iz=(ny*nx+2*nfnf)*(k-1+nnnn1o*(j-1))+ny*nx+nfnf
              do i=1,nfnf
                 z(i+iz)=z(i+iz)+zloc(i,j)
              end do
           end do
           deallocate(pall,zloc)
        end if

!    End of k loop over nlevs
     end do

! End of global block
  end if

  return
end subroutine sqrt_smoothrf_ad
! ------------------------------------------------------------------------------
subroutine sqrt_rfxyyx(z,p1,nx,ny,iix,jjx,dssx,totwgt)
  
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    sqrt_rfxyyx   sqrt perform horizontal smoothing
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: smoothrf perform self-adjoint horizontal smoothing. nsloop
!           smoothing fields.
!
! program history log:
!   2000-03-15  wu
!   2004-08-24  derber - change indexing add rfhyt to speed things up
!
!   input argument list:
!     p1       - horizontal field to be smoothed
!     nx       - first dimension of p1
!     ny       - second dimension of p1
!     iix      - array of pointers for smoothing table (first dimension)
!     jjx      - array of pointers for smoothing table (second dimension)
!     dssx     - renormalization constants including variance
!     wgt      - weight (empirical*expected)
!
!   output argument list:
!                 all after horizontal smoothing
!     p1       - horizontal field which has been smoothed
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only:  zero
  use berror, only: be,table,ndeg,nhscrf
  implicit none

! Declare passed variables
  integer(i_kind)                     ,intent(in   ) :: nx,ny
  integer(i_kind),dimension(nx,ny,nhscrf),intent(in   ) :: iix,jjx
  real(r_kind)   ,dimension(nx,ny,3)  ,intent(in   ) :: z
  real(r_kind)   ,dimension(nx,ny)    ,intent(  out) :: p1
  real(r_kind)   ,dimension(nx,ny)    ,intent(in   ) :: dssx
  real(r_kind)   ,dimension(nhscrf)      ,intent(in   ) :: totwgt

! Declare local variables
  integer(i_kind) ix,iy,i,j,im,n

  real(r_kind),dimension(nx,ny):: p2,p1t
  real(r_kind),dimension(nx,ny,ndeg):: alx,aly

! Zero local arrays
  do iy=1,ny
     do ix=1,nx
        p1(ix,iy)=zero
     enddo
  enddo

! Loop over number of scales
 
  do n=1,nhscrf

     do im=1,ndeg
        do j=1,ny
           do i=1,nx
              alx(i,j,im)=table(iix(i,j,n),im)
              aly(i,j,im)=table(jjx(i,j,n),im)
           enddo
        enddo
     enddo

     do iy=1,ny
        do ix=1,nx
           p1t(ix,iy)=z(ix,iy,n)*sqrt(dssx(ix,iy))*totwgt(n)
        enddo
     enddo

     call rfhy(p1t,p2,nx,ny,ndeg,aly,be)

     call rfhx0(p2,p1,nx,ny,ndeg,alx,be)

! end loop over number of horizontal scales
  end do

  return
end subroutine sqrt_rfxyyx
! ------------------------------------------------------------------------------
subroutine sqrt_rfxyyx_ad(z,p1,nx,ny,iix,jjx,dssx,totwgt)
  
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rfxyyx      perform horizontal smoothing
!   prgmmr: wu               org: np22                date: 2000-03-15
!
! abstract: smoothrf perform self-adjoint horizontal smoothing. nsloop
!           smoothing fields.
!
! program history log:
!   2000-03-15  wu
!   2004-08-24  derber - change indexing add rfhyt to speed things up
!
!   input argument list:
!     p1       - horizontal field to be smoothed
!     nx       - first dimension of p1
!     ny       - second dimension of p1
!     iix      - array of pointers for smoothing table (first dimension)
!     jjx      - array of pointers for smoothing table (second dimension)
!     dssx     - renormalization constants including variance
!     wgt      - weight (empirical*expected)
!
!   output argument list:
!                 all after horizontal smoothing
!     p1       - horizontal field which has been smoothed
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only:  zero
  use berror, only: be,table,ndeg,nhscrf
  implicit none

! Declare passed variables
  integer(i_kind)                     ,intent(in   ) :: nx,ny
  integer(i_kind),dimension(nx,ny,nhscrf),intent(in   ) :: iix,jjx
  real(r_kind)   ,dimension(nx,ny,3)  ,intent(  out) :: z
  real(r_kind)   ,dimension(nx,ny)    ,intent(inout) :: p1
  real(r_kind)   ,dimension(nx,ny)    ,intent(in   ) :: dssx
  real(r_kind)   ,dimension(nhscrf)      ,intent(in   ) :: totwgt

! Declare local variables
  integer(i_kind) ix,iy,i,j,im,n

  real(r_kind),dimension(nx,ny):: p2,p1t
  real(r_kind),dimension(nx,ny,ndeg):: alx,aly

! Loop over number of scales
 
  do n=1,nhscrf

     do iy=1,ny
        do ix=1,nx
           p2(ix,iy)=zero
        enddo
     enddo
     do im=1,ndeg
        do j=1,ny
           do i=1,nx
              alx(i,j,im)=table(iix(i,j,n),im)
              aly(i,j,im)=table(jjx(i,j,n),im)
           enddo
        enddo
     enddo

     call rfhx0(p1,p2,nx,ny,ndeg,alx,be)

     call rfhyt(p2,p1t,nx,ny,ndeg,aly,be)


     do iy=1,ny
        do ix=1,nx
           z(ix,iy,n)=p1t(ix,iy)*sqrt(dssx(ix,iy))*totwgt(n)
        enddo
     enddo

! end loop over number of horizontal scales
  end do

  return
end subroutine sqrt_rfxyyx_ad
! ------------------------------------------------------------------------------