subroutine calc_cpp3d(cppNI,rcpNI,t,p)

!==================================================================
!     Purpose
!     -------
!     Compute the atmospheric specific heat capacity as a 
!     function of pressure and temperature (for CO2 gas only)
!
!     Authors
!     ------- 
!     Robin Wordsworth (2009)
!
!==================================================================

        ! you still need to validate this equation vs something....

      implicit none

#include "comcstfi.h"

      real cp0, dB2dT2
      real cppNI      ! specific heat capacity at const. pressure
      real rcpNI      ! R / specific heat capacity
      real t
      real p
      real calmol
      real bar

      calmol = 94.9784 ! 4.18/(mco2/1000)
      bar    = 100000

      cp0    = (7.70+5.3e-3*t - 8.3e-7*t**2)*calmol
      dB2dT2 = 2.69e-5*t - 0.0098
      cppNI  = cp0 - t*(p/bar)*dB2dT2
      !cppNI  = 1000*(t/460)**0.35 ! Sebastian's version

      rcpNI  = R/cppNI

      return
    end subroutine calc_cpp3d