source: trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90 @ 2757

      subroutine calc_cpp_mugaz

!==================================================================
!     Purpose
!     -------
!     Computes atmospheric specific heat capacity and
!     mean molar mass for the gas mixture defined in gases.def
!     with standard values.
!     In 1d, one can force the model to adopt these standard values
!     by setting cpp_mugaz_mode=2 in callfis.def.
!
!     Authors
!     ------- 
!     Robin Wordsworth (2009)
!     A. Spiga: make the routine OK with latest changes in rcm1d
!     Jeremy Leconte (2022)
!
!==================================================================

      use gases_h
      use comcstfi_mod, only: cpp, mugaz
      use callkeys_mod, only: cpp_mugaz_mode
      use mod_phys_lmdz_para, only : is_master
      implicit none

      real cpp_c   
      real mugaz_c

      integer igas

      cpp_c   = 0.0
      mugaz_c = 0.0


! compute mugaz
      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(igas.eq.igas_CO2)then
               mugaz_c = mugaz_c + 44.01*gfrac(igas)
            elseif(igas.eq.igas_N2)then
               mugaz_c = mugaz_c + 28.01*gfrac(igas)
            elseif(igas.eq.igas_H2)then
               mugaz_c = mugaz_c + 2.01*gfrac(igas)
            elseif(igas.eq.igas_He)then
               mugaz_c = mugaz_c + 4.003*gfrac(igas)
            elseif(igas.eq.igas_H2O)then
               mugaz_c = mugaz_c + 18.02*gfrac(igas)
            elseif(igas.eq.igas_SO2)then
               mugaz_c = mugaz_c + 64.066*gfrac(igas)
            elseif(igas.eq.igas_H2S)then
               mugaz_c = mugaz_c + 34.08*gfrac(igas)
            elseif(igas.eq.igas_CH4)then
               mugaz_c = mugaz_c + 16.04*gfrac(igas)
            elseif(igas.eq.igas_NH3)then
               mugaz_c = mugaz_c + 17.03*gfrac(igas)
            elseif(igas.eq.igas_C2H6)then 
               ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
               mugaz_c = mugaz_c + 30.07*gfrac(igas)
            elseif(igas.eq.igas_C2H2)then
               ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
               mugaz_c = mugaz_c + 26.04*gfrac(igas)
            else
               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
               call abort
            endif
         endif

      enddo

!compute cpp
      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(igas.eq.igas_CO2)then
               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for
               !Mars conditions) 
               cpp_c   = cpp_c   + 0.846*gfrac(igas)*44.01/mugaz_c
            elseif(igas.eq.igas_N2)then
               cpp_c   = cpp_c   + 1.040*gfrac(igas)*28.01/mugaz_c
            elseif(igas.eq.igas_H2)then
               cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
            elseif(igas.eq.igas_He)then
               cpp_c   = cpp_c   + 5.19*gfrac(igas)*4.003/mugaz_c
            elseif(igas.eq.igas_H2O)then
               cpp_c   = cpp_c   + 1.864*gfrac(igas)*18.02/mugaz_c
            elseif(igas.eq.igas_SO2)then
               cpp_c   = cpp_c   + 0.64*gfrac(igas)*64.066/mugaz_c
            elseif(igas.eq.igas_H2S)then
               cpp_c   = cpp_c   + 1.003*gfrac(igas)*34.08/mugaz_c ! from wikipedia...
            elseif(igas.eq.igas_CH4)then
               cpp_c   = cpp_c   + 2.226*gfrac(igas)*16.04/mugaz_c
            elseif(igas.eq.igas_NH3)then
               cpp_c   = cpp_c   + 2.175*gfrac(igas)*17.03/mugaz_c
               print*,'WARNING, cpp for NH3 may be for liquid'
            elseif(igas.eq.igas_C2H6)then
               ! C2H6
               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
               cpp_c   = cpp_c   + 1.763*gfrac(igas)*30.07/mugaz_c
            elseif(igas.eq.igas_C2H2)then
               ! C2H2
               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
               cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
            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'

      if(((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
           (abs(1.-mugaz/mugaz_c).gt.1.e-6)).and. is_master ) then
         ! Ehouarn: tolerate a small mismatch between computed/stored values
         print*,'--> Values do not match with the predefined one ! (',cpp,',',mugaz,')'
         print*,'    Because cp varies with temperature and that some gases may not appear in gases.def,'
         print*,'    a small discrepancy might be completely normal.'
         print*,'    But you might want to check that!'
         print*,'    If you want to use the values calculated here, adjust cpp / mugaz in the dynamics via newstart (3d)'
         print*,'    or use cpp_mugaz_mode=2 (if you are in 1d).'
      endif

      if (cpp_mugaz_mode==2) then
          if (is_master) print*,'*** cpp_mugaz_mode==2, so setting cpp & mugaz to computations in calc_cpp_mugaz.'
          mugaz = mugaz_c
          cpp = cpp_c
      else
          if (is_master) print*,'*** Leaving cpp & mugaz equal to predefined values'
          if (is_master) print*,'(either from dynamics (cpp_mugaz_mode=0) or callfis (cpp_mugaz_mode=1)).'
      endif


      return
    end subroutine calc_cpp_mugaz