subroutine anbkerror(gradx,grady)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    anbkerror 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
!   2006-11-30  todling - add fpsproj as arg to (t)balance routine(s)
!   2008-10-10  derber - add strong constraint to background error
!   2008-12-29  todling - update interface to strong_bk/bk_ad
!   2009-04-13  derber - move strong_bk into balance
!   2009-07-01  sato - update for global mode
!   2010-03-09  zhu    - change anbkgcov interface for generalized control variables
!                      - make changes to interfaces of sub2grid and grid2sub
!
!   input argument list:
!     gradx    - input field  
!
!   output
!     grady    - background structure * gradx 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2
  use jfunc, only: nsclen,npclen
  use balmod, only: balance,tbalance
  use berror, only: varprd,fpsproj
  use constants, only: izero,zero
  use control_vectors
  use gsi_4dvar, only: nsubwin
  implicit none

! Declare passed variables
  type(control_vector),intent(inout) :: gradx
  type(control_vector),intent(inout) :: grady

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

! Put things in grady first since operations change input variables
  grady=gradx

! 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

! Loop on control steps
  do ii=1,nsubwin

!    Transpose of balance equation
     call tbalance(grady%step(ii)%t ,grady%step(ii)%p , &
                   grady%step(ii)%st,grady%step(ii)%vp,fpsproj)

!    Apply variances, as well as vertical & horizontal parts of background error
     call anbkgcov(grady%step(ii),sst,slndt,sicet)

!    Balance equation
     call balance(grady%step(ii)%t ,grady%step(ii)%p ,&
                  grady%step(ii)%st,grady%step(ii)%vp,fpsproj)

  end do

! Take care of background error for bias correction terms
  if(nsclen>izero)then
     do i=1,nsclen
        grady%predr(i)=grady%predr(i)*varprd(i)
     end do
  end if
  if(npclen>izero)then
     do i=1,npclen
        grady%predp(i)=grady%predp(i)*varprd(nsclen+i)
     end do
  end if

end subroutine anbkerror


subroutine anbkgcov(cstate,sst,slndt,sicet)
!$$$  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
!   2009-07-01  sato - update for global mode
!   2010-03-09  zhu     - make changes for generalizing control vectors
!                       - replace explicit use of each control variable
!                         by a control_state 'cstate'
!                       - use nrf* for generalized control variables
!
!   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,latlon11
  use anberror, only: rtma_subdomain_option,nsmooth, nsmooth_shapiro
  use control_vectors, only: nrf_levb,control_state,nrf3,nrf3_loc,nrf2_sst
  use constants, only: izero,ione,zero
  implicit none

! Passed Variables
  type(control_state),intent(inout) :: cstate
  real(r_kind),dimension(lat2,lon2)     ,intent(inout) :: sst,slndt,sicet

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


! break up skin temp into components
  if (nrf2_sst>izero) call anbkgvar(cstate%sst,sst,slndt,sicet,izero)

! 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.

  
  do k=1,nrf3
     loc=nrf3_loc(k)
     nk=(nrf_levb(loc)-ione)*latlon11+ione
     call vert_smther(cstate%values(nk),nsmooth,nsmooth_shapiro)
  end do

  if(rtma_subdomain_option) then

     cstate%oz=zero
     cstate%cw=zero
     sst=zero
     slndt=zero
     sicet=zero
     call ansmoothrf_reg_subdomain_option(cstate%values)

  else

! Convert from subdomain to full horizontal field distributed among processors
     iflg=ione
     call sub2grid(hwork,cstate,sst,slndt,sicet,iflg) 

! Apply horizontal smoother for number of horizontal scales
     call ansmoothrf(hwork)

! Put back onto subdomains
     call grid2sub(hwork,cstate,sst,slndt,sicet)

  end if

  do k=nrf3,1,-1
     loc=nrf3_loc(k)
     nk=(nrf_levb(loc)-ione)*latlon11+ione
     call tvert_smther(cstate%values(nk),nsmooth,nsmooth_shapiro)
  end do

! combine sst,sldnt, and sicet into skin temperature field
  if (nrf2_sst>izero) call anbkgvar(cstate%sst,sst,slndt,sicet,ione)

end subroutine anbkgcov


subroutine anbkgvar(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
!   2008-06-05  safford - rm unused uses
!
!   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 constants, only: izero,ione
  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2
  use guess_grids, only: isli2
  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 == izero) then
! Break skin temperature into components
!          If land point
                if(isli2(i,j) == ione) then
                   slndt(i,j)=skint(i,j)
!          If ice
                else if(isli2(i,j) == 2_i_kind) 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==ione) then
! Combine sst,slndt, and sicet into skin temperature field
!          Land point, load land sfc t into skint
                if(isli2(i,j) == ione) then
                   skint(i,j)=slndt(i,j)
!          Ice, load ice temp into skint
                else if(isli2(i,j) == 2_i_kind) 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


subroutine ansmoothrf(work)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ansmoothrf  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
!   2008-12-04  sato - update for global mode
!   2009-01-02  todling - get mype from mpimod directly
!
!   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: indices,indices_p,ngauss,pf2aP1, &
                      filter_all,filter_p2,filter_p3
  use patch2grid_mod, only: patch2grid, tpatch2grid
  use mpimod, only:  npe
  use constants, only: zero
  use gridmod, only: nlat,nlon,regional
  use fgrid2agrid_mod, only: fgrid2agrid,tfgrid2agrid
  use raflib, only: raf4_ad_wrap,raf4_wrap
  implicit none

! Declare passed variables
  real(r_kind),dimension(nlat,nlon,indices%kps:indices%kpe),intent(inout) :: work

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

  real(r_kind)  ,dimension( indices%ips:indices%ipe, &
                            indices%jps:indices%jpe, &
                            indices%kps:indices%kpe ):: worka
  real(r_single),dimension(ngauss, &
                           indices%ips:indices%ipe, &
                           indices%jps:indices%jpe, &
                           indices%kps:indices%kpe ):: workb

  real(r_kind)  ,allocatable,dimension(:,:,:)  :: workanp,workasp
  real(r_single),allocatable,dimension(:,:,:,:):: workbnp,workbsp

  if(.not.regional) then
     allocate(workanp(indices_p%ips:indices_p%ipe, &
                      indices_p%jps:indices_p%jpe, &
                      indices_p%kps:indices_p%kpe ))
     allocate(workbnp(ngauss, &
                      indices_p%ips:indices_p%ipe, &
                      indices_p%jps:indices_p%jpe, &
                      indices_p%kps:indices_p%kpe ))
     allocate(workasp(indices_p%ips:indices_p%ipe, &
                      indices_p%jps:indices_p%jpe, &
                      indices_p%kps:indices_p%kpe ))
     allocate(workbsp(ngauss, &
                      indices_p%ips:indices_p%ipe, &
                      indices_p%jps:indices_p%jpe, &
                      indices_p%kps:indices_p%kpe ))
  end if

!  adjoint of coarse to fine grid
  do k=indices%kps,indices%kpe
     if(regional) then
        call tfgrid2agrid(pf2aP1,work(1,1,k),worka(indices%ips,indices%jps,k))
     else
        call tpatch2grid(work(1,1,k), &
                         worka  (indices%ips,  indices%jps,k),   &
                         workanp(indices_p%ips,indices_p%jps,k), &
                         workasp(indices_p%ips,indices_p%jps,k))
     end if
  end do

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

  if(.not.regional) then
     do k=indices_p%kps,indices_p%kpe
        do j=indices_p%jps,indices_p%jpe
           do i=indices_p%ips,indices_p%ipe
              do igauss=1,ngauss
                 workbnp(igauss,i,j,k)=workanp(i,j,k)
                 workbsp(igauss,i,j,k)=workasp(i,j,k)
              end do
           end do
        end do
     end do
  end if

!   apply recursive filter

  call raf4_wrap(   workb,filter_all,ngauss,indices,npe)
  call raf4_ad_wrap(workb,filter_all,ngauss,indices,npe)

  if(.not.regional) then
     call raf4_wrap(   workbnp,filter_p2,ngauss,indices_p,npe)
     call raf4_ad_wrap(workbnp,filter_p2,ngauss,indices_p,npe)

     call raf4_wrap(   workbsp,filter_p3,ngauss,indices_p,npe)
     call raf4_ad_wrap(workbsp,filter_p3,ngauss,indices_p,npe)
  end if

!  add together ngauss copies
  worka=zero
  do k=indices%kps,indices%kpe
     do j=indices%jps,indices%jpe
        do i=indices%ips,indices%ipe
           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

  if(.not.regional) then
     workanp=zero
     workasp=zero
     do k=indices_p%kps,indices_p%kpe
        do j=indices_p%jps,indices_p%jpe
           do i=indices_p%ips,indices_p%ipe
              do igauss=1,ngauss
                 workanp(i,j,k)=workanp(i,j,k)+workbnp(igauss,i,j,k)
                 workasp(i,j,k)=workasp(i,j,k)+workbsp(igauss,i,j,k)
              end do
           end do
        end do
     end do
  end if

!  coarse to fine grid
  do k=indices%kps,indices%kpe
     if(regional) then
        call fgrid2agrid(pf2aP1,worka(indices%ips,indices%jps,k),work(1,1,k))
     else
        call patch2grid(work(1,1,k), &
                        worka  (indices%ips  ,indices%jps  ,k),&
                        workanp(indices_p%ips,indices_p%jps,k),&
                        workasp(indices_p%ips,indices_p%jps,k))
     end if
  end do

  if(.not.regional) then
     deallocate(workanp)
     deallocate(workbnp)
     deallocate(workasp)
     deallocate(workbsp)
  end if

end subroutine ansmoothrf

subroutine vert_smther(g,nsmooth,nsmooth_shapiro)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    vert_smther
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-15  lueken - added subprogram doc block
!
!   input argument list:
!    nsmooth
!    nsmooth_shapiro
!    g
!
!   output argument list:
!    g
!
! Notes:  nsmooth > 0 ==> apply 1-2-1 smoother
!         nsmooth_shapiro 0 ==> apply second moment preserving
!                               "shapiro smoother"
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use constants, only: izero,ione,quarter,half
  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==ione)return

  allocate(gaux(1:nsig))

  if (nsmooth > izero ) 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+ione,nsig) ; km=max(ione,k-ione)
                 g(i,j,k)=quarter*(gaux(kp)+gaux(km))+half*gaux(k)
              enddo
           enddo
        enddo
     enddo
  endif

  if (nsmooth_shapiro > izero .and. nsmooth <= izero) 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+ione,nsig) ; km=max(ione,k-ione)
                 kp3=min(k+3_i_kind,nsig) ; km3=max(ione,k-3_i_kind)
                 g(i,j,k)=.28125_r_kind*(gaux(kp)+gaux(km))+half*gaux(k)-.03125_r_kind*(gaux(kp3)+gaux(km3))
              enddo
           enddo
        enddo
     enddo
  endif
  deallocate(gaux)

  return
end subroutine vert_smther


subroutine tvert_smther(g,nsmooth,nsmooth_shapiro)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    tvert_smther
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-15  lueken - added subprogram doc block
!
!   input argument list:
!    nsmooth
!    nsmooth_shapiro
!    g
!
!   output argument list:
!    g
!
! Notes:  nsmooth > 0 ==>  1-2-1 smoother
!         nsmooth_shapiro 0 ==>  second moment preserving
!                               "shapiro smoother"
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use constants, only: izero,ione,zero,quarter,half
  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==ione)return

  allocate(gaux(1:nsig))

  if (nsmooth > izero ) then
     do i=1,lat2
        do j=1,lon2
           do l=1,nsmooth
              gaux(1:nsig)=zero
              do k=1,nsig
                 kp=min(k+ione,nsig) ; km=max(ione,k-ione)
                 gaux(k)=gaux(k)+half*g(i,j,k)
                 gaux(km)=gaux(km)+quarter*g(i,j,k)
                 gaux(kp)=gaux(kp)+quarter*g(i,j,k)
              enddo
              g(i,j,1:nsig)=gaux(1:nsig)
           enddo
        enddo
     enddo
  endif

  if (nsmooth_shapiro > izero .and. nsmooth <= izero) then
     do i=1,lat2
        do j=1,lon2
           do l=1,nsmooth_shapiro
              gaux(1:nsig)=zero
              do k=1,nsig
                 kp=min(k+ione,nsig) ; km=max(ione,k-ione)
                 kp3=min(k+3_i_kind,nsig) ; km3=max(ione,k-3_i_kind)
                 gaux(km3)=gaux(km3)-.03125_r_kind*g(i,j,k)
                 gaux(kp3)=gaux(kp3)-.03125_r_kind*g(i,j,k)
                 gaux(k)=gaux(k)+half*g(i,j,k)
                 gaux(km)=gaux(km)+.28125_r_kind*g(i,j,k)
                 gaux(kp)=gaux(kp)+.28125_r_kind*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(grady)
!$$$  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
!   2010-03-16  zhu - use an array grady for generalized control variables
!
!   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: indices, filter_all,ngauss,halo_update_reg
  use mpimod, only: mype,npe
  use constants, only: izero,ione,zero
  use gridmod, only: lat2,lon2,istart,jstart,nsig
  use raflib, only: raf4_ad_wrap,raf4_wrap
  use jfunc, only: nval_levs
  use control_vectors, only: nrf2,nrf3,nvars
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nval_levs):: grady

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

  integer(i_kind):: ids,ide,jds,jde,kds,kde,ips,ipe,jps,jpe,kps,kpe

  ids=indices%ids; ide=indices%ide
  jds=indices%jds; jde=indices%jde
  kds=indices%kds; kde=indices%kde
  ips=indices%ips; ipe=indices%ipe
  jps=indices%jps; jpe=indices%jpe
  kps=indices%kps; kpe=indices%kpe

  mm1=mype+ione

!  transfer variables to ngauss copies
  kk=izero
  do k=1,nval_levs
     kk=kk+ione
     do j=jps,jpe
        jloc=j-jstart(mm1)+2_i_kind
        do i=ips,ipe
           iloc=i-istart(mm1)+2_i_kind
           do igauss=1,ngauss
              workb(igauss,i,j,kk)=grady(iloc,jloc,k)
           end do
        end do
     end do
  end do

!   apply recursive filter

  call raf4_wrap(workb,filter_all,ngauss,indices,npe)
  call raf4_ad_wrap(workb,filter_all,ngauss,indices,npe)

!  add together ngauss copies
  kk=izero
  do k=1,nval_levs
     kk=kk+ione
     do j=jps,jpe
        jloc=j-jstart(mm1)+2_i_kind
        do i=ips,ipe
           iloc=i-istart(mm1)+2_i_kind
           grady(iloc,jloc,k)=zero
           do igauss=1,ngauss
              grady(iloc,jloc,k)=grady(iloc,jloc,k)+workb(igauss,i,j,kk)
           end do
        end do
     end do
  end do

!   update halos:
  call halo_update_reg(grady,nval_levs)

end subroutine ansmoothrf_reg_subdomain_option