!############################# Change Log ##################################
! 5.0.0
!
!###########################################################################
!  Copyright (C)  1990, 1995, 1999, 2000, 2003 - All Rights Reserved
!  Regional Atmospheric Modeling System - RAMS
!###########################################################################
module ModAcoust
contains



  subroutine coefz(mzp,mxp,myp,ia,iz,ja,jz,  &
       acoc,acof,acog,dn0,pi0,th0,rtgt,a1da2,amoe,amof,acoaa)

    use mem_grid, only: impl, dts, sspct, dzm, dzt, nzp, nz
    use mem_scratch, only: vctr11, vctr12
    use rconstants, only: rgas, cv

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(out) :: acoc(mzp,mxp,myp)
    real, intent(out) :: acof(mzp,mxp,myp)
    real, intent(out) :: acog(mzp,mxp,myp)
    real, intent(in) :: dn0(mzp,mxp,myp)
    real, intent(in) :: pi0(mzp,mxp,myp)
    real, intent(in) :: th0(mzp,mxp,myp)
    real, intent(in) :: rtgt(mxp,myp)
    real, intent(out) :: a1da2
    real, intent(out) :: amoe(mzp,mxp,myp)
    real, intent(out) :: amof(mzp,mxp,myp)
    real, intent(out) :: acoaa(mzp,mxp,myp)



    integer :: i,j,k
    real :: dt2al2,rdto2cv,dt2al2r,rdtr
    real :: acobb(mzp)
    real :: acocc(mzp)

    ! +--------------------------------------------------------------------+
    ! \   Calculate coefficients for the vertical pressure gradient        \
    ! \     and divergence terms.  These will be combined later for the    \
    ! \     implicit computations.                                         \
    ! +--------------------------------------------------------------------+

    acof=0.0 !LFR because acof(k+1,:,:) isn't made
    if (impl .eq. 1) then
       dt2al2 = dts * .75
       a1da2 = 1. / 3.
    else
       dt2al2 = dts
       a1da2 = 1.
    endif
    rdto2cv = sspct ** 2 * rgas * dts / (2.0 * cv)

    do j = ja,jz
       do i = ia,iz

          !         Coefficient for the vertical pressure gradient term

          dt2al2r = .5 * dt2al2 / rtgt(i,j)
          do k = 1,mzp-1
             acoc(k,i,j) = dt2al2r * dzm(k) * (th0(k,i,j) + th0(k+1,i,j))
          enddo

          !         Coefficients for the vertical divergence term

          rdtr = rdto2cv / rtgt(i,j)
          do k = 2,mzp
             vctr12(k) = dn0(k,i,j) * th0(k,i,j)
             vctr11(k) = rdtr * pi0(k,i,j) * dzt(k) / vctr12(k)
          enddo
          vctr12(1) = dn0(1,i,j) * th0(1,i,j)
          do k = 2,mzp-1
             acof(k,i,j) = -vctr11(k) * (vctr12(k) + vctr12(k+1))
             acog(k,i,j) = vctr11(k) * (vctr12(k) + vctr12(k-1))
          enddo
          acog(mzp,i,j) = vctr11(nzp) * (vctr12(nzp) + vctr12(nz))

          do k = 2,mzp-1
             acoaa(k,i,j) = acoc(k,i,j) * acog(k,i,j)
             acobb(k) = acoc(k,i,j) * (acof(k,i,j) - acog(k+1,i,j)) - 1.
             acocc(k) = -acoc(k,i,j) * acof(k+1,i,j)
          enddo
          acobb(1) = -1.
          acocc(1) = 0.
          acoaa(mzp,i,j) = 0.
          acobb(mzp) = -1.

          amof(1,i,j) = acobb(1)
          amoe(1,i,j) = acocc(1) / amof(1,i,j)
          do k = 2,mzp
             amof(k,i,j) = acobb(k) - acoaa(k,i,j) * amoe(k-1,i,j)
             amoe(k,i,j) = acocc(k) / amof(k,i,j)
          enddo

       enddo
    enddo
    return
  end subroutine coefz


  subroutine rayf_u(mzp,mxp,myp,ia,izu,ja,jz,up,rtgx,topx)
    use mem_grid, only : nfpt, distim, nnz, zmn, ztop, dts, nzp, zt, ngrid, nz
    use mem_scratch, only : vctr2, vctr5
    use ref_sounding, only : u01dn

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: izu
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: up(mzp,mxp,myp)
    real, intent(in) :: rtgx(mxp,myp)
    real, intent(in) :: topx(mxp,myp)

    real :: zmkf,c1,c2
    integer :: kf,i,j,k

    !  This routine calculates rayleigh friction terms velocity and theta_il

    if (nfpt /= 0 .and. distim > 0) then
       kf = nnz(1) - nfpt
       zmkf = zmn(kf,1)
       c1 = 1. / (distim * (ztop - zmkf))
       c2 = dts * c1

       ! u friction

       do j = ja,jz
          do i = ia,izu
             do k = 1,nzp
                vctr2(k) = zt(k) * rtgx(i,j) + topx(i,j)
             enddo
             call htint(nzp,u01dn(1,ngrid),zt,nzp,vctr5,vctr2)
             do k = nz,2,-1
                if (vctr2(k) .le. zmkf) exit
                up(k,i,j) = up(k,i,j) + c2 * (vctr2(k) - zmkf)  &
                     * (vctr5(k) - up(k,i,j))
             enddo
          enddo
       enddo
    end if
  end subroutine rayf_u



  subroutine prdctu(mzp,mxp,myp,ia,izu,ja,jz,  &
       up,ut,pp,th0,f13u,rtgu,rtgt,dxu,topu)

    use mem_grid, only : hw4, dzt, distim, dts, nstbot, itopo

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: izu
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: up(mzp,mxp,myp)
    real, intent(in) :: ut(mzp,mxp,myp)
    real, intent(in) :: pp(mzp,mxp,myp)
    real, intent(in) :: th0(mzp,mxp,myp)
    real, intent(in) :: f13u(mxp,myp)
    real, intent(in) :: rtgu(mxp,myp)
    real, intent(in) :: rtgt(mxp,myp)
    real, intent(in) :: dxu(mxp,myp)
    real, intent(in) :: topu(mxp,myp)

    integer :: i,j,k
    real :: dxl
    real :: vt3da(mzp,mxp,myp)
    real :: dpdx(mzp,mxp,myp)

    !     U prediction

    dpdx = 0.

    !     Calculate acoustic tendency (horizontal pressure gradient)

    do j = ja,jz
       do i = ia,izu
          do k = 1,mzp-1
             vt3da(k,i,j) = (pp(k,i,j) + pp(k+1,i,j)  &
                  + pp(k,i+1,j) + pp(k+1,i+1,j)) * hw4(k)
          enddo
       enddo
    enddo

    do j = ja,jz
       do i = ia,izu
          dxl = dxu(i,j) / rtgu(i,j)
          do k = 2,mzp-1
             dpdx(k,i,j) = -(th0(k,i,j) + th0(k,i+1,j)) * .5  &
                  * ((pp(k,i+1,j) * rtgt(i+1,j) - pp(k,i,j) * rtgt(i,j)) * dxl  &
                  + (vt3da(k,i,j) - vt3da(k-1,i,j)) * dzt(k) * f13u(i,j))
          enddo
       enddo
    enddo

    if (distim .ne. 0.) then
       call rayf_u(mzp,mxp,myp,ia,izu,ja,jz,up,rtgu,topu)
    endif

    do j = 1,myp
       do i = 1,mxp
          do k = 1,mzp
             up(k,i,j) = up(k,i,j) + dts * (dpdx(k,i,j) + ut(k,i,j))
          enddo
       enddo
    enddo

    if (nstbot .eq. 1 .and. itopo .eq. 1) then
       do i = 1,mxp
          do j = 1,myp
             up(1,i,j) = up(2,i,j)
          enddo
       enddo
    end if
  end subroutine prdctu




  subroutine rayf_v(mzp,mxp,myp,ia,iz,ja,jz,vp,rtgx,topx)

    use mem_grid, only : nfpt, distim, nnz, zmn, ztop, &
         dts, nzp, zt, ngrid, nz, jdim, icorflg
    use mem_scratch, only : vctr2, vctr5
    use ref_sounding, only: v01dn

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: vp(mzp,mxp,myp)
    real, intent(in) :: rtgx(mxp,myp)
    real, intent(in) :: topx(mxp,myp)

    real :: zmkf,c1,c2
    integer :: kf,i,j,k

    !  This routine calculates rayleigh friction terms velocity and theta_il

    if (nfpt /= 0 .and. distim > 0 .and. (jdim /= 0 .or. icorflg /= 0)) then
       kf = nnz(1) - nfpt
       zmkf = zmn(kf,1)
       c1 = 1. / (distim * (ztop - zmkf))
       c2 = dts * c1

       ! v friction

       do j = ja,jz
          do i = ia,iz
             do k = 1,nzp
                vctr2(k) = zt(k) * rtgx(i,j) + topx(i,j)
             enddo
             call htint(nzp,v01dn(1,ngrid),zt,nzp,vctr5,vctr2)
             do k = nz,2,-1
                if (vctr2(k) .le. zmkf) exit
                vp(k,i,j) = vp(k,i,j) + c2 * (vctr2(k) - zmkf)  &
                     * (vctr5(k) - vp(k,i,j))
             enddo
          enddo
       enddo
    end if
  end subroutine rayf_v





  subroutine prdctv(mzp,mxp,myp,ia,iz,ja,jzv,  &
       vp,vt,pp,th0,f23v,rtgv,rtgt,dyv,topv)

    use mem_grid, only : hw4, dzt, distim, dts, nstbot, itopo, jdim

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jzv
    real, intent(inout) :: vp(mzp,mxp,myp)
    real, intent(in) :: vt(mzp,mxp,myp)
    real, intent(in) :: pp(mzp,mxp,myp)
    real, intent(in) :: th0(mzp,mxp,myp)
    real, intent(in) :: f23v(mxp,myp)
    real, intent(in) :: rtgv(mxp,myp)
    real, intent(in) :: rtgt(mxp,myp)
    real, intent(in) :: dyv(mxp,myp)
    real, intent(in) :: topv(mxp,myp)

    integer :: i,j,k
    real :: dyl
    real :: vt3da(mzp,mxp,myp)
    real :: dpdy(mzp,mxp,myp)

    !     V prediction

    dpdy = 0.

    if (jdim .eq. 1) then

       ! calculate acoustic tendency (horizontal pressure gradient)

       do j = ja,jzv
          do i = ia,iz
             do k = 1,mzp-1
                vt3da(k,i,j) =(pp(k,i,j) + pp(k+1,i,j)  &
                     + pp(k,i,j+1) + pp(k+1,i,j+1)) * hw4(k)
             enddo
          enddo
       enddo

       do j = ja,jzv
          do i = ia,iz
             dyl = dyv(i,j) / rtgv(i,j)
             do k = 2,mzp-1
                dpdy(k,i,j) = -(th0(k,i,j) + th0(k,i,j+1)) * .5  &
                     * ((pp(k,i,j+1)*rtgt(i,j+1) - pp(k,i,j)*rtgt(i,j)) * dyl  &
                     + (vt3da(k,i,j) - vt3da(k-1,i,j)) * dzt(k) * f23v(i,j))
             enddo
          enddo
       enddo
    endif

    if (distim .ne. 0.) then
       call rayf_v(mzp,mxp,myp,ia,iz,ja,jzv,vp,rtgv,topv)
    endif

    do j = 1,myp
       do i = 1,mxp
          do k = 1,mzp
             vp(k,i,j) = vp(k,i,j) + dts * (dpdy(k,i,j) + vt(k,i,j))
          enddo
       enddo
    enddo

    if (nstbot .eq. 1 .and. itopo .eq. 1) then
       do i = 1,mxp
          do j = 1,myp
             vp(1,i,j) = vp(2,i,j)
          enddo
       enddo
    end if
  end subroutine prdctv




  subroutine rayf_w(mzp,mxp,myp,ia,iz,ja,jz,wp,rtgx,topx)

    use mem_grid, only : nfpt, distim, nnz, zmn, ztop, dts, zt, nz
    use mem_scratch, only : vctr2

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: wp(mzp,mxp,myp)
    real, intent(in) :: rtgx(mxp,myp)
    real, intent(in) :: topx(mxp,myp)

    real :: zmkf,c1,c2
    integer :: kf,i,j,k

    !  This routine calculates rayleigh friction terms velocity and theta_il

    if (nfpt /= 0 .and. distim > 0) then
       kf = nnz(1) - nfpt
       zmkf = zmn(kf,1)
       c1 = 1. / (distim * (ztop - zmkf))
       c2 = dts * c1

       !     W friction

       do j = ja,jz
          do i = ia,iz
             do k = nz,2,-1
                vctr2(k) = zt(k) * rtgx(i,j) + topx(i,j)
                if (vctr2(k) .le. zmkf) exit
                wp(k,i,j) = wp(k,i,j) - c2 * (vctr2(k) - zmkf) * wp(k,i,j)
             enddo
          enddo
       enddo
    end if
  end subroutine rayf_w





  subroutine prdctw1(mzp,mxp,myp,ia,iz,ja,jz,  &
       wp,wt,pp,acoc,a1da2,rtgt,topt)

    use mem_grid, only : distim, dts

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: wp(mzp,mxp,myp)
    real, intent(in) :: wt(mzp,mxp,myp)
    real, intent(in) :: pp(mzp,mxp,myp)
    real, intent(in) :: acoc(mzp,mxp,myp)
    real, intent(in) :: a1da2
    real, intent(in) :: rtgt(mxp,myp)
    real, intent(in) :: topt(mxp,myp)

    integer :: i,j,k

    !     First part of prediction at I,J point

    !     Compute forward part of Crank-Nickelson scheme. This will be total
    !     W prediction for explicit case.

    if (distim .ne. 0.) then
       call rayf_w(mzp,mxp,myp,ia,iz,ja,jz,wp,rtgt,topt)
    endif

    do j = 1,myp
       do i = 1,mxp
          do k = 1,mzp-2
             wp(k,i,j) = wp(k,i,j) + dts * wt(k,i,j)
          enddo
       enddo
    enddo

    do j = ja,jz
       do i = ia,iz
          do k = 1,mzp-2
             wp(k,i,j) = wp(k,i,j) + &
                  a1da2 * acoc(k,i,j) * (pp(k,i,j)-pp(k+1,i,j))
          enddo
       enddo
    enddo
  end subroutine prdctw1


  subroutine prdctp1(mzp,mxp,myp,ia,iz,ja,jz,  &
       pp,up,vp,pi0,dn0,th0,pt,f13t,f23t,rtgt,rtgu,rtgv,  &
       heatfx1,fmapui,fmapvi,dxt,dyt,fmapt)

    use mem_grid, only : sspct, & !intent(in)
         jdim,                  & !intent(in)
         hw4,                   & !intent(in)
         dzt,                   & !intent(in)
         dts                      !intent(in)
    use rconstants, only : rocv   !intent(in)

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: pp(mzp,mxp,myp)
    real, intent(in) :: up(mzp,mxp,myp)
    real, intent(in) :: vp(mzp,mxp,myp)
    real, intent(in) :: pi0(mzp,mxp,myp)
    real, intent(in) :: dn0(mzp,mxp,myp)
    real, intent(in) :: th0(mzp,mxp,myp)
    real, intent(in) :: pt(mzp,mxp,myp)
    real, intent(in) :: f13t(mxp,myp)
    real, intent(in) :: f23t(mxp,myp)
    real, intent(in) :: rtgt(mxp,myp)
    real, intent(in) :: rtgu(mxp,myp)
    real, intent(in) :: rtgv(mxp,myp)
    real, intent(out) :: heatfx1(mxp,myp)
    real, intent(in) :: fmapui(mxp,myp)
    real, intent(in) :: fmapvi(mxp,myp)
    real, intent(in) :: dxt(mxp,myp)
    real, intent(in) :: dyt(mxp,myp)
    real, intent(in) :: fmapt(mxp,myp)

    ! Local Variables
    integer :: i, j, k
    real :: rocvpct
    real :: heatdv(mzp,mxp,myp)
    real :: heatfx(mzp,mxp,myp)

    heatdv = 0.
    rocvpct =rocv *sspct ** 2

    !     Divergence calculations for topographical transformation

    !     First calculate vertically transformed heat flux

    do j = ja,jz
       do i = ia,iz
          do k = 1,mzp
             heatfx(k,i,j) = ((up(k,i,j) + up(k,i-1,j)) * f13t(i,j)  &
                  + (vp(k,i,j) + vp(k,i,j-jdim)) * f23t(i,j)  &
                  ) * dn0(k,i,j) * th0(k,i,j)
          enddo
       enddo
    enddo
    do j = ja,jz
       do i = ia,iz
          do k = 1,mzp-1
             heatfx(k,i,j) = (heatfx(k,i,j) + heatfx(k+1,i,j)) * hw4(k)
          enddo
          heatfx1(i,j) = heatfx(1,i,j) / (.5 * (dn0(1,i,j) * th0(1,i,j)  &
               + dn0(2,i,j) * th0(2,i,j)))
       enddo
    enddo

    do j = ja,jz
       do i = ia,iz
          do k = 2,mzp-1
             heatdv(k,i,j) = (heatfx(k,i,j) - heatfx(k-1,i,j)) * dzt(k)
          enddo
       enddo
    enddo
    do j = 1,myp
       do i = 1,mxp
          do k = 1,mzp
             heatfx(k,i,j) = dn0(k,i,j) * th0(k,i,j)
          enddo
       enddo
    enddo
    do j = ja,jz
       do i = ia,iz
          do k = 2,mzp-1

             heatdv(k,i,j) = -rocvpct * pi0(k,i,j) / heatfx(k,i,j)  &
                  * (heatdv(k,i,j) + fmapt(i,j)  &
                  * ((up(k,i,j) * rtgu(i,j) * fmapui(i,j)  &
                  * (heatfx(k,i,j) + heatfx(k,i+1,j))  &
                  - up(k,i-1,j) * rtgu(i-1,j) * fmapui(i-1,j)  &
                  * (heatfx(k,i,j) + heatfx(k,i-1,j))) * dxt(i,j) * .5  &
                  + (vp(k,i,j) * rtgv(i,j) * fmapvi(i,j)  &
                  * (heatfx(k,i,j) + heatfx(k,i,j+jdim))  &
                  - vp(k,i,j-jdim) * rtgv(i,j-jdim)  &
                  * fmapvi(i,j-jdim)  &
                  * (heatfx(k,i,j) + heatfx(k,i,j-jdim)))  &
                  * dyt(i,j) * .5) / rtgt(i,j))

          enddo
       enddo
    enddo


    do j = ja,jz
       do i = ia,iz
          do k = 1,mzp
             pp(k,i,j) = pp(k,i,j) + (pt(k,i,j) + heatdv(k,i,j)) * dts
          enddo
       enddo
    enddo
  end subroutine prdctp1




  subroutine  prdctw2(mzp,mxp,myp,ia,iz,ja,jz,  &
       wp,pp,acoc,amof,amog,acoaa,rtgt,heatfx1)

    use mem_grid, only: nstbot, nsttop, impl, nz

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: wp(mzp,mxp,myp)
    real, intent(in) :: pp(mzp,mxp,myp)
    real, intent(in) :: acoc(mzp,mxp,myp)
    real, intent(in) :: amof(mzp,mxp,myp)
    real, intent(out) :: amog(mzp,mxp,myp)
    real, intent(in) :: acoaa(mzp,mxp,myp)
    real, intent(in) :: rtgt(mxp,myp)
    real, intent(in) :: heatfx1(mxp,myp)

    integer :: i,j,k

    if (nstbot .eq. 1) then
       do j = ja,jz
          do i = ia,iz
             wp(1,i,j) = -heatfx1(i,j) * rtgt(i,j)
          enddo
       enddo
    endif

    if (nsttop .eq. 1) then
       do j = ja,jz
          do i = ia,iz
             wp(nz,i,j) = 0.
          enddo
       enddo
    endif

    if (impl .eq. 1) then

       !  First implicit part of the w prediction

       do j = ja,jz
          do i = ia,iz
             do k = 2,mzp-2
                wp(k,i,j) = wp(k,i,j) - (pp(k+1,i,j) - pp(k,i,j)) * acoc(k,i,j)
             enddo
          enddo
       enddo

       do j = ja,jz
          do i = ia,iz
             amog(1,i,j) = -wp(1,i,j) / amof(1,i,j)
          enddo
          do k = 2,mzp-2
             do i = ia,iz
                amog(k,i,j) = (-wp(k,i,j) - acoaa(k,i,j) * amog(k-1,i,j))  &
                     / amof(k,i,j)
             enddo
          enddo
       enddo
    endif
  end subroutine prdctw2



  subroutine prdctw3(mzp,mxp,myp,ia,iz,ja,jz,wp,amog,amoe,impl)

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: wp(mzp,mxp,myp)
    real, intent(in) :: amog(mzp,mxp,myp)
    real, intent(in) :: amoe(mzp,mxp,myp)
    integer, intent(in) :: impl

    integer :: i,j,k

    !  Conclusion of implicit w prediction

    if (impl .eq. 1) then
       do k = mzp-2,2,-1
          do j = ja,jz
             do i = ia,iz
                wp(k,i,j) = amog(k,i,j) - amoe(k,i,j) * wp(k+1,i,j)
             enddo
          enddo
       enddo
    endif
  end subroutine prdctw3




  subroutine prdctp2(mzp,mxp,myp,ia,iz,ja,jz,pp,wp,acof,acog)

    use mem_grid, only: nstbot, dzm

    implicit none
    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    integer, intent(in) :: ia
    integer, intent(in) :: iz
    integer, intent(in) :: ja
    integer, intent(in) :: jz
    real, intent(inout) :: pp(mzp,mxp,myp)
    real, intent(in) :: wp(mzp,mxp,myp)
    real, intent(in) :: acof(mzp,mxp,myp)
    real, intent(in) :: acog(mzp,mxp,myp)

    integer :: i,j,k
    real :: dzmr

    !  Finish pressure prediction

    do j = ja,jz
       do i = ia,iz
          do k = 2,mzp-1
             pp(k,i,j) = pp(k,i,j)  &
                  + (wp(k,i,j) * acof(k,i,j) + wp(k-1,i,j) * acog(k,i,j))
          enddo
       enddo
    enddo

    if (nstbot .eq. 1) then
       dzmr = dzm(2) / dzm(1)
       do i = 1,mxp
          do j = 1,myp
             pp(1,i,j) = pp(2,i,j) + (pp(2,i,j) - pp(3,i,j)) * dzmr
          enddo
       enddo
    end if
  end subroutine prdctp2




  subroutine acoustic_new(OneGrid)
    !--------------------------------------------------------------------
    !                  Acoustic terms small time-step driver
    !
    !     This routine calls all the necessary routines to march the model
    !     through the small timesteps.
    !-----------------------------------------------------------------------
    use mem_tend
    use mem_grid
    use mem_basic
    use mem_scratch
    use node_mod
    use ModGrid, only: &
         Grid
    use ModAcoust_adap, only: acoust_adap

    implicit none

    type(Grid), pointer :: OneGrid
    real :: scr2(mzp*mxp*myp)

    scr2=0.0

    if (if_adap == 0) then

       call acoust_new(OneGrid, &
            mzp,mxp,myp, &
            basic_g(ngrid)%dn0(1,1,1),basic_g(ngrid)%pi0(1,1,1),  &
            basic_g(ngrid)%th0(1,1,1),basic_g(ngrid)%up(1,1,1),  &
            basic_g(ngrid)%vp(1,1,1),basic_g(ngrid)%wp(1,1,1),  &
            basic_g(ngrid)%pp(1,1,1),  &
            tend%ut(1),tend%vt(1),tend%wt(1),tend%pt(1),  &
            grid_g(ngrid)%topt(1,1),grid_g(ngrid)%topu(1,1),  &
            grid_g(ngrid)%topv(1,1),grid_g(ngrid)%rtgt(1,1),  &
            grid_g(ngrid)%rtgu(1,1),grid_g(ngrid)%f13u(1,1),  &
            grid_g(ngrid)%dxu(1,1),grid_g(ngrid)%rtgv(1,1),  &
            grid_g(ngrid)%dyv(1,1),  &
            grid_g(ngrid)%f23v(1,1),grid_g(ngrid)%f13t(1,1),  &
            grid_g(ngrid)%f23t(1,1),grid_g(ngrid)%fmapui(1,1),  &
            grid_g(ngrid)%fmapvi(1,1),grid_g(ngrid)%dxt(1,1),  &
            grid_g(ngrid)%dyt(1,1),grid_g(ngrid)%fmapt(1,1))



    else

       call acoust_adap(OneGrid, &
            mzp,mxp,myp   &
            ,grid_g(ngrid)%lpu   (1,1)   ,grid_g(ngrid)%lpv   (1,1)    &
            ,grid_g(ngrid)%lpw   (1,1)   ,scratch%scr1        (1)      &
            ,scr2        (1)     ,scratch%vt3da       (1)      &
            ,scratch%vt3db       (1)     ,scratch%vt3dc       (1)      &
            ,scratch%vt3dd       (1)     ,scratch%vt3de       (1)      &
            ,scratch%vt3df       (1)     ,scratch%vt3dg       (1)      &
            ,scratch%vt3dh       (1)     ,scratch%vt2da       (1)      &
            ,basic_g(ngrid)%dn0  (1,1,1) ,basic_g(ngrid)%pi0  (1,1,1)  &
            ,basic_g(ngrid)%th0  (1,1,1) ,basic_g(ngrid)%up   (1,1,1)  &
            ,basic_g(ngrid)%vp   (1,1,1) ,basic_g(ngrid)%wp   (1,1,1)  &
            ,basic_g(ngrid)%pp   (1,1,1) ,tend%ut             (1)      &
            ,tend%vt             (1)     ,tend%wt             (1)      &
            ,tend%pt             (1)     ,grid_g(ngrid)%dxu   (1,1)    &
            ,grid_g(ngrid)%dyv   (1,1)   ,grid_g(ngrid)%fmapui(1,1)    &
            ,grid_g(ngrid)%fmapvi(1,1)   ,grid_g(ngrid)%dxt   (1,1)    &
            ,grid_g(ngrid)%dyt   (1,1)   ,grid_g(ngrid)%fmapt (1,1)    &
            ,grid_g(ngrid)%aru   (1,1,1) ,grid_g(ngrid)%arv   (1,1,1)  &
            ,grid_g(ngrid)%arw   (1,1,1) ,grid_g(ngrid)%volt  (1,1,1)  &
            ,grid_g(ngrid)%volu  (1,1,1) ,grid_g(ngrid)%volv  (1,1,1)  &
            ,grid_g(ngrid)%volw  (1,1,1)                               )

    endif

    return
  end subroutine acoustic_new



  subroutine acoust_new(OneGrid, mzp,mxp,myp,  &
       dn0,pi0,th0,up,vp,wp,pp,ut,vt,wt,pt,  &
       topt,topu,topv,rtgt,rtgu,f13u,dxu,rtgv,  &
       dyv,f23v,f13t,f23t,fmapui,fmapvi,dxt,dyt,fmapt)
    !--------------------------------------------------------------------
    !                  Acoustic terms small time-step driver
    !
    !     This routine calls all the necessary routines to march the model
    !     through the small timesteps.
    !-----------------------------------------------------------------------

    use ModParallelEnvironment, only: MsgDump

    use ModGrid, only: Grid, DumpGrid

    use ModMessageSet, only: &
         PostRecvSendMsgs, &
         WaitRecvMsgs

    use mem_grid, only : nnacoust, & ! intent(in)
         ngrid,                    & ! intent(in)
         dts,                      & ! intent(out)
         dtlt,                     & ! intent(in)
         nxtnest,                  & ! intent(in)
         nzp, nxp, nyp,            & ! intent(in)
         impl                        ! intent(in)

    use node_mod, only :  nmachs,  & ! intent(in)
         ia,                       & ! intent(in)
         iz,                       & ! intent(in)
         ja,                       & ! intent(in)
         jz,                       & ! intent(in)
         izu,                      & ! intent(in)
         jzv                         ! intent(in)

    implicit none

    type(Grid), pointer :: OneGrid

    integer, intent(in) :: mzp
    integer, intent(in) :: mxp
    integer, intent(in) :: myp
    real, intent(in) :: dn0(mzp,mxp,myp)
    real, intent(in) :: pi0(mzp,mxp,myp)
    real, intent(in) :: th0(mzp,mxp,myp)
    real, intent(inout) :: up(mzp,mxp,myp)
    real, intent(inout) :: vp(mzp,mxp,myp)
    real, intent(inout) :: wp(mzp,mxp,myp)
    real, intent(inout) :: pp(mzp,mxp,myp)
    real, intent(in) :: topt(mxp,myp)
    real, intent(in) :: topu(mxp,myp)
    real, intent(in) :: topv(mxp,myp)
    real, intent(in) :: rtgt(mxp,myp)
    real, intent(in) :: rtgu(mxp,myp)
    real, intent(in) :: f13u(mxp,myp)
    real, intent(in) :: dxu(mxp,myp)
    real, intent(in) :: rtgv(mxp,myp)
    real, intent(in) :: dyv(mxp,myp)
    real, intent(in) :: f23v(mxp,myp)
    real, intent(in) :: f13t(mxp,myp)
    real, intent(in) :: f23t(mxp,myp)
    real, intent(in) :: fmapui(mxp,myp)
    real, intent(in) :: fmapvi(mxp,myp)
    real, intent(in) :: dxt(mxp,myp)
    real, intent(in) :: dyt(mxp,myp)
    real, intent(in) :: fmapt(mxp,myp)
    real, intent(in) :: ut(mzp,mxp,myp) 
    real, intent(in) :: vt(mzp,mxp,myp)
    real, intent(in) :: wt(mzp,mxp,myp)
    real, intent(in) :: pt(mzp,mxp,myp)

    include "tsNames.h"

    integer :: iter
    integer :: lastIter
    logical :: singleProcRun
    logical :: outermostGrid
    real :: a1da2
    real :: acoaa(mzp,mxp,myp)
    real :: acoc(mzp,mxp,myp)
    real :: acof(mzp,mxp,myp)
    real :: acog(mzp,mxp,myp)
    real :: amoe(mzp,mxp,myp)
    real :: amof(mzp,mxp,myp)
    real :: amog(mzp,mxp,myp)
    real :: heatfx1(mxp,myp)
    character(len=*), parameter :: h="**(acoust_new)**"
    logical, parameter :: dumpLocal=.false.

    lastIter = nnacoust(ngrid)
    singleProcRun = nmachs == 1
    outermostGrid = nxtnest(ngrid) == 0

    ! get coefficients for computations

    dts = 2. * dtlt / lastIter
    ! computes acoc, acof, acog, a1da2, amoe, amof, acoaa
    ! uses dn0, pi0, th0, rtgt
    call coefz(mzp,mxp,myp,ia,iz,ja,jz,  &
         acoc,acof,acog,dn0,pi0,th0,rtgt,a1da2,amoe,amof,acoaa)

    ! run small time steps

    do iter=1,lastIter

       ! if not first loop iteration,
       !   receives pp(i+1,j) and pp(i,j+1)
       !   if outermost grid, 
       !     receives pp on full grid boundaries

       if (iter > 1) then
          call WaitRecvMsgs(OneGrid%AcouSendP, OneGrid%AcouRecvP)
       end if

       ! up = f(up, pp); uses pp(i+1,j)
       ! remaining arguments are loop invariant

       call prdctu(mzp,mxp,myp,ia,izu,ja,jz,  &
            up,ut,pp,th0,f13u,rtgu,rtgt,dxu,topu)

       ! if not last loop iteration,
       !   sends up(i-1,j)

       if (iter < lastIter) then
          call PostRecvSendMsgs(OneGrid%AcouSendU, OneGrid%AcouRecvU)
       end if

       ! vp = f(vp, pp); uses pp(i,j+1)
       ! remaining arguments are loop invariant

       call prdctv(mzp,mxp,myp,ia,iz,ja,jzv,&
            vp,vt,pp,th0,f23v,rtgv,rtgt,dyv,topv)

       ! on parallel runs,
       !   if not last loop iteration
       !      sends vp(i,j-1)
       !   else if last loop iteration,
       !      sends up, vp to update all neighbour processes ghost zone

       if (.not. singleProcRun) then
          if (iter < lastIter) then
             call PostRecvSendMsgs(OneGrid%AcouSendV, OneGrid%AcouRecvV)
          else
             call PostRecvSendMsgs(OneGrid%AcouSendUV, OneGrid%AcouRecvUV)
          end if
       end if

       ! wp = f(wp, pp); uses node inner cells only
       ! remaining arguments are loop invariant

       call prdctw1(mzp,mxp,myp,ia,iz,ja,jz, &
            wp,wt,pp,acoc,a1da2,rtgt,topt)

       ! on parallel runs,
       !   if not last loop iteration
       !      receives up(i-1,j)
       !      if outermost grid, 
       !        receives up on full grid boundaries
       !      receives vp(i,j-1)
       !      if outermost grid, 
       !        receives vp on full grid boundaries
       !   else if last loop iteration,
       !      receives up, vp and updates this process ghost zone
       !      if outermost grid, 
       !        receives up, vp on full grid boundaries

       if (.not. singleProcRun) then
          if (iter < lastIter) then
             call WaitRecvMsgs(OneGrid%AcouSendU, OneGrid%AcouRecvU)
             call WaitRecvMsgs(OneGrid%AcouSendV, OneGrid%AcouRecvV)
          else
             call WaitRecvMsgs(OneGrid%AcouSendUV, OneGrid%AcouRecvUV)
          end if
       end if

       ! pp = f(pp, up, vp); uses up(i-1,j), vp(i,j-1)
       ! also computes heatfx1
       ! remaining arguments are loop invariant

       call prdctp1(mzp,mxp,myp,ia,iz,ja,jz,  &
            pp,up,vp,pi0,dn0,th0,pt,f13t,f23t,rtgt,rtgu,rtgv, &
            heatfx1,fmapui,fmapvi,dxt,dyt,fmapt)

       ! wp = f(wp, pp); uses node inner cells only
       ! uses heatfx1
       ! also computes amog
       ! remaining arguments are loop invariant

       call prdctw2(mzp,mxp,myp,ia,iz,ja,jz, &
            wp,pp,acoc,amof,amog,acoaa,rtgt,heatfx1)

       ! finishes updating wp=f(wp); uses node inner cells only
       ! uses amog
       ! remaining arguments are loop invariant

       call prdctw3(mzp,mxp,myp,ia,iz,ja,jz,wp,amog,amoe,impl)

       ! finishes updating pp=f(pp,wp); uses node inner cells only
       ! remaining arguments are loop invariant

       call prdctp2(mzp,mxp,myp,ia,iz,ja,jz,pp,wp,acof,acog)

       ! on parallel runs,
       !   if not last loop iteration
       !      sends pp(i,j+1) and pp(i+1,j)
       !   else if last loop iteration,
       !      sends wp, pp to update neighbour processes ghost zone
       !      receives wp and pp and updates this process ghost zone

       if (.not. singleProcRun) then
          if (iter < lastIter) then
             call PostRecvSendMsgs(OneGrid%AcouSendP, OneGrid%AcouRecvP)
          else
             call PostRecvSendMsgs(OneGrid%AcouSendWP, OneGrid%AcouRecvWP)
             call WaitRecvMsgs(OneGrid%AcouSendWP, OneGrid%AcouRecvWP)
          end if
       end if

    enddo

  end subroutine acoust_new
end module ModAcoust



!******************************************************************************

subroutine buoyancy()

  use mem_tend
  use mem_basic
  use mem_scratch
  use mem_grid
  use node_mod
  use micphys

  implicit none

  call boyanc(mzp,mxp,myp,ia,iz,ja,jz,level                   &
       ,grid_g(ngrid)%lpw   (1,1)   ,tend%wt           (1)      &
       ,basic_g(ngrid)%theta(1,1,1) ,basic_g(ngrid)%rtp(1,1,1)  &
       ,basic_g(ngrid)%rv   (1,1,1) ,basic_g(ngrid)%th0(1,1,1)  &
       ,scratch%vt3da       (1)     ,mynum                      )

  return
end subroutine buoyancy

!******************************************************************************

subroutine boyanc(mzp,mxp,myp,ia,iz,ja,jz,level,lpw  &
     ,wt,theta,rtc,rv,th0,vtemp,mynum)

  use rconstants

  implicit none

  integer :: mzp,mxp,myp,ia,iz,ja,jz,level,mynum
  integer, dimension(mxp,myp) :: lpw
  real, dimension(mzp,mxp,myp) :: wt,theta,rtc,rv,th0,vtemp

  integer :: i,j,k

  if (level .ge. 1) then
     do j = ja,jz
        do i = ia,iz
           do k = lpw(i,j),mzp-1
              vtemp(k,i,j) = gg * ((theta(k,i,j) * (1. + .61 * rv(k,i,j))  &
                   - th0(k,i,j)) / th0(k,i,j) - (rtc(k,i,j) - rv(k,i,j)) )
           enddo
        enddo
     enddo
  else
     do j = ja,jz
        do i = ia,iz
           do k = lpw(i,j),mzp-1
              vtemp(k,i,j) = gg * (theta(k,i,j) / th0(k,i,j) - 1.)
           enddo
        enddo
     enddo
  endif

  do j = ja,jz
     do i = ia,iz
        do k = lpw(i,j),mzp-2
           wt(k,i,j) = wt(k,i,j) + vtemp(k,i,j) + vtemp(k+1,i,j)
        enddo
     enddo
  enddo

  return
end subroutine boyanc