#include <misc.h>
#include <params.h>


subroutine omcalc(rcoslat ,d       ,u       ,v       ,dpsl    , & 2,3
                  dpsm    ,pmid    ,pdel    ,rpmid   ,pbot    , &
                  omga    ,nlon    )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Calculate vertical pressure velocity (omga = dp/dt)
! 
! Method: 
! First evaluate the expressions for omega/p, then rescale to omega at
! the end.
! 
! Author: CCM1
! 
!-----------------------------------------------------------------------
!
! $Id: omcalc.F90,v 1.1.2.1 2002/06/15 13:46:58 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use shr_kind_mod, only: r8 => shr_kind_r8
  use pmgrid
  use pspect
  implicit none
#include <comhyb.h>

!------------------------------Arguments--------------------------------
  integer , intent(in) :: nlon                 ! lonitude dimension
  real(r8), intent(in) :: rcoslat(nlon)        ! 1 / cos(lat)
  real(r8), intent(in) :: d(plond,plev)        ! divergence
  real(r8), intent(in) :: u(plond,plev)        ! zonal wind * cos(lat)
  real(r8), intent(in) :: v(plond,plev)        ! meridional wind * cos(lat)
  real(r8), intent(in) :: dpsl(plond)          ! longitudinal component of grad ln(ps)
  real(r8), intent(in) :: dpsm(plond)          ! latitudinal component of grad ln(ps)
  real(r8), intent(in) :: pmid(plond,plev)     ! mid-level pressures
  real(r8), intent(in) :: pdel(plond,plev)     ! layer thicknesses (pressure)
  real(r8), intent(in) :: rpmid(plond,plev)    ! 1./pmid
  real(r8), intent(in) :: pbot(plond)          ! bottom interface pressure
  real(r8), intent(out):: omga(plond,plev)     ! vertical pressure velocity
!-----------------------------------------------------------------------

!---------------------------Local workspace-----------------------------
  integer i,k               ! longitude, level indices
  real(r8) hkk(plond)       ! diagonal element of hydrostatic matrix
  real(r8) hlk(plond)       ! super diagonal element
  real(r8) suml(plond)      ! partial sum over l = (1, k-1)
  real(r8) vgpk             ! v dot grad ps
  real(r8) tmp              ! vector temporary
!-----------------------------------------------------------------------
!
! Zero partial sum
!
  do i=1,nlon
     suml(i) = 0.
  end do
!
! Pure pressure part: top level
!
  do i=1,nlon
     hkk(i) = 0.5*rpmid(i,1)
     omga(i,1) = -hkk(i)*d(i,1)*pdel(i,1)
     suml(i) = suml(i) + d(i,1)*pdel(i,1)
  end do
!
! sum(k)(v(j)*ps*grad(lnps)*db(j)) part. Not normally invoked since 
! the top layer is normally a pure pressure layer.
!
  if (1>=nprlev) then
     do i=1,nlon
        vgpk = rcoslat(i)*(u(i,1)*dpsl(i) + v(i,1)*dpsm(i))*pbot(i)
        tmp = vgpk*hybd(1)
        omga(i,1) = omga(i,1) + hybm(1)*rpmid(i,1)*vgpk - hkk(i)*tmp
        suml(i) = suml(i) + tmp
     end do
  end if
!
! Integrals to level above bottom
!
  do k=2,plev-1
!
! Pure pressure part
!
     do i=1,nlon
        hkk(i) = 0.5*rpmid(i,k)
        hlk(i) =     rpmid(i,k)
        omga(i,k) = -hkk(i)*d(i,k)*pdel(i,k) - hlk(i)*suml(i)
        suml(i) = suml(i) + d(i,k)*pdel(i,k)
     end do
!
! v(j)*grad(lnps) part
!
     if (k>=nprlev) then
        do i=1,nlon
           vgpk = rcoslat(i)*(u(i,k)*dpsl(i) + v(i,k)*dpsm(i))*pbot(i)
           tmp = vgpk*hybd(k)
           omga(i,k) = omga(i,k) + hybm(k)*rpmid(i,k)*vgpk - hkk(i)*tmp
           suml(i) = suml(i) + tmp
        end do
     end if
  end do
!
! Pure pressure part: bottom level
!
  do i=1,nlon
     hkk(i) = 0.5*rpmid(i,plev)
     hlk(i) =     rpmid(i,plev)
     omga(i,plev) = -hkk(i)*d(i,plev)*pdel(i,plev) - hlk(i)*suml(i)
  end do
!
! v(j)*grad(lnps) part. Normally invoked, but omitted if the model is
! running in pure pressure coordinates throughout (e.g. stratospheric 
! mechanistic model).
!
  if (plev>=nprlev) then
     do i=1,nlon
        vgpk = rcoslat(i)*(u(i,plev)*dpsl(i) + v(i,plev)*dpsm(i))* pbot(i)
        omga(i,plev) = omga(i,plev) + hybm(plev)*rpmid(i,plev)*vgpk - &
             hkk(i)*vgpk*hybd(plev)
     end do
  end if
!
! The above expressions give omega/p. Rescale to omega.
!
  do k=1,plev
     do i=1,nlon
        omga(i,k) = omga(i,k)*pmid(i,k)
     end do
  end do
!
  return
end subroutine omcalc