subroutine getprs(ps,prs)
! subprogram:    getprs       get 3d pressure or 3d pressure deriv
!   prgmmr: kleist           org: np20                date: 2005-09-29
!
! abstract: calculate 3d pressure and its horizontal derivatives
!
! program history log:
!   2005-09-29  kleist
!   2006-04-12  treadon - replace sigi with bk5
!   2006-07-31  kleist  - analysis variable change from ln(ps) to ps
!   2007-05-08  kleist  - add generalized vert coord and derivative call
!   2007-07-26  cucurull- compute 3d pressure and derivatives in different subroutines
!                       - remove gues_tv from argument list; clean up code
!   2008-06-04  safford - rm unused uses
!   2008-09-05  lueken  - merged ed's changes into q1fy09 code
!   2010-09-15  pagowski  - added cmaq
!   2013-10-19  todling - metguess now holds background
!   2017-03-23  Hu      - add code to use hybrid vertical coodinate in WRF MASS core
!   2018-02-15  wu      - add code for fv3_regional
!
!   input argument list:
!     ps       - surface pressure
!
!   output argument list:
!     prs        - 3d pressure
!
! attributes:
!   language:  f90
!   machine:   ibm/RS6000 SP
!
!$$$ end documentation block

  use kinds,only: r_kind,i_kind
  use constants,only: zero,half,one_tenth,rd_over_cp,one
  use gridmod,only: nsig,lat2,lon2,ak5,bk5,ck5,tref5,idvc5
  use gridmod,only: wrf_nmm_regional,nems_nmmb_regional,eta1_ll,eta2_ll,pdtop_ll,pt_ll,&
       regional,wrf_mass_regional,cmaq_regional,twodvar_regional,fv3_regional
  use guess_grids, only: ntguessig
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2)       ,intent(in   ) :: ps
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(  out) :: prs

! Declare local variables
  real(r_kind) kapr,trk
  real(r_kind),dimension(:,:,:),pointer::ges_tv_it=>NULL()
  integer(i_kind) i,j,k,k2,it,istatus

! Declare local parameter
  real(r_kind),parameter:: ten = 10.0_r_kind

! prs=zero 
  it=ntguessig

  if (regional) then
     if(wrf_nmm_regional.or.nems_nmmb_regional.or.&
          cmaq_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=one_tenth* &
                      (eta1_ll(k)*pdtop_ll + &
                      eta2_ll(k)*(ten*ps(i,j)-pdtop_ll-pt_ll) + &
                      pt_ll)
              end do
           end do
        end do
     elseif (fv3_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=eta1_ll(k)+ eta2_ll(k)*ps(i,j)
              end do
           end do
        end do

     elseif (twodvar_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=one_tenth*(eta1_ll(k)*(ten*ps(i,j)-pt_ll) + pt_ll)
              end do
           end do
        end do
     elseif (wrf_mass_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=one_tenth*(eta1_ll(k)*(ten*ps(i,j)-pt_ll) + &
                                       eta2_ll(k) + pt_ll)
              end do
           end do
        end do
     endif
  else
     k=1
     k2=nsig+1
     do j=1,lon2
        do i=1,lat2
           prs(i,j,k)=ps(i,j)
           prs(i,j,k2)=zero
        end do
     end do
     if (idvc5 /= 3) then
        do k=2,nsig
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=ak5(k)+bk5(k)*ps(i,j)
              end do
           end do
        end do
     else
        kapr=one/rd_over_cp
        do k=2,nsig
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=ak5(k)+(bk5(k)*ps(i,j))
              end do
           end do
        end do
        call gsi_bundlegetpointer(gsi_metguess_bundle(it),'tv',ges_tv_it,istatus)
        if(istatus==0) then
           do k=2,nsig
              do j=1,lon2
                 do i=1,lat2
                    trk=(half*(ges_tv_it(i,j,k-1)+ges_tv_it(i,j,k))/tref5(k))**kapr
                    prs(i,j,k)=prs(i,j,k)+(ck5(k)*trk)
                 end do
              end do
           end do
        end if
     end if
  end if

  return
end subroutine getprs


subroutine getprs_horiz(ps_x,ps_y,prs,prs_x,prs_y)
!$$$  subprogram docuentation block
!                .     .    .                     .
! subprogram:    getprs_horiz
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-06-04  safford - complete documentation block, rm unused var k2
!   2008-09-05  lueken  - merged ed's changes into q1fy09 code
!   2012-06-12  parrish - replace sub2grid2/grid2sub2 with general_sub2grid/general_grid2sub
!
!   input argument list:
!     prs      - 3d pressure
!     ps_x     - dps/dx
!     ps_y     - dps/dy
!
!   output argument list:
!     prs_x      - dp/dx
!     prs_y      - dp/dy
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$
  
  use kinds,only: r_kind,i_kind
  use constants,only: zero
  use gridmod,only: nsig,lat2,lon2
  use gridmod,only: regional,wrf_nmm_regional,nems_nmmb_regional,eta2_ll,&
       cmaq_regional,fv3_regional
  use compact_diffs, only: compact_dlat,compact_dlon
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  use general_commvars_mod, only: s2g2
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2)       ,intent(in   ) :: ps_x,ps_y
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(in   ) :: prs
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(  out) :: prs_x,prs_y

! Declare local variables
  integer(i_kind) i,j,k
  real(r_kind),allocatable,dimension(:,:,:,:):: hwork,hwork_g
  real(r_kind),dimension(1,lat2,lon2,nsig+1):: p4

  allocate(hwork(s2g2%inner_vars,s2g2%nlat,s2g2%nlon,s2g2%kbegin_loc:s2g2%kend_alloc))
  allocate(hwork_g(s2g2%inner_vars,s2g2%nlat,s2g2%nlon,s2g2%kbegin_loc:s2g2%kend_alloc))

  if(regional)then
     if(wrf_nmm_regional.or.nems_nmmb_regional.or.cmaq_regional&
         .or. fv3_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 prs_x(i,j,k)=eta2_ll(k)*ps_x(i,j)
                 prs_y(i,j,k)=eta2_ll(k)*ps_y(i,j)
              end do
           end do
        end do
     else
        prs_x=zero ; prs_y=zero
     end if
  else
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              p4(1,i,j,k)=prs(i,j,k)
           end do
        end do
     end do
     call general_sub2grid(s2g2,p4,hwork)
     do k=s2g2%kbegin_loc,s2g2%kend_loc
        call compact_dlon(hwork(1,:,:,k),hwork_g(1,:,:,k),.false.)
     end do
     call general_grid2sub(s2g2,hwork_g,p4)
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              prs_x(i,j,k)=p4(1,i,j,k)
           end do
        end do
     end do
     do k=s2g2%kbegin_loc,s2g2%kend_loc
        call compact_dlat(hwork(1,:,:,k),hwork_g(1,:,:,k),.false.)
     end do
     call general_grid2sub(s2g2,hwork_g,p4)
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              prs_y(i,j,k)=p4(1,i,j,k)
           end do
        end do
     end do
  end if

!  clean work space
  deallocate(hwork,hwork_g)

  return
end subroutine getprs_horiz


subroutine getprs_tl(ps,t,prs)
!$$$ subprogram documentation block
!               .      .    .                       .
! subprogram:    getprs_tl    TLM of getprs
!   prgmmr: kleist           org: np20                date: 2005-09-29
!
! abstract: TLM of routine that gets 3d pressure and its derivatives
!
! program history log:
!   2005-09-29  kleist
!   2006-04-12  treadon - replace sigi with bk5
!   2006-07-31  kleist  - analysis variable change from ln(ps) to ps
!   2007-05-08  kleist  - add generalized vert coord and derivative call
!   2007-07-26  cucurull- compute 3d pressure and derivatives in different subroutines
!                       - remove gues_tv from argument list; clean up code;
!                       - fix buf for t dimension
!   2008-06-04  safford - complete doc block, rm unused uses
!   2008-09-05  lueken  - merged ed's changes into q1fy09 code
!   2013-10-19  todling - metguess now holds background
!
!   input argument list:
!     ps       - surface pressure
!
!   output argument list:
!     prs        - 3d pressure
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block
  
  use kinds,only: r_kind,i_kind
  use constants,only: zero,one,rd_over_cp,half
  use gridmod,only: nsig,lat2,lon2,bk5,ck5,idvc5,tref5
  use gridmod,only: wrf_nmm_regional,nems_nmmb_regional,eta2_ll,eta1_ll,regional,wrf_mass_regional,cmaq_regional,&
       twodvar_regional,fv3_regional
  use guess_grids, only: ntguessig
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2)       ,intent(in   ) :: ps
  real(r_kind),dimension(lat2,lon2,nsig)  ,intent(in   ) :: t
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(  out) :: prs

! Declare local variables
  real(r_kind) kapr,kaprm1,trk,tc1,t9trm
  real(r_kind),dimension(:,:,:),pointer::ges_tv_it=>NULL()
  integer(i_kind) i,j,k,k2,it,istatus

  if (regional) then
     if(wrf_nmm_regional.or.nems_nmmb_regional.or.&
        fv3_regional .or.  cmaq_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=eta2_ll(k)*ps(i,j)
              end do
           end do
        end do
     elseif (wrf_mass_regional .or. twodvar_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=eta1_ll(k)*ps(i,j)
              end do
           end do
        end do
     endif
  else
     k=1
     k2=nsig+1
     do j=1,lon2
        do i=1,lat2
           prs(i,j,k)=ps(i,j)
           prs(i,j,k2)=zero
        end do
     end do
     if (idvc5 /= 3) then
        do k=2,nsig
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=bk5(k)*ps(i,j)
              end do
           end do
        end do
     else
        kapr=one/rd_over_cp
        kaprm1=kapr-one
        it=ntguessig
        do k=2,nsig
           do j=1,lon2
              do i=1,lat2
                 prs(i,j,k)=bk5(k)*ps(i,j)
              end do
           end do
        end do
        call gsi_bundlegetpointer(gsi_metguess_bundle(it),'tv',ges_tv_it,istatus)
        if(istatus==0) then
           do k=2,nsig
              do j=1,lon2
                 do i=1,lat2
                    t9trm=half*(ges_tv_it(i,j,k-1)+ges_tv_it(i,j,k))/tref5(k)
                    tc1=half/tref5(k)
                    trk=kapr*tc1*(t(i,j,k-1)+t(i,j,k))*(t9trm**kaprm1)
                    prs(i,j,k)=prs(i,j,k) + ck5(k)*trk
                 end do
              end do
           end do
        end if
     end if
  end if

  return
end subroutine getprs_tl


subroutine getprs_horiz_tl(ps_x,ps_y,prs,prs_x,prs_y)
!$$$ subprogram documentation block
!               .      .    .                     .
! subprogram:   getprs_horiz_tl
!
!  prgrmmr:
!
! abstract:
!
! program history log:
!   2008-06-04  safford - complete doc block
!   2008-09-05  lueken  - merged ed's changes into q1fy09 code
!   2010-05-23  todling - unwired location of ps in control array
!   2012-06-12  parrish - replace sub2grid2/grid2sub2 with general_sub2grid/general_grid2sub
!
!   input argument list:
!     prs      - 3d pressure
!     ps_x     - dps/dx
!     ps_y     - dps/dy
!
!   output argument list:
!     prs_x      - dp/dx
!     prs_y      - dp/dy
!
! attributes:
!   language:  f90
!   machine    ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds,only: r_kind,i_kind
  use constants,only: zero
  use gridmod,only: nsig,lat2,lon2
  use gridmod,only: regional,wrf_nmm_regional,nems_nmmb_regional,eta2_ll,&
       cmaq_regional,fv3_regional
  use compact_diffs, only: compact_dlat,compact_dlon
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  use general_commvars_mod, only: s2g2
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2)       ,intent(in   ) :: ps_x,ps_y
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(in   ) :: prs
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(  out) :: prs_x,prs_y

! Declare local variables
  integer(i_kind) i,j,k
  real(r_kind),allocatable,dimension(:,:,:,:):: hwork,hwork_g
  real(r_kind),dimension(1,lat2,lon2,nsig+1):: p4

  allocate(hwork(s2g2%inner_vars,s2g2%nlat,s2g2%nlon,s2g2%kbegin_loc:s2g2%kend_alloc))
  allocate(hwork_g(s2g2%inner_vars,s2g2%nlat,s2g2%nlon,s2g2%kbegin_loc:s2g2%kend_alloc))

  if(regional)then
     if(wrf_nmm_regional.or.nems_nmmb_regional.or.cmaq_regional&
         .or. fv3_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 prs_x(i,j,k)=eta2_ll(k)*ps_x(i,j)
                 prs_y(i,j,k)=eta2_ll(k)*ps_y(i,j)
              end do
           end do
        end do
     else
        prs_x=zero
        prs_y=zero
     end if
  else
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              p4(1,i,j,k)=prs(i,j,k)
           end do
        end do
     end do
     call general_sub2grid(s2g2,p4,hwork)
     do k=s2g2%kbegin_loc,s2g2%kend_loc
        call compact_dlon(hwork(1,:,:,k),hwork_g(1,:,:,k),.false.)
     end do
     call general_grid2sub(s2g2,hwork_g,p4)
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              prs_x(i,j,k)=p4(1,i,j,k)
           end do
        end do
     end do
     do k=s2g2%kbegin_loc,s2g2%kend_loc
        call compact_dlat(hwork(1,:,:,k),hwork_g(1,:,:,k),.false.)
     end do
     call general_grid2sub(s2g2,hwork_g,p4)
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              prs_y(i,j,k)=p4(1,i,j,k)
           end do
        end do
     end do
  end if

!  clean work space
  deallocate(hwork,hwork_g)

  return
end subroutine getprs_horiz_tl


subroutine getprs_ad(ps,t,prs)
!$$$ subprogram documentation block
!               .      .    .                    .
! subprogram:    getprs_ad    adjoint of getprs_tl
!   prgmmr: kleist           org: np20                date: 2005-09-29
!
! abstract: adjoint of linear routine that gets 3d pressure and derivs
!
! program history log:
!   2005-09-29  kleist
!   2006-04-12  treadon - replace sigi with bk5
!   2006-07-31  kleist  - analysis variable change from ln(ps) to ps
!   2007-05-08  kleist  - add generalized vert coord and derivative call
!   2007-07-26  cucurull- compute 3d pressure and derivatives in different subroutines
!                       - remove gues_tv from argument list; clean up code
!   2008-06-04  safford - complete doc block, rm unused uses
!   2008-09-05  lueken  - merged ed's changes into q1fy09 code
!   2013-10-19  todling - metguess now holds background
!
!   input argument list:
!     prs        - 3d pressure
!
!   output argument list:
!     ps       - log surface pressure
!
!  notes:
!     Adjoint check performed and verified on 2005-08-29 by d. kleist
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block
  
  use kinds,only: r_kind,i_kind
  use gridmod,only: nsig,lat2,lon2,bk5,ck5,tref5,idvc5
  use gridmod,only: wrf_nmm_regional,nems_nmmb_regional,eta2_ll,regional,wrf_mass_regional,cmaq_regional,eta1_ll,&
       twodvar_regional,fv3_regional
  use guess_grids, only: ntguessig 
  use constants,only: zero,half,one,rd_over_cp
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer

  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(inout) :: prs
  real(r_kind),dimension(lat2,lon2,nsig)  ,intent(inout) :: t
  real(r_kind),dimension(lat2,lon2)       ,intent(inout) :: ps

! Declare local variables
  real(r_kind) kapr,kaprm1,trk,tc1,t9trm
  real(r_kind),dimension(:,:,:),pointer::ges_tv_it=>NULL()
  integer(i_kind) i,j,k,it,istatus


  if (regional) then
     if(wrf_nmm_regional.or.nems_nmmb_regional.or.&
          cmaq_regional .or. fv3_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 ps(i,j) = ps(i,j) + eta2_ll(k)*prs(i,j,k)
              end do
           end do
        end do
     elseif (wrf_mass_regional .or. twodvar_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 ps(i,j) = ps(i,j) + eta1_ll(k)*prs(i,j,k)
              end do
           end do
        end do
     endif
  else
     if (idvc5 /= 3) then
        do k=2,nsig
           do j=1,lon2
              do i=1,lat2
                 ps(i,j) = ps(i,j) + bk5(k)*prs(i,j,k)
              end do
           end do
        end do
     else
        kapr=one/rd_over_cp
        kaprm1=kapr-one
        it=ntguessig
        call gsi_bundlegetpointer(gsi_metguess_bundle(it),'tv',ges_tv_it,istatus)
        if(istatus==0) then
           do k=2,nsig
              do j=1,lon2
                 do i=1,lat2
                    t9trm=half*(ges_tv_it(i,j,k-1)+ges_tv_it(i,j,k))/tref5(k)
                    tc1=half/tref5(k)
                    ps(i,j) = ps(i,j) + bk5(k)*prs(i,j,k)
                    trk = ck5(k)*prs(i,j,k)
                    t(i,j,k-1) = t(i,j,k-1) + kapr*tc1*trk*(t9trm**kaprm1)
                    t(i,j,k  ) = t(i,j,k  ) + kapr*tc1*trk*(t9trm**kaprm1)
                 end do
              end do
           end do
        else
           do k=2,nsig
              do j=1,lon2
                 do i=1,lat2
                    ps(i,j) = ps(i,j) + bk5(k)*prs(i,j,k)
                 end do
              end do
           end do
        end if
     end if
     k=1
     do j=1,lon2
        do i=1,lat2
           ps(i,j)=ps(i,j) + prs(i,j,k)
        end do
     end do
  end if

  do k=1,nsig+1
     do j=1,lon2
        do i=1,lat2
           prs(i,j,k)=zero
        end do
     end do
  end do 
 
  return
end subroutine getprs_ad


subroutine getprs_horiz_ad(ps_x,ps_y,prs,prs_x,prs_y)
!$$$ subprogram documentation block
!               .      .    .              .
! subprogram:  getprs_horiz_ad
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-06-04  safford - complete doc block
!   2008-09-05  lueken  - merged ed's changes into q1fy09
!   2010-05-23  todling - unwired location of ps in control array
!   2012-06-12  parrish - replace sub2grid2/grid2sub2 with general_sub2grid/general_grid2sub
!
!   input argument list:
!     prs_x      - dp/dx
!     prs_y      - dp/dy
!
!   output argument list:
!     prs      - 3d pressure
!     ps_x     - d(ln(ps))/dx
!     ps_y     - d(ln(ps))/dy
!
!  notes:
!     Adjoint check performed and verified on 2005-08-29 by d. kleist
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds,only: r_kind,i_kind
  use constants,only: zero
  use gridmod,only: nsig,lat2,lon2
  use gridmod,only: regional,wrf_nmm_regional,nems_nmmb_regional,eta2_ll,&
       cmaq_regional,fv3_regional
  use compact_diffs, only: tcompact_dlat,tcompact_dlon
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  use general_commvars_mod, only: s2g2

  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(in   ) :: prs_x,prs_y
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(inout) :: prs
  real(r_kind),dimension(lat2,lon2)       ,intent(inout) :: ps_x,ps_y

! Declare local variables
  integer(i_kind) i,j,k
  real(r_kind),allocatable,dimension(:,:,:,:):: hwork,hwork_g
  real(r_kind),dimension(1,lat2,lon2,nsig+1):: p4

  allocate(hwork(s2g2%inner_vars,s2g2%nlat,s2g2%nlon,s2g2%kbegin_loc:s2g2%kend_alloc))
  allocate(hwork_g(s2g2%inner_vars,s2g2%nlat,s2g2%nlon,s2g2%kbegin_loc:s2g2%kend_alloc))

! Adjoint of horizontal derivatives
  if (regional) then
     if(wrf_nmm_regional.or.nems_nmmb_regional.or.cmaq_regional&
         .or. fv3_regional) then
        do k=1,nsig+1
           do j=1,lon2
              do i=1,lat2
                 ps_y(i,j) = ps_y(i,j) + eta2_ll(k)*prs_y(i,j,k)
                 ps_x(i,j) = ps_x(i,j) + eta2_ll(k)*prs_x(i,j,k)
              end do
           end do
        end do
     end if
  else
     hwork=zero
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              p4(1,i,j,k)=prs_x(i,j,k)
           end do
        end do
     end do
     call general_sub2grid(s2g2,p4,hwork_g)
     do k=s2g2%kbegin_loc,s2g2%kend_loc
        call tcompact_dlon(hwork(1,:,:,k),hwork_g(1,:,:,k),.false.)
     end do
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              p4(1,i,j,k)=prs_y(i,j,k)
           end do
        end do
     end do
     call general_sub2grid(s2g2,p4,hwork_g)
     do k=s2g2%kbegin_loc,s2g2%kend_loc
        call tcompact_dlat(hwork(1,:,:,k),hwork_g(1,:,:,k),.false.)
     end do
     call general_grid2sub(s2g2,hwork,p4)
     do k=1,nsig+1
        do j=1,lon2
           do i=1,lat2
              prs(i,j,k)=prs(i,j,k)+p4(1,i,j,k)
           end do
        end do
     end do
  end if

!  clean work space
  deallocate(hwork,hwork_g)

  return
end subroutine getprs_horiz_ad