module patch2grid_mod
!$$$   module documentation block
!                .      .    .                                       .
! module:    patch2grid_mod
!   prgmmr: sato             org: np23                date: 2007-09-20
!
! abstract: contains the stuff to convert the full latlon grid data to
!           three filter patch spaces.
!           The three patches are north/south polar, and zonal patch.
!
! program history log:
!   2007-09-20  sato
!
! subroutines included:
!   sub setup_patch2grid
!   sub setup_blend
!   sub destroy_patch2grid
!   sub grid2patch
!   sub patch2grid
!   sub tpatch2grid
!   sub vpatch2grid
!
! Variable Definitions:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_kind,i_kind
  use constants,only: zero,one
  use anberror, only: pf2aP1,pf2aP2,pf2aP3,&
                      nx,ny,mr,nr,nf
  use gridmod,  only: nlat,nlon

  implicit none

! set default to private
  private
! set subroutines to public
  public :: setup_patch2grid
  public :: setup_blend
  public :: destroy_patch2grid
  public :: grid2patch
  public :: patch2grid
  public :: tpatch2grid
  public :: vpatch2grid

  integer(i_kind):: ndx,ndy,ndx2,nmix,nrmxb,nmixp,nymx,nlatxb

  real(r_kind),allocatable,dimension(:,:):: p1all
  real(r_kind),allocatable,dimension(:,:):: p2all,p3all
  real(r_kind),allocatable,dimension(:,:):: p2pol,p3pol

  real(r_kind),allocatable,dimension(:):: bl,bl2

  contains

subroutine setup_patch2grid
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    setup_patch2grid
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-22  lueken - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
  implicit none

  integer(i_kind):: ier

  ndx=(nx-nlon)/2
  ndy=(nlat-ny)/2
  ndx2=2*ndx
  nmix=nr+1+(ny-nlat)/2
  nrmxb=ndy-1
  nmixp=nmix+1
  nymx=ny-nmix
  nlatxb=nlat-nrmxb

  allocate(p1all(ny,nx))
  allocate(p2all(nlon+1,mr:nr))
  allocate(p3all(nlon+1,mr:nr))
  allocate(p2pol(nf*2+1,nf*2+1))
  allocate(p3pol(nf*2+1,nf*2+1),stat=ier)

  if( ier /= 0 ) then
     write(6,*) 'setup_patch2grid: could not allocate memories'
     call stop2(99)
  end if

!   initialize blend array (bl,bl2)
  allocate( bl(nx-nlon), bl2(nr+1+(ny-nlat)/2) )
  call setup_blend

end subroutine

subroutine setup_blend
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    setup_blend
!   prgmmr:
!
! abstract: setup blend array bl & bl2
!
! program history log:
!   2009-09-22  lueken - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use blendmod, only: blend
  implicit none

  integer(i_kind),dimension(0:40):: iblend
  integer(i_kind):: mm,nolp,nbuf,nmixbl,ndxbl
  real(r_kind):: dxx,x,y

  integer(i_kind):: i,j,k

! Setup blending
  mm=4
  call blend(mm,iblend)

  nolp=nr+1+(ny-nlat)/2
  nbuf=0
  nmixbl=nolp-nbuf*2
  dxx=one/(nmixbl+1)
  bl2=zero
  k=0
  do i=1,nmixbl
     k=k+1
     x=i*dxx
     y=iblend(mm)
     do j=mm-1,0,-1
        y=x*y+iblend(j)
     end do
     y=y*x**(mm+1)
     bl2(k)=one-y
  end do
  do k=1,nmixbl
     bl2(k)=sqrt(bl2(k))
  end do

  nmixbl=(nx-nlon)
  dxx=one/(nmixbl+1)
  ndxbl=(nx-nlon)
  bl=zero
  k=ndxbl-nmixbl
  do i=1,nmixbl
     k=k+1
     x=i*dxx
     y=iblend(mm)
     do j=mm-1,0,-1
        y=x*y+iblend(j)
     end do
     y=y*x**(mm+1)
     bl(k)=one-y
  enddo
  do k=1,nmixbl
     bl(k)=sqrt(bl(k))
  end do

end subroutine setup_blend

subroutine destroy_patch2grid
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroy_patch2grid
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-22  lueken - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
  implicit none

  deallocate(p1all)
  deallocate(p2all,p3all)
  deallocate(p2pol,p3pol)
  deallocate(bl,bl2)

end subroutine

!=======================================================================
!=======================================================================
subroutine grid2patch(grid_wrk,hflt_all,hflt_p2,hflt_p3)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grid2patch
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-22  lueken - added subprogram doc block
!
!   input argument list:
!    grid_wrk
!
!   output argument list:
!    hflt_all
!    hflt_p2
!    hflt_p3
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use smooth_polcarf, only: smooth_caspol
  use fgrid2agrid_mod, only: agrid2fgrid

  implicit none

  real(r_kind),dimension(nlat,nlon),intent(in   ) :: grid_wrk

  real(r_kind)                     ,intent(  out) :: hflt_all(pf2aP1%nlatf,pf2aP1%nlonf)
  real(r_kind)                     ,intent(  out) :: hflt_p2 (pf2aP2%nlatf ,pf2aP2%nlonf )
  real(r_kind)                     ,intent(  out) :: hflt_p3 (pf2aP3%nlatf ,pf2aP3%nlonf )

  integer(i_kind):: i,i1,i2,j,j1

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

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

! Handle polar patches
  do i=1,nlon
     do j=mr,nrmxb+nmix
        p2all(i,j)=grid_wrk(nlat-j,i)
        p3all(i,j)=grid_wrk(j+1,i)
     end do
  end do

  call agrid2fgrid(pf2aP1,p1all,hflt_all) !analysis to filter grid_wrk

  call smooth_caspol(p2pol,p2all)
  call agrid2fgrid(pf2aP2 ,p2pol,hflt_p2)  !analysis to filter grid_wrk

  call smooth_caspol(p3pol,p3all)
  call agrid2fgrid(pf2aP3 ,p3pol,hflt_p3)  !analysis to filter grid_wrk

end subroutine grid2patch
!=======================================================================
!=======================================================================
subroutine tpatch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    tpatch2grid
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-22  lueken - added subprogram doc block
!
!   input argument list:
!    grid_wrk
!
!   output argument list:
!    hflt_all
!    hflt_p2
!    hflt_p3
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use smooth_polcarf, only: smooth_polcasa
  use fgrid2agrid_mod, only: tfgrid2agrid

  implicit none

  real(r_kind),dimension(nlat,nlon),intent(in   ) :: grid_wrk

  real(r_kind)                     ,intent(  out) :: hflt_all(pf2aP1%nlatf,pf2aP1%nlonf)
  real(r_kind)                     ,intent(  out) :: hflt_p2 (pf2aP2%nlatf ,pf2aP2%nlonf )
  real(r_kind)                     ,intent(  out) :: hflt_p3 (pf2aP3%nlatf ,pf2aP3%nlonf )

  integer(i_kind):: i,i1,i2,j,j1

! Extract central patch (band) from full grid_wrk (work --> p1)
! Blending zones
  p1all=zero

  do i=1,ndx
     i1=i-ndx+nlon
     i2=nx-ndx+i
     do j=1,ny
        j1=j+ndy
        p1all(j,i) =grid_wrk(j1,i1)      ! left (west) blending zone
        p1all(j,i2)=grid_wrk(j1,i)       ! right (east) blending zone
     end do
  end do

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

! 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
     end do
  end do

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

! Handle polar patches
  p2all=zero
  p3all=zero

  do j=mr,nrmxb+nmix
     do i=1,nlon
        p2all(i,j)=grid_wrk(nlat-j,i)
        p3all(i,j)=grid_wrk(j+1,i)
     end do
  end do

! Apply blending coefficients
  do j=nrmxb+1,nrmxb+nmix
     j1=j-nrmxb
     do i=1,nlon
        p2all(i,j)=p2all(i,j)*bl2(j1)
        p3all(i,j)=p3all(i,j)*bl2(j1)
     end do
  end do

  call tfgrid2agrid(pf2aP1,p1all,hflt_all) !analysis to filter grid_wrk

  call smooth_polcasa(p2pol,p2all)
  call tfgrid2agrid(pf2aP2,p2pol,hflt_p2) !analysis to filter grid_wrk

  call smooth_polcasa(p3pol,p3all)
  call tfgrid2agrid(pf2aP3,p3pol,hflt_p3) !analysis to filter grid_wrk

end subroutine tpatch2grid
!=======================================================================
!=======================================================================
subroutine patch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    patch2grid
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-22  lueken - added subprogram doc block
!
!   input argument list:
!    hflt_all
!    hflt_p2
!    hflt_p3
!
!   output argument list:
!    grid_wrk
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use smooth_polcarf, only: smooth_polcas
  use fgrid2agrid_mod, only: fgrid2agrid

  implicit none

  real(r_kind),dimension(nlat,nlon),intent(  out) :: grid_wrk

  real(r_kind)                     ,intent(in   ) :: hflt_all(pf2aP1%nlatf,pf2aP1%nlonf)
  real(r_kind)                     ,intent(in   ) :: hflt_p2 (pf2aP2%nlatf ,pf2aP2%nlonf )
  real(r_kind)                     ,intent(in   ) :: hflt_p3 (pf2aP3%nlatf ,pf2aP3%nlonf )

  integer(i_kind):: i,i1,i2,j,j1

  call fgrid2agrid(pf2aP1,hflt_all,p1all) !analysis to filter grid_wrk

  call fgrid2agrid(pf2aP2,hflt_p2,p2pol) !analysis to filter grid_wrk
  call smooth_polcas(p2pol,p2all)

  call fgrid2agrid(pf2aP3,hflt_p3,p3pol) !analysis to filter grid_wrk
  call smooth_polcas(p3pol,p3all)

! zero output array
  grid_wrk=zero

! 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
     end do
  end do

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

! 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
        grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)  ! left (west) blending zone
        grid_wrk(j1,i )=grid_wrk(j1,i )+p1all(j,i2) ! right (east) blending zone
     end do
  end do

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

! Adjoint of North pole patch(p2) -- blending and transfer to grid
! Adjoint of South pole patch(p3) -- blending and transfer to grid

  do j=nlatxb-nmix,nlatxb-1
! Adjoint of blending
     do i=1,nlon
        p2all(i,nlat-j)=p2all(i,nlat-j)*bl2(nlatxb-j)
     end do
  end do

  do j=nrmxb+1,nrmxb+nmix
!    Adjoint of blending
     do i=1,nlon
        p3all(i,j)=p3all(i,j)*bl2(j-nrmxb)
     end do
  end do

  do i=1,nlon
!    Adjoint of transfer
     do j=mr,nrmxb+nmix
        grid_wrk(j+1,i)=grid_wrk(j+1,i)+p3all(i,j)
     end do
     do j=nlatxb-nmix,nlat-mr
        grid_wrk(j  ,i)=grid_wrk(j  ,i)+p2all(i,nlat-j)
     end do
  end do

end subroutine patch2grid

!=======================================================================
!=======================================================================
subroutine vpatch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    vpatch2grid
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-22  lueken - added subprogram doc block
!
!   input argument list:
!    hflt_all
!    hflt_p2
!    hflt_p3
!
!   output argument list:
!    grid_wrk
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use smooth_polcarf, only: smooth_polcasv
  use fgrid2agrid_mod, only: fgrid2agrid

  implicit none

  real(r_kind),dimension(nlat,nlon),intent(  out) :: grid_wrk

  real(r_kind)                     ,intent(in   ) :: hflt_all(pf2aP1%nlatf,pf2aP1%nlonf)
  real(r_kind)                     ,intent(in   ) :: hflt_p2 (pf2aP2%nlatf ,pf2aP2%nlonf )
  real(r_kind)                     ,intent(in   ) :: hflt_p3 (pf2aP3%nlatf ,pf2aP3%nlonf )

  real(r_kind),allocatable,dimension(:,:) :: wgt_wrk
  real(r_kind),allocatable,dimension(:,:) :: p1wgt
  real(r_kind),allocatable,dimension(:,:) :: p2wgt,p3wgt

  integer(i_kind):: i,i1,i2,j,j1

  allocate(wgt_wrk(nlat,nlon))

  call fgrid2agrid(pf2aP1,hflt_all,p1all) !analysis to filter grid_wrk

  call fgrid2agrid(pf2aP2,hflt_p2,p2pol) !analysis to filter grid_wrk
  call smooth_polcasv(p2pol,p2all)

  call fgrid2agrid(pf2aP3,hflt_p3,p3pol) !analysis to filter grid_wrk
  call smooth_polcasv(p3pol,p3all)

  wgt_wrk=zero

! Equatorial patch
! Adjoint of central patch blending on left/right sides of patch
  allocate(p1wgt(ny,nx)); p1wgt=one

  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
        p1wgt(j,i) =p1wgt(j,i) *bl(i1)
        p1wgt(j,i2)=p1wgt(j,i2)*bl(i)
     end do
  end do

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

! zero output array
  grid_wrk=zero
  wgt_wrk=zero

! 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
        grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)  ! left  (west) blending zone
        grid_wrk(j1,i )=grid_wrk(j1,i )+p1all(j,i2) ! right (east) blending zone
        wgt_wrk (j1,i1)=wgt_wrk (j1,i1)+p1wgt(j,i)  ! left  (west) blending zone
        wgt_wrk (j1,i )=wgt_wrk (j1,i )+p1wgt(j,i2) ! right (east) blending zone
     end do
  end do

! Middle zone (no blending)
  do i=ndx+1,nx-ndx
     i1=i-ndx
     do j=1,ny
        j1=j+ndy
        grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)
        wgt_wrk (j1,i1)=wgt_wrk (j1,i1)+p1wgt(j,i)
     end do
  end do

  deallocate(p1wgt)

! Adjoint of North pole patch(p2) -- blending and transfer to grid
! Adjoint of South pole patch(p3) -- blending and transfer to grid
  allocate(p2wgt(nlon+1,mr:nr)); p2wgt=one
  allocate(p3wgt(nlon+1,mr:nr)); p3wgt=one

  do j=nlatxb-nmix,nlatxb-1
! Adjoint of blending
     do i=1,nlon
        p2all(i,nlat-j)=p2all(i,nlat-j)*bl2(nlatxb-j)
        p2wgt(i,nlat-j)=p2wgt(i,nlat-j)*bl2(nlatxb-j)
     end do
  end do

  do j=nrmxb+1,nrmxb+nmix
!    Adjoint of blending
     do i=1,nlon
        p3all(i,j)=p3all(i,j)*bl2(j-nrmxb)
        p3wgt(i,j)=p3wgt(i,j)*bl2(j-nrmxb)
     end do
  end do

  do i=1,nlon
!    Adjoint of transfer
     do j=mr,nrmxb+nmix
        grid_wrk(j+1,i)=grid_wrk(j+1,i)+p3all(i,j)
        wgt_wrk (j+1,i)=wgt_wrk (j+1,i)+p3wgt(i,j)
     end do
     do j=nlatxb-nmix,nlat-mr
        grid_wrk(j  ,i)=grid_wrk(j  ,i)+p2all(i,nlat-j)
        wgt_wrk (j  ,i)=wgt_wrk (j  ,i)+p2wgt(i,nlat-j)
     end do
  end do

  deallocate(p2wgt,p3wgt)

  where(wgt_wrk>zero) grid_wrk=grid_wrk/wgt_wrk

  deallocate(wgt_wrk)

end subroutine vpatch2grid

end module patch2grid_mod