subroutine anbkerror_reg(gradx,grady,mype)
!-----------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    anbkerror_reg apply anisotropic background error covariance            
!   prgmmr: parrish          org: np22                date: 2005-02-03
!
! abstract: apply regional anisotropic background error.
!
! program history log:
!   2005-02-08  parrish
!   2005-04-29  parrish - replace coarse2fine with fgrid2agrid;
!                         remove ansmoothrf_reg_d
!
!   input argument list:
!     gradx    - input field  
!
!   output
!     grady    - background structure * gradx 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2
  use jfunc, only: nclen,nt,np,nq,nst,nvp,noz,nsst,ncw,nrclen,nclen1
  use balmod, only: balance,tbalance
  use berror, only: varprd
  use constants, only: zero
  implicit none

! Declare passed variables
  integer(i_kind),intent(in):: mype
  real(r_kind),dimension(nclen),intent(inout):: gradx
  real(r_kind),dimension(nclen),intent(out):: grady

! Declare local variables
  integer(i_kind) i,j
  real(r_kind),dimension(lat2,lon2):: sst,slndt,sicet

! Put things in grady first since operations change input variables
  do i=1,nclen
     grady(i)=gradx(i)
  end do

! Zero arrays for land, ocean, ice skin (surface) temperature.
  do j=1,lon2
     do i=1,lat2
        slndt(i,j)=zero
        sst(i,j)  =zero
        sicet(i,j)=zero
     end do
  end do

! Transpose of balance equation
  call tbalance(grady(nt),grady(np),grady(nq),grady(nst),grady(nvp))

! Apply background error covariance
  call anbkgcov_reg(grady(nst),grady(nvp),grady(nt),grady(np),grady(nq),&
       grady(noz),grady(nsst),sst,slndt,sicet,grady(ncw),mype)

! Balance equation
  call balance(grady(nt),grady(np),grady(nq),grady(nst),grady(nvp))

! Take care of background error for bias correction terms
  do i=1,nrclen
     grady(i+nclen1)=grady(i+nclen1)*varprd(i)
  end do

end subroutine anbkerror_reg

subroutine anbkgcov_reg(st,vp,t,p,q,oz,skint,sst,slndt,sicet,cwmr,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    anbkgcov    apply anisotropic background error covar
!   prgmmr: parrish        org: np22                date: 2005-02-14
!
! abstract: apply regional anisotropic background error covariance
!
! program history log:
!   2005-02-14  parrish
!
!   input argument list:
!     t        - t on subdomain
!     p        - p surface pressure on subdomain
!     q        - q on subdomain
!     oz       - ozone on subdomain
!     skint    - skin temperature on subdomain
!     sst      - sea surface temperature on subdomain
!     slndt    - land surface temperature on subdomain
!     sicet    - ice surface temperature on subdomain
!     cwmr     - cloud water mixing ratio on subdomain
!     st       - streamfunction on subdomain
!     vp       - velocity potential on subdomain
!
!   output argument list:
!                 all after smoothing, combining scales
!     t        - t on subdomain
!     p        - p surface pressure on subdomain
!     q        - q on subdomain
!     oz       - ozone on subdomain
!     skint    - skin temperature on subdomain
!     sst      - sea surface temperature on subdomain
!     slndt    - land surface temperature on subdomain
!     sicet    - ice surface temperature on subdomain
!     cwmr     - cloud water mixing ratio on subdomain
!     st       - streamfunction on subdomain
!     vp       - velocity potential on subdomain
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2,nlat,nlon,nsig,nsig1o
  use anberror, only: rtma_subdomain_option
  use constants, only: zero
  implicit none

! Passed Variables
  integer(i_kind),intent(in):: mype
  real(r_kind),dimension(lat2,lon2),intent(inout):: p,skint,sst,slndt,sicet
  real(r_kind),dimension(lat2,lon2,nsig),intent(inout):: t,q,cwmr,oz,st,vp

! Local Variables
  integer(i_kind) iflg
  real(r_kind),dimension(nlat,nlon,nsig1o):: hwork

  integer(i_kind) nsmooth,nsmooth_shapiro
  nsmooth=0! hardwired for
  nsmooth_shapiro=1 !testing purposes. Will change to namelist variable

! break up skin temp into components
  call anbkgvar_reg(skint,sst,slndt,sicet,0)

! Perform simple vertical smoothing while fields are in sudomain mode.
! The accompanying smoothing in the horizontal is performed inside the
! recursive filter. Motivation: Reduce possible high frequency noise in
! the analysis that would arise from the use of a "non-blending" RF algorithm.
                   

  call vert_smther(t,nsmooth,nsmooth_shapiro)                 
  call vert_smther(q,nsmooth,nsmooth_shapiro)                 
  call vert_smther(oz,nsmooth,nsmooth_shapiro)                 
  call vert_smther(cwmr,nsmooth,nsmooth_shapiro)                 
  call vert_smther(st,nsmooth,nsmooth_shapiro)                 
  call vert_smther(vp,nsmooth,nsmooth_shapiro)                 

  if(rtma_subdomain_option) then

    oz=zero
    cwmr=zero
    sst=zero
    slndt=zero
    sicet=zero
    call ansmoothrf_reg_subdomain_option(t,p,q,st,vp,mype)

  else

!   Convert from subdomain to full horizontal field distributed among processors
    iflg=1
    call sub2grid(hwork,t,p,q,oz,sst,slndt,sicet,cwmr,st,vp,iflg)

!   Apply horizontal smoother for number of horizontal scales
    call ansmoothrf_reg(hwork,mype)

!   Put back onto subdomains
    call grid2sub(hwork,t,p,q,oz,sst,slndt,sicet,cwmr,st,vp)

  end if

  call tvert_smther(vp,nsmooth,nsmooth_shapiro)                 
  call tvert_smther(st,nsmooth,nsmooth_shapiro)                 
  call tvert_smther(cwmr,nsmooth,nsmooth_shapiro)                 
  call tvert_smther(oz,nsmooth,nsmooth_shapiro)                 
  call tvert_smther(q,nsmooth,nsmooth_shapiro)                 
  call tvert_smther(t,nsmooth,nsmooth_shapiro)                 

! combine sst,sldnt, and sicet into skin temperature field
  call anbkgvar_reg(skint,sst,slndt,sicet,1)

end subroutine anbkgcov_reg

subroutine anbkgvar_reg(skint,sst,slndt,sicet,iflg)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    anbkgvar_reg      manipulate skin temp
!   prgmmr: parrish          org: np22                date: 1990-10-06
!
! abstract: manipulate skin temp <--> sst,sfc temp, and ice temp fields
!
! program history log:
!   2005-01-22  parrish
!
!   input argument list:
!     skint    - skin temperature grid values
!     sst      - sst grid values
!     slndt    - land surface temperature grid values
!     sicet    - snow/ice covered surface temperature grid values
!     iflg     - flag for skin temperature manipulation
!                0: skint --> sst,slndt,sicet
!                1: sst,slndt,sicet --> skint
!
!   output argument list:
!     skint    - skin temperature grid values
!     sst      - sst grid values
!     slndt    - land surface temperature grid values
!     sicet    - snow/ice covered surface temperature grid values
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only:  one
  use gridmod, only: lat2,lon2
  use guess_grids, only: ntguessfc,isli
  implicit none

! Declare passed variables
  integer(i_kind),intent(in):: iflg
  real(r_kind),dimension(lat2,lon2),intent(inout):: skint,sst,slndt,sicet

! Declare local variables
  integer(i_kind) i,j

       do j=1,lon2
         do i=1,lat2
           if(iflg == 0) then
! Break skin temperature into components
!          If land point
             if(isli(i,j,ntguessfc) == 1) then
               slndt(i,j)=skint(i,j)
!          If ice
             else if(isli(i,j,ntguessfc) == 2) then
               sicet(i,j)=skint(i,j)
!          Else treat as a water point
             else
               sst(i,j)=skint(i,j)
             end if

           else if (iflg.eq.1) then
! Combine sst,slndt, and sicet into skin temperature field
!          Land point, load land sfc t into skint
             if(isli(i,j,ntguessfc) == 1) then
               skint(i,j)=slndt(i,j)
!          Ice, load ice temp into skint
             else if(isli(i,j,ntguessfc) == 2) then
               skint(i,j)=sicet(i,j)
!          Treat as a water point, load sst into skint
             else
               skint(i,j)=sst(i,j)
             end if
           end if
         end do
       end do

  return
end subroutine anbkgvar_reg

subroutine ansmoothrf_reg(work,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ansmoothrf_reg  anisotropic rf for regional mode
!   prgmmr: parrish          org: np22                date: 2005-02-14
!
! abstract: apply anisotropic rf for regional mode
!
! program history log:
!   2005-02-14  parrish
!
!   input argument list:
!     work     - fields to be smoothed
!
!   output argument list:
!     work     - smoothed fields
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind,r_single
  use anberror, only: ids,ide,jds,jde,kds,kde,ips,ipe,jps,jpe,kps,kpe, &
                                   filter_all,ngauss
  use mpimod, only: npe
  use constants, only: zero
  use gridmod, only: nlat,nlon
  use fgrid2agrid_mod, only: fgrid2agrid,tfgrid2agrid
  use raflib, only: raf4_ad,raf4
  implicit none

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

! Declare local variables
  integer(i_kind) i,igauss,j,k
  real(r_kind),dimension(ips:ipe,jps:jpe,kps:kpe):: worka
  real(r_single),dimension(ngauss,ips:ipe,jps:jpe,kps:kpe):: workb


!  adjoint of coarse to fine grid
  do k=kps,kpe
    call tfgrid2agrid(work(1,1,k),worka(ips,jps,k))

  end do

!  transfer coarse grid fields to ngauss copies
  do k=kps,kpe
    do j=jps,jpe
      do i=ips,ipe
        do igauss=1,ngauss
          workb(igauss,i,j,k)=worka(i,j,k)
        end do
      end do
    end do
  end do

!   apply recursive filter

  call raf4(workb,filter_all,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npe)
  call raf4_ad(workb,filter_all,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npe)

!  add together ngauss copies
  do k=kps,kpe
    do j=jps,jpe
      do i=ips,ipe
        worka(i,j,k)=zero
        do igauss=1,ngauss
          worka(i,j,k)=worka(i,j,k)+workb(igauss,i,j,k)
        end do
      end do
    end do
  end do

!  coarse to fine grid
  do k=kps,kpe
    call fgrid2agrid(worka(ips,jps,k),work(1,1,k))
  end do


end subroutine ansmoothrf_reg

subroutine vert_smther(g,nsmooth,nsmooth_shapiro)

! nsmooth > 0 ==> apply 1-2-1 smoother
! nsmooth_shapiro 0 ==> apply second moment preserving  
!                       "shapiro smoother"

  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2,nsig
  implicit none


! Declare passed variables
  integer(i_kind),intent(in):: nsmooth,nsmooth_shapiro
  real(r_kind),dimension(1:lat2,1:lon2,1:nsig),intent(inout):: g
 
! Declare local variables
  integer(i_kind) i,j,l,k,kp,km,kp3,km3
  real(r_kind), allocatable:: gaux(:)
 
  if (nsig==1)return

  allocate(gaux(1:nsig))

  if (nsmooth > 0 ) then
    do i=1,lat2
    do j=1,lon2
     do l=1,nsmooth
       gaux(1:nsig)=g(i,j,1:nsig)
       do k=1,nsig
         kp=min(k+1,nsig) ; km=max(1,k-1)
         g(i,j,k)=.25*(gaux(kp)+gaux(km))+.5*gaux(k)                
       enddo
     enddo
    enddo 
    enddo 
  endif

  if (nsmooth_shapiro > 0 .and. nsmooth <= 0) then
    do i=1,lat2
    do j=1,lon2
     do l=1,nsmooth_shapiro
       gaux(1:nsig)=g(i,j,1:nsig)
       do k=1,nsig
         kp=min(k+1,nsig) ; km=max(1,k-1)
         kp3=min(k+3,nsig) ; km3=max(1,k-3)
         g(i,j,k)=.28125*(gaux(kp)+gaux(km))+.5*gaux(k)-.03125*(gaux(kp3)+gaux(km3))
       enddo
     enddo
    enddo 
    enddo 
  endif
  deallocate(gaux)

  return
end subroutine vert_smther

subroutine tvert_smther(g,nsmooth,nsmooth_shapiro)

! adjoint of tvert_smther
! nsmooth > 0 ==>  1-2-1 smoother
! nsmooth_shapiro 0 ==>  second moment preserving  
!                       "shapiro smoother"

  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2,nsig
  implicit none


! Declare passed variables
  integer(i_kind),intent(in):: nsmooth,nsmooth_shapiro
  real(r_kind),dimension(1:lat2,1:lon2,1:nsig),intent(inout):: g
 
! Declare local variables
  integer(i_kind) i,j,l,k,kp,km,kp3,km3
  real(r_kind), allocatable:: gaux(:)
 
  if (nsig==1)return

  allocate(gaux(1:nsig))

  if (nsmooth > 0 ) then
    do i=1,lat2
    do j=1,lon2
     do l=1,nsmooth
       gaux(1:nsig)=0._r_kind
       do k=1,nsig
         kp=min(k+1,nsig) ; km=max(1,k-1)
         gaux(k)=gaux(k)+.5*g(i,j,k)
         gaux(km)=gaux(km)+.25*g(i,j,k)
         gaux(kp)=gaux(kp)+.25*g(i,j,k)
       enddo
       g(i,j,1:nsig)=gaux(1:nsig)
     enddo
    enddo 
    enddo 
  endif

  if (nsmooth_shapiro > 0 .and. nsmooth <= 0) then
    do i=1,lat2
    do j=1,lon2
     do l=1,nsmooth_shapiro
       gaux(1:nsig)=0._r_kind
       do k=1,nsig
         kp=min(k+1,nsig) ; km=max(1,k-1)
         kp3=min(k+3,nsig) ; km3=max(1,k-3)
         gaux(km3)=gaux(km3)-.03125*g(i,j,k)
         gaux(kp3)=gaux(kp3)-.03125*g(i,j,k)
         gaux(k)=gaux(k)+.5*g(i,j,k)
         gaux(km)=gaux(km)+.28125*g(i,j,k)
         gaux(kp)=gaux(kp)+.28125*g(i,j,k)
       enddo
       g(i,j,1:nsig)=gaux(1:nsig)
     enddo
    enddo 
    enddo 
  endif
  deallocate(gaux)

  return
end subroutine tvert_smther

subroutine ansmoothrf_reg_subdomain_option(t,p,q,st,vp,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ansmoothrf_reg_subdomain_option  anisotropic rf for regional mode
!   prgmmr: parrish          org: np22                date: 2005-02-14
!
! abstract: apply anisotropic rf for regional mode (using subdomains instead of slabs)
!              NOTE: only works if filter grid is same as analysis grid
!
! program history log:
!   2005-02-14  parrish
!
!   input argument list:
!     t,p,q,oz,st,slndt,sicet,cwmr,st,vp   -  fields to be smoothed
!
!   output argument list:
!     t,p,q,oz,st,slndt,sicet,cwmr,st,vp   -  smoothed fields
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind,r_single
  use anberror, only: ids,ide,jds,jde,kds,kde,ips,ipe,jps,jpe,kps,kpe, &
                        filter_all,ngauss,halo_update_reg
  use mpimod, only: npe
  use constants, only: zero
  use gridmod, only: lat2,lon2,istart,jstart,nsig
  use raflib, only: raf4_ad,raf4
  implicit none

! Declare passed variables
  integer(i_kind), intent(in):: mype
  real(r_kind),dimension(lat2,lon2),intent(inout):: p
  real(r_kind),dimension(lat2,lon2,nsig),intent(inout):: t,q,vp,st

! Declare local variables
  integer(i_kind) i,igauss,iloc,j,jloc,k,kk,mm1
  real(r_single),dimension(ngauss,ips:ipe,jps:jpe,kps:kpe):: workb

  mm1=mype+1

!  transfer variables to ngauss copies
  kk=0
  do k=1,nsig
    kk=kk+1
    do j=jps,jpe
      jloc=j-jstart(mm1)+2
      do i=ips,ipe
        iloc=i-istart(mm1)+2
        do igauss=1,ngauss
          workb(igauss,i,j,kk)=st(iloc,jloc,k)
        end do
      end do
    end do
  end do
  do k=1,nsig
    kk=kk+1
    do j=jps,jpe
      jloc=j-jstart(mm1)+2
      do i=ips,ipe
        iloc=i-istart(mm1)+2
        do igauss=1,ngauss
          workb(igauss,i,j,kk)=vp(iloc,jloc,k)
        end do
      end do
    end do
  end do
  kk=kk+1
  do j=jps,jpe
    jloc=j-jstart(mm1)+2
    do i=ips,ipe
      iloc=i-istart(mm1)+2
      do igauss=1,ngauss
        workb(igauss,i,j,kk)=p(iloc,jloc)
      end do
    end do
  end do
  do k=1,nsig
    kk=kk+1
    do j=jps,jpe
      jloc=j-jstart(mm1)+2
      do i=ips,ipe
        iloc=i-istart(mm1)+2
        do igauss=1,ngauss
          workb(igauss,i,j,kk)=t(iloc,jloc,k)
        end do
      end do
    end do
  end do
  do k=1,nsig
    kk=kk+1
    do j=jps,jpe
      jloc=j-jstart(mm1)+2
      do i=ips,ipe
        iloc=i-istart(mm1)+2
        do igauss=1,ngauss
          workb(igauss,i,j,kk)=q(iloc,jloc,k)
        end do
      end do
    end do
  end do

!   apply recursive filter

  call raf4(workb,filter_all,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npe)
  call raf4_ad(workb,filter_all,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npe)

!  add together ngauss copies
  kk=0
  do k=1,nsig
    kk=kk+1
    do j=jps,jpe
      jloc=j-jstart(mm1)+2
      do i=ips,ipe
        iloc=i-istart(mm1)+2
        st(iloc,jloc,k)=zero
        do igauss=1,ngauss
          st(iloc,jloc,k)=st(iloc,jloc,k)+workb(igauss,i,j,kk)
        end do
      end do
    end do
  end do
  do k=1,nsig
    kk=kk+1
    do j=jps,jpe
      jloc=j-jstart(mm1)+2
      do i=ips,ipe
        iloc=i-istart(mm1)+2
        vp(iloc,jloc,k)=zero
        do igauss=1,ngauss
          vp(iloc,jloc,k)=vp(iloc,jloc,k)+workb(igauss,i,j,kk)
        end do
      end do
    end do
  end do
  kk=kk+1
  do j=jps,jpe
    jloc=j-jstart(mm1)+2
    do i=ips,ipe
      iloc=i-istart(mm1)+2
      p(iloc,jloc)=zero
      do igauss=1,ngauss
        p(iloc,jloc)=p(iloc,jloc)+workb(igauss,i,j,kk)
      end do
    end do
  end do
  do k=1,nsig
    kk=kk+1
    do j=jps,jpe
      jloc=j-jstart(mm1)+2
      do i=ips,ipe
        iloc=i-istart(mm1)+2
        t(iloc,jloc,k)=zero
        do igauss=1,ngauss
          t(iloc,jloc,k)=t(iloc,jloc,k)+workb(igauss,i,j,kk)
        end do
      end do
    end do
  end do
  do k=1,nsig
    kk=kk+1
    do j=jps,jpe
      jloc=j-jstart(mm1)+2
      do i=ips,ipe
        iloc=i-istart(mm1)+2
        q(iloc,jloc,k)=zero
        do igauss=1,ngauss
          q(iloc,jloc,k)=q(iloc,jloc,k)+workb(igauss,i,j,kk)
        end do
      end do
    end do
  end do

!   update halos:
  call halo_update_reg(p,1)
  call halo_update_reg(t,nsig)
  call halo_update_reg(q,nsig)
  call halo_update_reg(vp,nsig)
  call halo_update_reg(st,nsig)

end subroutine ansmoothrf_reg_subdomain_option