subroutine calc_cpp_mugaz

!==================================================================
!     Purpose
!     -------
!     Check to see if the atmospheric specific heat capacity and
!     mean molar mass for the gas mixture defined in gases.def
!     corresponds to what we're using. If it doesn't, abort run
!     unless option 'check_cpp_match' is set to false in 
!     callphys.def.
!     NOTE: for now, in 1D we do as before. Jeremy, if you're 
!     re-writing rcm1d you may want to alter this.
!
!     Authors
!     ------- 
!     Robin Wordsworth (2009)
!
!==================================================================

      use gases_h
      implicit none

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"

      real cpp_c
      real mugaz_c

      integer igas

      cpp_c   = 0.0
      mugaz_c = 0.0

      do igas=1,ngasmx

         if(igas.eq.vgas)then
            ! ignore variable gas in cpp calculation
         else
            ! all values at 300 K from Engineering Toolbox
            if(gnom(igas).eq.'CO2')then
               cpp_c   = cpp_c   + 0.846*gfrac(igas)
               mugaz_c = mugaz_c + 44.01*gfrac(igas)
            elseif(gnom(igas).eq.'N2_')then
               cpp_c   = cpp_c   + 1.040*gfrac(igas)
               mugaz_c = mugaz_c + 28.01*gfrac(igas)
            elseif(gnom(igas).eq.'H2_')then
               cpp_c   = cpp_c   + 14.31*gfrac(igas)
               mugaz_c = mugaz_c + 2.01*gfrac(igas)
            elseif(gnom(igas).eq.'H2O')then
               cpp_c   = cpp_c   + 1.864*gfrac(igas)
               mugaz_c = mugaz_c + 18.02*gfrac(igas)
            elseif(gnom(igas).eq.'CH4')then
               cpp_c   = cpp_c   + 2.226*gfrac(igas)
               mugaz_c = mugaz_c + 16.04*gfrac(igas)
            elseif(gnom(igas).eq.'NH3')then
               cpp_c   = cpp_c   + 2.175*gfrac(igas)
               mugaz_c = mugaz_c + 17.03*gfrac(igas)
               print*,'WARNING, cpp for NH3 may be for liquid' 
            else
               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
               call abort
            endif
         endif

      enddo

      cpp_c = 1000.0*cpp_c

      print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1'
      print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu'
      print*,'Predefined Cp in physics is ',cpp,'J kg^-1 K^-1'
      print*,'Predefined Mg in physics is ',mugaz,'amu'

      if(ngridmx.eq.1)then
         print*,'Automatically setting cpp & mugaz to calculated values in calc_cpp_mugaz' 
         cpp   = cpp_c
         mugaz = mugaz_c
         R     = 8.314511E+0 *1000.E+0/mugaz
         rcp   = R/cpp
      elseif((cpp.ne.cpp_c) .or. (mugaz.ne.mugaz_c))then
         if(check_cpp_match)then
            print*,'Values do not match!'
            print*,'Either adjust cpp / mugaz via newstart to calculated values,'
            print*,'or set check_cpp_match to .false. in callphys.def.'
            stop
         endif
      endif

      return
    end subroutine calc_cpp_mugaz