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


subroutine difcor(klev    ,ztodt   ,delps   ,u       ,v       , & 1,3
                  qsave   ,pdel    ,pint    ,t       ,tdif    , &
                  udif    ,vdif    ,nlon    )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Add correction term to t and q horizontal diffusions and
! determine the implied heating rate due to momentum diffusion.
! 
! Method: 
! 1. Add correction term to t and q horizontal diffusions. This term
! provides a partial correction of horizontal diffusion on hybrid (sigma)
! surfaces to horizontal diffusion on pressure surfaces. The appropriate
! function of surface pressure (delps, which already contains the diffusion
! coefficient and the time step) is computed during the transform
! from  spherical harmonic coefficients to grid point values. This term
! can only be applied in the portion of the vertical domain in which
! biharmonic horizontal diffusion is employed. In addition, the term is
! unnecessary on pure pressure levels.
!
! 2. Determine the implied heating rate due to momentum diffusion in order
! to conserve total energy and add it to the temperature.
! Reduce complex matrix (ac) to upper Hessenburg matrix (ac)
!
! Author: D. Williamson
! 
!-----------------------------------------------------------------------
!
! $Id: difcor.F90,v 1.1.2.1 2002/06/15 13:46:56 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use shr_kind_mod, only: r8 => shr_kind_r8
  use pmgrid
  use physconst, only: cpair, cpvir

  implicit none

#include <comctl.h>
#include <comhyb.h>

!------------------------------Arguments--------------------------------

  integer , intent(in) :: klev                  ! k-index of top hybrid level
  integer , intent(in) :: nlon                  ! longitude dimension
  real(r8), intent(in) :: ztodt                 ! twice time step unless nstep = 0
  real(r8), intent(in) :: delps(plond)          ! srf press function for correction
  real(r8), intent(in) :: u(plond,plev)         ! u-wind
  real(r8), intent(in) :: v(plond,plev)         ! v-wind
  real(r8), intent(in) :: qsave(plond,plev)     ! moisture fm prv fcst
  real(r8), intent(in) :: pdel(plond,plev)      ! pdel(k) = pint(k+1) - pint(k)
  real(r8), intent(in) :: pint(plond,plevp)     ! pressure at model interfaces
  real(r8), intent(inout) :: t(plond,plev)      ! temperature
  real(r8), intent(inout) :: tdif(plond,plev)   ! initial/final temperature diffusion
  real(r8), intent(inout) :: udif(plond,plev)   ! initial/final u-momentum diffusion
  real(r8), intent(inout) :: vdif(plond,plev)   ! initial/final v-momentum diffusion

!---------------------------Local workspace-----------------------------

  integer i,k                ! longitude, level indices
  real(r8) tcor(plond,plev)  ! temperature correction term
!-----------------------------------------------------------------------
!
! Compute the pressure surface correction term for horizontal diffusion of
! temperature. 
!
  do k=klev,plev
     if (k==1) then
        do i=1,nlon
           tcor(i,k) = delps(i)*0.5/pdel(i,k)*(hybi(k+1)*(t(i,k+1)-t(i,k)))*pint(i,plevp)
        end do
     else if (k==plev) then
        do i=1,nlon
           tcor(i,k) = delps(i)*0.5/pdel(i,k)*(hybi(k)*(t(i,k)-t(i,k-1)))*pint(i,plevp)
        end do
     else
        do i=1,nlon
           tcor(i,k) = delps(i)*0.5/pdel(i,k)*(hybi(k+1)*(t(i,k+1)-t(i,k)) + &
                       hybi(k  )*(t(i,k)-t(i,k-1)))*pint(i,plevp)
        end do
     end if
  end do
!
! Add the temperture diffusion correction to the diffusive heating term 
! and to the temperature.
!
  if (.not.adiabatic .and. .not.ideal_phys) then
     do k=klev,plev
        do i=1,nlon
           tdif(i,k) = tdif(i,k) + tcor(i,k)/ztodt
           t(i,k) = t(i,k) + tcor(i,k)
        end do
     end do
!
! Convert momentum diffusion tendencies to heating rates in order to 
! conserve internal energy. Add the heating to the temperature and to 
! diffusive heating term.
!
     do k=1,plev
        do i=1,nlon
           t(i,k) = t(i,k) - ztodt * (u(i,k)*udif(i,k) + v(i,k)*vdif(i,k)) / &
                                     (cpair*(1. + cpvir*qsave(i,k)))
           tdif(i,k) = tdif(i,k) - (u(i,k)*udif(i,k) + v(i,k)*vdif(i,k)) / &
                                   (cpair*(1. + cpvir*qsave(i,k)))
        end do
     end do
  end if
 
  return
end subroutine difcor