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 character(len=20),parameter :: myname="calc_cpp_mugaz" 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) ! GG MODIF JAN2019 elseif(igas.eq.igas_CO)then mugaz_c = mugaz_c + 28.01*gfrac(igas) elseif(igas.eq.igas_OCS)then ! https://pubchem.ncbi.nlm.nih.gov/compound/carbonyl_sulfide mugaz_c = mugaz_c + 60.0751*gfrac(igas) elseif(igas.eq.igas_HCl)then mugaz_c = mugaz_c + 36.46*gfrac(igas) elseif(igas.eq.igas_HF)then mugaz_c = mugaz_c + 20.01*gfrac(igas) else if (is_master) then print*,'Error in calc_cpp_mugaz: Gas species not recognised!' endif call abort_physic(myname,'Gas species not recognised!',1) 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 !!!!! MODIF GG JAN 2019 (check source values !!) elseif(igas.eq.igas_CO)then cpp_c = cpp_c + 1.0425*gfrac(igas)*28.01/mugaz_c elseif(igas.eq.igas_OCS)then cpp_c = cpp_c + 0.6909*gfrac(igas)*60.07/mugaz_c elseif(igas.eq.igas_HCl)then cpp_c = cpp_c + 1.7087*gfrac(igas)*36.46/mugaz_c elseif(igas.eq.igas_HF)then cpp_c = cpp_c + 1.4553*gfrac(igas)*20.01/mugaz_c else if (is_master) then print*,'Error in calc_cpp_mugaz: Gas species not recognised!' endif call abort_physic(myname,'Gas species not recognised!',1) endif endif enddo cpp_c = 1000.0*cpp_c if (is_master) then print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1' print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu' endif 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 end subroutine calc_cpp_mugaz