subroutine frfhvo(p1,iv)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    frfhvo      performs vertical smoothing of fields
!   prgmmr: derber           org: np23                date: 2004-05-13
!
! abstract: performs vertical smoothing of fields
!
! program history log:
!   2004-05-13  derber, document
!   2005-01-22  parrish - make use of balmod and rename variables
!   2005-07-14  wu - add max bound to l2
!   2008-04-11  safford - rm unsed vars
!   2010-10-08  derber - optimize and clean up
!
!   input argument list:
!     p1       - input field to be smoothed (lat2,lon2,nsig)
!     iv       - location in alv for smoothing coefficients
!              - iv = 1 streamfunction
!              - iv = 2 velocity potential
!              - iv = 3 temperature        
!              - iv = 4 specific humidity  
!              - iv = 5 ozone           
!              - iv = 6 cloud condensate
!
!   output argument list
!     p1       - output field after vertical smoothing
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind
  use constants, only: zero,one
  use gridmod, only: nsig,regional,lat2,lon2
  use balmod, only: rllat1,llmax
  use berror, only: alv,be,ndeg
  implicit none

! lat2 = number of latitudes (lat2)
! lon2 = number of longitudes (lon2)
! nsig  = number of model levels 
! ndeg  = degree of smoothing (ndeg=4)

  integer(i_kind)                       ,intent(in   ) :: iv
  real(r_kind),dimension(lat2,lon2,nsig),intent(inout) :: p1

  integer(i_kind) j,i,k,kr,ki,l,l2
  real(r_kind) alvl,alvr,alvi
  real(r_kind) gaki,gakr,dekr,deki,bekr,beki
  real(r_kind),dimension(lat2,lon2,nsig):: p2
  real(r_kind),dimension(lat2,lon2,ndeg):: ga,de
  real(r_kind),dimension(lat2,lon2):: dl1,dl2


! Zero output array
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           p2(i,j,k)=p1(i,j,k)
           p1(i,j,k)=zero
        end do
     end do
  end do

! Zero local work arrays
  do k=1,ndeg
     do j=1,lon2
        do i=1,lat2
           ga(i,j,k)=zero
           de(i,j,k)=zero
        end do
     end do
  end do

  if (regional)then
     
     do j=1,lon2
        do k=1,lat2
           l=int(rllat1(k,j))
           l2=min0(l+1,llmax)
           dl2(k,j)=rllat1(k,j)-float(l)
           dl1(k,j)=one-dl2(k,j)
        end do
     end do
! Regional mode, odd degree
     if (mod(ndeg,2) == 1) then

!       advancing filter:
        do i=1,nsig
           do j=1,lon2
              do k=1,lat2
                 l=int(rllat1(k,j))
                 l2=min0(l+1,llmax)
                 alvl=dl1(k,j)*alv(l,1,i,iv)+dl2(k,j)*alv(l2,1,i,iv)
                 ga(k,j,1)=alvl*ga(k,j,1)+be(1)*p2(k,j,i)
                 p1(k,j,i)=ga(k,j,1)
              end do
           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 j=1,lon2
                 do k=1,lat2
                    l=int(rllat1(k,j))
                    l2=min0(l+1,llmax)
                    alvr=dl1(k,j)*alv(l,kr,i,iv)+dl2(k,j)*alv(l2,kr,i,iv)
                    alvi=dl1(k,j)*alv(l,ki,i,iv)+dl2(k,j)*alv(l2,ki,i,iv)
                    gakr=ga(k,j,kr)
                    gaki=ga(k,j,ki)
                    ga(k,j,kr)=alvr*gakr-alvi*gaki+bekr*p2(k,j,i)
                    ga(k,j,ki)=alvi*gakr+alvr*gaki+beki*p2(k,j,i)
                    p1(k,j,i)=p1(k,j,i)+ga(k,j,kr)
                 end do
              enddo
           enddo
        enddo
 
!       backing filter:
        do i=nsig,1,-1
           do j=1,lon2
              do k=1,lat2
                 p1(k,j,i)=p1(k,j,i)+de(k,j,1)
                 l=int(rllat1(k,j))
                 l2=min0(l+1,llmax)
                 alvl=dl1(k,j)*alv(l,1,i,iv)+dl2(k,j)*alv(l2,1,i,iv)
                 de(k,j,1)=alvl*(de(k,j,1)+be(1)*p2(k,j,i))
              end do
           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 j=1,lon2
                 do k=1,lat2
                    p1(k,j,i)=p1(k,j,i)+de(k,j,kr)
                    l=int(rllat1(k,j))
                    l2=min0(l+1,llmax)
                    alvr=dl1(k,j)*alv(l,kr,i,iv)+dl2(k,j)*alv(l2,kr,i,iv)
                    alvi=dl1(k,j)*alv(l,ki,i,iv)+dl2(k,j)*alv(l2,ki,i,iv)
                    dekr=de(k,j,kr)+bekr*p2(k,j,i)
                    deki=de(k,j,ki)+beki*p2(k,j,i)
                    de(k,j,kr)=alvr*dekr-alvi*deki
                    de(k,j,ki)=alvi*dekr+alvr*deki
                 end do
              enddo
           enddo
        enddo
 
! Regional mode, even degree
     else

        ! advancing filter:
        do i=1,nsig
           !       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 j=1,lon2
                 do k=1,lat2
                    l=int(rllat1(k,j))
                    l2=min0(l+1,llmax)
                    alvr=dl1(k,j)*alv(l,kr,i,iv)+dl2(k,j)*alv(l2,kr,i,iv)
                    alvi=dl1(k,j)*alv(l,ki,i,iv)+dl2(k,j)*alv(l2,ki,i,iv)
                    gakr=ga(k,j,kr)
                    gaki=ga(k,j,ki)
                    ga(k,j,kr)=alvr*gakr-alvi*gaki+bekr*p2(k,j,i)
                    ga(k,j,ki)=alvi*gakr+alvr*gaki+beki*p2(k,j,i)
                    p1(k,j,i)=p1(k,j,i)+ga(k,j,kr)
                 end do
              enddo
           enddo
        enddo
 
        !    backing filter:
        do i=nsig,1,-1
           !       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 j=1,lon2
                 do k=1,lat2
                    p1(k,j,i)=p1(k,j,i)+de(k,j,kr)
                    l=int(rllat1(k,j))
                    l2=min0(l+1,llmax)
                    alvr=dl1(k,j)*alv(l,kr,i,iv)+dl2(k,j)*alv(l2,kr,i,iv)
                    alvi=dl1(k,j)*alv(l,ki,i,iv)+dl2(k,j)*alv(l2,ki,i,iv)
                    dekr=de(k,j,kr)+bekr*p2(k,j,i)
                    deki=de(k,j,ki)+beki*p2(k,j,i)
                    de(k,j,kr)=alvr*dekr-alvi*deki
                    de(k,j,ki)=alvi*dekr+alvr*deki
                 end do
              enddo
           enddo
        enddo
     endif

! Global branch
  else
!    Global mode, odd degree
     if (mod(ndeg,2) == 1) then

!       advancing filter:
        do i=1,nsig
           do j=1,lon2
              do k=1,lat2
                 ga(k,j,1)=alv(k,1,i,iv)*ga(k,j,1)+be(1)*p2(k,j,i)
                 p1(k,j,i)=ga(k,j,1)
              end do
           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 j=1,lon2
                do k=1,lat2
                    gakr=ga(k,j,kr)
                    gaki=ga(k,j,ki)
                    ga(k,j,kr)=alv(k,kr,i,iv)*gakr-alv(k,ki,i,iv)*gaki+bekr*p2(k,j,i)
                    ga(k,j,ki)=alv(k,ki,i,iv)*gakr+alv(k,kr,i,iv)*gaki+beki*p2(k,j,i)
                    p1(k,j,i)=p1(k,j,i)+ga(k,j,kr)
                 end do
              enddo
           enddo
        enddo

!       backing filter:
        do i=nsig,1,-1
           do j=1,lon2
              do k=1,lat2
                 p1(k,j,i)=p1(k,j,i)+de(k,j,1)
                 de(k,j,1)=alv(k,1,i,iv)*(de(k,j,1)+be(1)*p2(k,j,i))
              end do
           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 j=1,lon2
                 do k=1,lat2
                    p1(k,j,i)=p1(k,j,i)+de(k,j,kr)
                    dekr=de(k,j,kr)+bekr*p2(k,j,i)
                    deki=de(k,j,ki)+beki*p2(k,j,i)
                    de(k,j,kr)=alv(k,kr,i,iv)*dekr-alv(k,ki,i,iv)*deki
                    de(k,j,ki)=alv(k,ki,i,iv)*dekr+alv(k,kr,i,iv)*deki
                 end do
              enddo
           enddo
        enddo

! Global branch, even degree     
     else

        ! advancing filter:
        do i=1,nsig
           !       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 j=1,lon2
                 do k=1,lat2
                    gakr=ga(k,j,kr)
                    gaki=ga(k,j,ki)

                    ga(k,j,kr)=alv(k,kr,i,iv)*gakr-alv(k,ki,i,iv)*gaki+bekr*p2(k,j,i)
                    ga(k,j,ki)=alv(k,ki,i,iv)*gakr+alv(k,kr,i,iv)*gaki+beki*p2(k,j,i)
                    p1(k,j,i)= p1(k,j,i)+ga(k,j,kr)
                 end do
              enddo
           enddo
        enddo
     
        !    backing filter:
        do i=nsig,1,-1
           !       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 j=1,lon2
                 do k=1,lat2
                    p1(k,j,i)=p1(k,j,i)+de(k,j,kr)
                    dekr=de(k,j,kr)+bekr*p2(k,j,i)
                    deki=de(k,j,ki)+beki*p2(k,j,i)
                    de(k,j,kr)=alv(k,kr,i,iv)*dekr-alv(k,ki,i,iv)*deki
                    de(k,j,ki)=alv(k,ki,i,iv)*dekr+alv(k,kr,i,iv)*deki
                 end do
              enddo
           enddo
        enddo
     endif
  end if
  return
end subroutine frfhvo

subroutine smoothzo(vx,samp,rate,iv,jx,dsv)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    smoothzo    initializes and renormalizes vertical smoothing coefs.
!   prgmmr: derber           org: np22                date: 2004-05-13
!
! abstract: initializes and renormalizes vertical smoothing coefficients
!           initializes dssv and alv
!
! program history log:
!   2004-05-13  derber, document
!   2004-11-30  treadon - add longitude dimension to variance array dssv
!   2010-03-01  zhu     - decide the location in alv and dsv based on anavinfo
!                       - add dsv in the interface, rm dssv in berror
!   2010-10-08  derber - optimize and clean up
!
!   input argument list:
!     vx       - vertical smoothing scales
!     samp     - parameter for smoothing        
!     rate     - parameter for smoothing       
!     iv       - location in alv and dssv for smoothing coefficients
!     jx       - latitude index
!
!   output argument list
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only:  zero
  use gridmod, only: nsig,lon2
  use berror, only: alv,ndeg
  implicit none
 
  integer(i_kind)             ,intent(in   ) :: jx,iv
  real(r_kind),dimension(nsig),intent(in   ) :: vx
  real(r_kind)                ,intent(in   ) :: samp
  real(r_kind),dimension(ndeg),intent(in   ) :: rate
  real(r_kind),dimension(lon2,nsig),intent(out):: dsv
 
  integer(i_kind) i,k,m
  real(r_kind),dimension(nsig):: dss
  real(r_kind),dimension(nsig,nsig):: p1
  real(r_kind),dimension(nsig,ndeg):: al

  call rfdparv(vx,rate,al,nsig,ndeg)
  do m=1,ndeg
     do k=1,nsig
        alv(jx,m,k,iv)=al(k,m)
     end do
  end do
  p1=zero
  do k=1,nsig
     dss(k)=sqrt(samp*vx(k))
     p1(k,k)=dss(k)
  end do

  call rfhvo(p1,nsig,nsig,al)
  
  call rfhvo(p1,nsig,nsig,al)

  do k=1,nsig
     do i=1,lon2
        dsv(i,k)=sqrt(dss(k)/p1(k,k))
     end do
  end do

  return
end subroutine smoothzo

subroutine rfhvo(p1,nc,n,al)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rfhvo   performs vertical smoothing for renormalization
!   prgmmr: derber           org: np22                date: 2004-05-13
!
! abstract: performs vertical smoothing of identity matrix for use with 
!           renormalization
!
! program history log:
!   2004-05-13  derber, document
!   2010-10-08  derber - optimize and clean up
!
!   input argument list:
!     p1       - input field to be smoothed (nc,n)
!     nc       - first array dimension for p1
!     n        - second array dimension for p1
!     al       - smoothing coefficients
!
!   output argument list
!     p1       - output field after vertical smoothing
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: zero
  use berror, only: be,ndeg
  implicit none
 
  integer(i_kind)               ,intent(in   ) :: n,nc
  real(r_kind),dimension(n,ndeg),intent(in   ) :: al
  real(r_kind),dimension(nc,n)  ,intent(inout) :: p1
 
  integer(i_kind) i,j,kr,ki
  real(r_kind) gaki,dekr,deki,gakr
  real(r_kind),dimension(nc,n):: p2
  real(r_kind),dimension(nc,ndeg):: ga,de


! Zero local work arrays.
  do j=1,ndeg
     do i=1,nc
        ga(i,j)=zero
        de(i,j)=zero
     end do
  end do

  do j=1,n
     do i=1,nc
        p2(i,j)=p1(i,j)
        p1(i,j)=zero
     end do
  end do

  if (mod(ndeg,2) == 1) then
!    advancing filter:
     do i=1,n
        do j=1,nc
           ga(j,1)=al(i,1)*ga(j,1)+be(1)*p2(j,i)
           p1(j,i)=p1(j,i)+ga(j,1)
        enddo
                           ! treat remaining complex roots:
        do kr=2,ndeg,2   ! <-- index of "real" components
           ki=kr+1            ! <-- index of "imag" components
           do j=1,nc
              gakr=ga(j,kr)
              gaki=ga(j,ki)
              ga(j,kr)=al(i,kr)*gakr-al(i,ki)*gaki+be(kr)*p2(j,i)
              ga(j,ki)=al(i,ki)*gakr+al(i,kr)*gaki+be(ki)*p2(j,i)
              p1(j,i)=p1(j,i)+ga(j,kr)
           enddo
        enddo
     enddo

!    backing filter:
     do i=n,1,-1
        do j=1,nc
           p1(j,i)=p1(j,i)+de(j,1)
           de(j,1)=al(i,1)*(de(j,1)+be(1)*p2(j,i))
        enddo
                           ! treat remaining complex roots:
        do kr=2,ndeg,2   ! <-- index of "real" components
           ki=kr+1            ! <-- index of "imag" components
           do j=1,nc
              p1(j,i)=p1(j,i)+de(j,kr)
              dekr=de(j,kr)+be(kr)*p2(j,i)
              deki=de(j,ki)+be(ki)*p2(j,i)
              de(j,kr)=al(i,kr)*dekr-al(i,ki)*deki
              de(j,ki)=al(i,ki)*dekr+al(i,kr)*deki
           enddo
        enddo
     enddo
     
  else
! advancing filter:
     do i=1,n
                           ! treat remaining complex roots:
        do kr=1,ndeg,2   ! <-- index of "real" components
           ki=kr+1            ! <-- index of "imag" components
           do j=1,nc
              gakr=ga(j,kr)
              gaki=ga(j,ki)
              ga(j,kr)=al(i,kr)*gakr-al(i,ki)*gaki+be(kr)*p2(j,i)
              ga(j,ki)=al(i,ki)*gakr+al(i,kr)*gaki+be(ki)*p2(j,i)
              p1(j,i)=p1(j,i)+ga(j,kr)
           enddo
        enddo
     enddo

! backing filter:
     do i=n,1,-1
        ! treat remaining complex roots:
        do kr=1,ndeg,2   ! <-- index of "real" components
           ki=kr+1            ! <-- index of "imag" components
           do j=1,nc
              p1(j,i)=p1(j,i)+de(j,kr)
              dekr=de(j,kr)+be(kr)*p2(j,i)
              deki=de(j,ki)+be(ki)*p2(j,i)
              de(j,kr)=al(i,kr)*dekr-al(i,ki)*deki
              de(j,ki)=al(i,ki)*dekr+al(i,kr)*deki
           enddo
        enddo
     enddo
     
  endif
  return
end subroutine rfhvo

subroutine smoothzo1(vx,samp,rate,dsv)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    smoothzo    initializes and renormalizes vertical smoothing
! coefs.
!   prgmmr: derber           org: np22                date: 2004-05-13
!
! abstract: initializes and renormalizes vertical smoothing coefficients
!           initializes dssv  (no alv)
!
! program history log:
!   2004-05-13  derber, document
!   2004-11-30  treadon - add longitude dimension to variance array dssv
!   2010-03-01  zhu     - decide the location in alv and dsv based on anavinfo
!                       - add dsv in the interface, rm dssv in berror
!   2010-10-08  derber - optimize and clean up
!
!   2017-02-09  CAPS(G. Zhao) - remove initialization of alv (it is defined on
!                               llmin:llmax for regional, might cause problem.)
!
!   input argument list:
!     vx       - vertical smoothing scales
!     samp     - parameter for smoothing
!     rate     - parameter for smoothing
!     iv       - location in alv and dssv for smoothing coefficients
!     jx       - latitude index
!
!   output argument list
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only:  zero
  use gridmod, only: nsig,lon2
  use berror, only: ndeg
  implicit none

  real(r_kind),dimension(nsig),intent(in   ) :: vx
  real(r_kind)                ,intent(in   ) :: samp
  real(r_kind),dimension(ndeg),intent(in   ) :: rate
  real(r_kind),dimension(lon2,nsig),intent(out):: dsv

  integer(i_kind) i,k
  real(r_kind),dimension(nsig):: dss
  real(r_kind),dimension(nsig,nsig):: p1
  real(r_kind),dimension(nsig,ndeg):: al

  call rfdparv(vx,rate,al,nsig,ndeg)
  p1=zero
  do k=1,nsig
     dss(k)=sqrt(samp*vx(k))
     p1(k,k)=dss(k)
  end do

  call rfhvo(p1,nsig,nsig,al)

  call rfhvo(p1,nsig,nsig,al)

  do k=1,nsig
     do i=1,lon2
        dsv(i,k)=sqrt(dss(k)/p1(k,k))
     end do
  end do

  return
end subroutine smoothzo1