Index: /trunk/LMDZ.GENERIC/README =================================================================== --- /trunk/LMDZ.GENERIC/README (revision 2654) +++ /trunk/LMDZ.GENERIC/README (revision 2655) @@ -1715,2 +1715,11 @@ Minor bug fix for aerave_new.F when input wavelenght in data file are in descending order. + +== 30/03/2022 == GM +Major changes to CIA interpolation: +1) Add contribution from CH4 (H2-CH4,He-CH4,CH4-CH4) ; +2) Add some tests before interpolation for H2, He and CH4 ; +3) Add the possibility to choose between a normal or equilibrium ortho:para +fraction for CIA H2; +4) Change "strictboundH2H2cia" to the generic "strictboundcia" for H2,He,CH4. +It can be added for others CIA (N2,H,CO2...) if you want. Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 (revision 2654) +++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 (revision 2655) @@ -14,6 +14,6 @@ logical,save :: strictboundcorrk !$OMP THREADPRIVATE(strictboundcorrk) - logical,save :: strictboundH2H2cia -!$OMP THREADPRIVATE(strictboundH2H2cia) + logical,save :: strictboundcia +!$OMP THREADPRIVATE(strictboundcia) logical,save :: enertest @@ -76,6 +76,7 @@ integer,save :: startype integer,save :: versH2H2cia + character(64),save :: H2orthopara_mixture integer,save :: nlayaero -!$OMP THREADPRIVATE(iddist,iaervar,iradia,startype,versH2H2cia,nlayaero) +!$OMP THREADPRIVATE(iddist,iaervar,iradia,startype,versH2H2cia,H2orthopara_mixture,nlayaero) integer,dimension(:),allocatable,save :: aeronlay_choice !$OMP THREADPRIVATE(aeronlay_choice) Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 (revision 2654) +++ /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 (revision 2655) @@ -259,9 +259,9 @@ if (is_master) write(*,*) trim(rname)//& - ": prohibit calculations outside H2H2 CIA T grid?" - strictboundH2H2cia=.true. ! default value - call getin_p("strictboundH2H2cia",strictboundH2H2cia) - if (is_master) write(*,*) trim(rname)//& - ": strictboundH2H2cia = ",strictboundH2H2cia + ": prohibit calculations outside CIA T grid?" + strictboundcia=.true. ! default value + call getin_p("strictboundcia",strictboundcia) + if (is_master) write(*,*) trim(rname)//& + ": strictboundcia = ",strictboundcia if (is_master) write(*,*) trim(rname)//& @@ -304,5 +304,12 @@ if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then call abort_physic(rname,"Error: Choose a correct value (2011 or 2018) for versH2H2cia !",1) - endif + endif + + if (is_master) write(*,*) trim(rname)//& + ": CIA - normal or equilibrium ortho-para mixture for H2?" + H2orthopara_mixture = 'normal' + call getin_p("H2orthopara_mixture",H2orthopara_mixture) + if (is_master) write(*,*)trim(rname)//& + ": H2orthopara_mixture = ",trim(H2orthopara_mixture) if (is_master) write(*,*) trim(rname)//& Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateCH4CH4.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/interpolateCH4CH4.F90 (revision 2655) +++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateCH4CH4.F90 (revision 2655) @@ -0,0 +1,155 @@ +subroutine interpolateCH4CH4(wn,temp,pres,abcoef,firstcall,ind) + + !================================================================== + ! + ! Purpose + ! ------- + ! Calculates the CH4-CH4 CIA opacity, using a lookup table from + ! HITRAN (2011) + ! + ! Authors + ! ------- + ! R. Wordsworth (2011) + ! + !================================================================== + + use callkeys_mod, only: strictboundcia + use datafile_mod, only: datadir + implicit none + + ! input + double precision wn ! wavenumber (cm^-1) + double precision temp ! temperature (Kelvin) + double precision pres ! CH4 partial pressure (Pascals) + + ! output + double precision abcoef ! absorption coefficient (m^-1) + + integer nS,nT + parameter(nS=1018) + parameter(nT=10) + + double precision, parameter :: losch = 2.6867774e19 + ! Loschmit's number (molecule cm^-3 at STP) + ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2 + ! see Richard et al. 2011, JQSRT for details + + double precision amagat + double precision wn_arr(nS) + double precision temp_arr(nT) + double precision abs_arr(nS,nT) + + integer k,iT + logical firstcall + + save wn_arr, temp_arr, abs_arr !read by master + + character*100 dt_file + integer strlen,ios + + character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" + + character*20 bleh + double precision blah, Ttemp + integer nres + integer ind + + if(temp.gt.400)then + if (strictboundcia) then + print*,'Your temperatures are too high for this CH4-CH4 CIA dataset.' + print*,'Please run mixed CH4-CH4 atmospheres below T = 400 K.' + stop + else + print*,'Your temperatures are too high for this CH4-CH4 CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 400 ' + print*,' ... for CH4-CH4 CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 400 + endif + elseif(temp.lt.40)then + if (strictboundcia) then + print*,'Your temperatures are too low for this CH4-CH4 CIA dataset.' + print*,'Please run mixed CH4-CH4 atmospheres above T = 40 K.' + stop + else + print*,'Your temperatures are too low for this CH4-CH4 CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 40 ' + print*,' ... for CH4-CH4 CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 40 + endif + endif + + amagat = (273.15/temp)*(pres/101325.0) + + if(firstcall)then ! called by sugas_corrk only + print*,'----------------------------------------------------' + print*,'Initialising CH4-CH4 continuum from HITRAN database...' + + ! 1.1 Open the ASCII files + dt_file=TRIM(datadir)//'/continuum_data/CH4-CH4_2011.cia' + +!$OMP MASTER + open(33,file=dt_file,form='formatted',status='old',iostat=ios) + if (ios.ne.0) then ! file not found + write(*,*) 'Error from interpolateCH4CH4' + write(*,*) 'Data file ',trim(dt_file),' not found.' + write(*,*) 'Check that your path to datagcm:',trim(datadir) + write(*,*) 'is correct. You can change it in callphys.def with:' + write(*,*) 'datadir = /absolute/path/to/datagcm' + write(*,*) 'Also check that the continuum data continuum_data/CH4-CH4_2011.cia is there.' + call abort + else + + do iT=1,nT + + read(33,fmat1) bleh,blah,blah,nres,Ttemp + if(nS.ne.nres)then + print*,'Resolution given in file: ',trim(dt_file) + print*,'is ',nres,', which does not match nS.' + print*,'Please adjust nS value in interpolateCH4CH4.F90' + stop + endif + temp_arr(iT)=Ttemp + + do k=1,nS + read(33,*) wn_arr(k),abs_arr(k,it) + end do + + end do + + endif + close(33) +!$OMP END MASTER +!$OMP BARRIER + + print*,'interpolateCH4CH4: At wavenumber ',wn,' cm^-1' + print*,' temperature ',temp,' K' + print*,' pressure ',pres,' Pa' + + endif + + call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) + +! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' +! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' + + abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1 + +! print*,'We have ',amagat,' amagats of CH4' +! print*,'So the absorption is ',abcoef,' m^-1' + + +! unlike for Rayleigh scattering, we do not currently weight by the BB function +! however our bands are normally thin, so this is no big deal. + + + return +end subroutine interpolateCH4CH4 + Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2CH4.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2CH4.F90 (revision 2655) +++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2CH4.F90 (revision 2655) @@ -0,0 +1,164 @@ +subroutine interpolateH2CH4(wn,temp,presH2,presCH4,abcoef,firstcall,ind) + + !================================================================== + ! + ! Purpose + ! ------- + ! Calculates the H2-CH4 CIA opacity, using a lookup table from + ! HITRAN (2011) + ! + ! Authors + ! ------- + ! R. Wordsworth (2011) + ! + !================================================================== + + use callkeys_mod, only: H2orthopara_mixture, strictboundcia + use datafile_mod, only: datadir + implicit none + + ! input + double precision wn ! wavenumber (cm^-1) + double precision temp ! temperature (Kelvin) + double precision presCH4 ! CH4 partial pressure (Pascals) + double precision presH2 ! H2 partial pressure (Pascals) + + ! output + double precision abcoef ! absorption coefficient (m^-1) + + integer nS,nT + parameter(nS=1974) + parameter(nT=10) + + double precision, parameter :: losch = 2.6867774e19 + ! Loschmit's number (molecule cm^-3 at STP) + ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2 + ! see Richard et al. 2011, JQSRT for details + + double precision amagatCH4 + double precision amagatH2 + double precision wn_arr(nS) + double precision temp_arr(nT) + double precision abs_arr(nS,nT) + + integer k,iT + logical firstcall + + save wn_arr, temp_arr, abs_arr !read by master + + character*100 dt_file + integer strlen,ios + + character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" + + character*20 bleh + double precision blah, Ttemp + integer nres + integer ind + + if(temp.gt.400)then + if (strictboundcia) then + print*,'Your temperatures are too high for this H2-CH4 CIA dataset.' + print*,'Please run mixed H2-CH4 atmospheres below T = 400 K.' + stop + else + print*,'Your temperatures are too high for this H2-CH4 CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 400 ' + print*,' ... for H2-CH4 CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 400 + endif + elseif(temp.lt.40)then + if (strictboundcia) then + print*,'Your temperatures are too low for this H2-CH4 CIA dataset.' + print*,'Please run mixed H2-CH4 atmospheres above T = 40 K.' + stop + else + print*,'Your temperatures are too low for this H2-CH4 CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 40 ' + print*,' ... for H2-CH4 CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 40 + endif + endif + + amagatCH4 = (273.15/temp)*(presCH4/101325.0) + amagatH2 = (273.15/temp)*(presH2/101325.0) + + if(firstcall)then ! called by sugas_corrk only + print*,'----------------------------------------------------' + print*,'Initialising H2-CH4 continuum from HITRAN database...' + + ! 1.1 Open the ASCII files + if (H2orthopara_mixture.eq."normal") then + dt_file=TRIM(datadir)//'/continuum_data/H2-CH4_norm_2011.cia' + else if (H2orthopara_mixture.eq."equilibrium") then + dt_file=TRIM(datadir)//'/continuum_data/H2-CH4_eq_2011.cia' + endif + +!$OMP MASTER + open(33,file=dt_file,form='formatted',status='old',iostat=ios) + if (ios.ne.0) then ! file not found + write(*,*) 'Error from interpolateH2CH4' + write(*,*) 'Data file ',trim(dt_file),' not found.' + write(*,*) 'Check that your path to datagcm:',trim(datadir) + write(*,*) 'is correct. You can change it in callphys.def with:' + write(*,*) 'datadir = /absolute/path/to/datagcm' + write(*,*) 'Also check that the continuum data continuum_data/H2-CH4_norm_2011.cia is there.' + call abort + else + + do iT=1,nT + + read(33,fmat1) bleh,blah,blah,nres,Ttemp + if(nS.ne.nres)then + print*,'Resolution given in file: ',trim(dt_file) + print*,'is ',nres,', which does not match nS.' + print*,'Please adjust nS value in interpolateH2CH4.F90' + stop + endif + temp_arr(iT)=Ttemp + + do k=1,nS + read(33,*) wn_arr(k),abs_arr(k,it) + end do + + end do + + endif + close(33) +!$OMP END MASTER +!$OMP BARRIER + + print*,'interpolateH2CH4: At wavenumber ',wn,' cm^-1' + print*,' temperature ',temp,' K' + print*,' CH4 partial pressure ',presCH4,' Pa' + print*,' and H2 partial pressure ',presH2,' Pa' + + endif + + call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) + +! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' +! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' + + abcoef=abcoef*losch**2*100.0*amagatCH4*amagatH2 ! convert to m^-1 + +! print*,'We have ',amagatCH4,' amagats of CH4' +! print*,'and ',amagatH2,' amagats of H2' +! print*,'So the absorption is ',abcoef,' m^-1' + + +! unlike for Rayleigh scattering, we do not currently weight by the BB function +! however our bands are normally thin, so this is no big deal. + + + return +end subroutine interpolateH2CH4 + Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90 (revision 2654) +++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90 (revision 2655) @@ -1,2 +1,3 @@ + subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,ind) @@ -17,5 +18,5 @@ !================================================================== - use callkeys_mod, only: versH2H2cia, strictboundH2H2cia + use callkeys_mod, only: versH2H2cia, strictboundcia, H2orthopara_mixture use datafile_mod, only: datadir @@ -60,21 +61,37 @@ integer ind - + if(temp.gt.400)then - if (strictboundH2H2cia) then - print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you ' - print*,'really want to run simulations with hydrogen at T > 400 K, contact' - print*,'Robin Wordsworth [rwordsworth@uchicago.edu].' - stop - else - print*,'Your temperatures are too high for this H2-H2 CIA dataset' - print*,'you have chosen strictboundH2H2cia = ', strictboundH2H2cia - print*,'*********************************************************' - print*,' we allow model to continue but with temp = 400 ' - print*,' ... for H2-H2 CIA dataset ... ' - print*,' ... we assume we know what you are doing ... ' - print*,'*********************************************************' - temp = 400 - endif + if (strictboundcia) then + print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you ' + print*,'really want to run simulations with hydrogen at T > 400 K, contact' + print*,'Robin Wordsworth [rwordsworth@uchicago.edu].' + stop + else + print*,'Your temperatures are too high for this H2-H2 CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 400 ' + print*,' ... for H2-H2 CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 400 + endif + elseif(temp.lt.40)then + if (strictboundcia) then + print*,'Your temperatures are too low for this H2-H2 CIA dataset. If you ' + print*,'really want to run simulations with hydrogen at T < 40 K, contact' + print*,'Robin Wordsworth [rwordsworth@uchicago.edu].' + stop + else + print*,'Your temperatures are too low for this H2-H2 CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 40 ' + print*,' ... for H2-H2 CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 40 + endif endif @@ -88,9 +105,17 @@ ! Only two possible versions for now : 2011 or 2018 (sanity check in inifis_mod) if (versH2H2cia.eq.2011) then - dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia' nS = 2428 + if (H2orthopara_mixture.eq."normal") then + dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia' + else if (H2orthopara_mixture.eq."equilibrium") then + dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2011.cia' + endif else if (versH2H2cia.eq.2018) then - dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2018.cia' nS = 9600 + if (H2orthopara_mixture.eq."normal") then + dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2018.cia' + else if (H2orthopara_mixture.eq."equilibrium") then + dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2018.cia' + endif endif @@ -106,5 +131,5 @@ write(*,*) 'is correct. You can change it in callphys.def with:' write(*,*) 'datadir = /absolute/path/to/datagcm' - write(*,*) 'Also check that the continuum data continuum_data/H2-H2_norm_2011.cia or H2-H2_norm_2018.cia is there.' + write(*,*) 'Also check that the continuum data is there.' call abort else @@ -113,5 +138,5 @@ write(*,*) '... You are using H2-H2 CIA from 2011 but you should use more recent data available on HITRAN !' write(*,*) '... (Especially if you are running a giant planet atmosphere)' - write(*,*) '... Just find out the H2-H2_norm_2018.cia, put it in your datadir and have a look at interpolateH2H2.F90 ! .' + write(*,*) '... Just find out the H2-H2 CIA from 2018, put it in your datadir and have a look at interpolateH2H2.F90 ! .' endif @@ -149,5 +174,7 @@ endif - + + + call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90 (revision 2654) +++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90 (revision 2655) @@ -14,4 +14,5 @@ !================================================================== + use callkeys_mod, only: H2orthopara_mixture, strictboundcia use datafile_mod, only: datadir @@ -57,9 +58,35 @@ integer ind - + if(temp.gt.400)then - print*,'Your temperatures are too high for this H2-He CIA dataset.' - print*,'Please run mixed H2-He atmospheres below T = 400 K.' - stop + if (strictboundcia) then + print*,'Your temperatures are too high for this H2-He CIA dataset.' + print*,'Please run mixed H2-He atmospheres below T = 400 K.' + stop + else + print*,'Your temperatures are too high for this H2-He CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 400 ' + print*,' ... for H2-He CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 400 + endif + elseif(temp.lt.40)then + if (strictboundcia) then + print*,'Your temperatures are too low for this H2-He CIA dataset.' + print*,'Please run mixed H2-He atmospheres above T = 40 K.' + stop + else + print*,'Your temperatures are too low for this H2-He CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 40 ' + print*,' ... for H2-He CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 40 + endif endif @@ -72,5 +99,9 @@ ! 1.1 Open the ASCII files - dt_file=TRIM(datadir)//'/continuum_data/H2-He_norm_2011.cia' + if (H2orthopara_mixture.eq."normal") then + dt_file=TRIM(datadir)//'/continuum_data/H2-He_norm_2011.cia' + else if (H2orthopara_mixture.eq."equilibrium") then + dt_file=TRIM(datadir)//'/continuum_data/H2-He_eq_2011.cia' + endif !$OMP MASTER @@ -82,5 +113,5 @@ write(*,*) 'is correct. You can change it in callphys.def with:' write(*,*) 'datadir = /absolute/path/to/datagcm' - write(*,*) 'Also check that the continuum data continuum_data/H2-He_norm_2011.cia is there.' + write(*,*) 'Also check that the continuum data is there.' call abort else Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateHeCH4.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/interpolateHeCH4.F90 (revision 2655) +++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateHeCH4.F90 (revision 2655) @@ -0,0 +1,160 @@ +subroutine interpolateHeCH4(wn,temp,presHe,presCH4,abcoef,firstcall,ind) + + !================================================================== + ! + ! Purpose + ! ------- + ! Calculates the He-CH4 CIA opacity, using a lookup table from + ! HITRAN (2018) + ! + ! Authors + ! ------- + ! R. Wordsworth (2011) + ! + !================================================================== + + use callkeys_mod, only: strictboundcia + use datafile_mod, only: datadir + implicit none + + ! input + double precision wn ! wavenumber (cm^-1) + double precision temp ! temperature (Kelvin) + double precision presCH4 ! CH4 partial pressure (Pascals) + double precision presHe ! He partial pressure (Pascals) + + ! output + double precision abcoef ! absorption coefficient (m^-1) + + integer nS,nT + parameter(nS=1000) + parameter(nT=10) + + double precision, parameter :: losch = 2.6867774e19 + ! Loschmit's number (molecule cm^-3 at STP) + ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2 + ! see Richard et al. 2011, JQSRT for details + + double precision amagatCH4 + double precision amagatHe + double precision wn_arr(nS) + double precision temp_arr(nT) + double precision abs_arr(nS,nT) + + integer k,iT + logical firstcall + + save wn_arr, temp_arr, abs_arr !read by master + + character*100 dt_file + integer strlen,ios + + character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" + + character*20 bleh + double precision blah, Ttemp + integer nres + integer ind + + if(temp.gt.350)then + if (strictboundcia) then + print*,'Your temperatures are too high for this He-CH4 CIA dataset.' + print*,'Please run mixed He-CH4 atmospheres below T = 350 K.' + stop + else + print*,'Your temperatures are too high for this He-CH4 CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 350 ' + print*,' ... for He-CH4 CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 350 + endif + elseif(temp.lt.40)then + if (strictboundcia) then + print*,'Your temperatures are too low for this He-CH4 CIA dataset.' + print*,'Please run mixed He-CH4 atmospheres above T = 40 K.' + stop + else + print*,'Your temperatures are too low for this He-CH4 CIA dataset' + print*,'you have chosen strictboundcia = ', strictboundcia + print*,'*********************************************************' + print*,' we allow model to continue but with temp = 40 ' + print*,' ... for He-CH4 CIA dataset ... ' + print*,' ... we assume we know what you are doing ... ' + print*,'*********************************************************' + temp = 40 + endif + endif + + amagatCH4 = (273.15/temp)*(presCH4/101325.0) + amagatHe = (273.15/temp)*(presHe/101325.0) + + if(firstcall)then ! called by sugas_corrk only + print*,'----------------------------------------------------' + print*,'Initialising He-CH4 continuum from HITRAN database...' + + ! 1.1 Open the ASCII files + dt_file=TRIM(datadir)//'/continuum_data/CH4-He_2018.cia' + +!$OMP MASTER + open(33,file=dt_file,form='formatted',status='old',iostat=ios) + if (ios.ne.0) then ! file not found + write(*,*) 'Error from interpolateHeCH4' + write(*,*) 'Data file ',trim(dt_file),' not found.' + write(*,*) 'Check that your path to datagcm:',trim(datadir) + write(*,*) 'is correct. You can change it in callphys.def with:' + write(*,*) 'datadir = /absolute/path/to/datagcm' + write(*,*) 'Also check that the continuum data continuum_data/He-CH4_2018.cia is there.' + call abort + else + + do iT=1,nT + + read(33,fmat1) bleh,blah,blah,nres,Ttemp + if(nS.ne.nres)then + print*,'Resolution given in file: ',trim(dt_file) + print*,'is ',nres,', which does not match nS.' + print*,'Please adjust nS value in interpolateHeCH4.F90' + stop + endif + temp_arr(iT)=Ttemp + + do k=1,nS + read(33,*) wn_arr(k),abs_arr(k,it) + end do + + end do + + endif + close(33) +!$OMP END MASTER +!$OMP BARRIER + + print*,'interpolateHeCH4: At wavenumber ',wn,' cm^-1' + print*,' temperature ',temp,' K' + print*,' CH4 partial pressure ',presCH4,' Pa' + print*,' and He partial pressure ',presHe,' Pa' + + endif + + call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) + +! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' +! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' + + abcoef=abcoef*losch**2*100.0*amagatCH4*amagatHe ! convert to m^-1 + +! print*,'We have ',amagatCH4,' amagats of CH4' +! print*,'and ',amagatHe,' amagats of He' +! print*,'So the absorption is ',abcoef,' m^-1' + + +! unlike for Rayleigh scattering, we do not currently weight by the BB function +! however our bands are normally thin, so this is no big deal. + + + return +end subroutine interpolateHeCH4 + Index: /trunk/LMDZ.GENERIC/libf/phystd/optci.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/optci.F90 (revision 2654) +++ /trunk/LMDZ.GENERIC/libf/phystd/optci.F90 (revision 2655) @@ -12,5 +12,5 @@ L_NLEVRAD, L_REFVAR, naerkind use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig - use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2 + use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, igas_CH4 use comcstfi_mod, only: g, r, mugaz use callkeys_mod, only: kastprof,continuum,graybody @@ -208,4 +208,27 @@ enddo + elseif(igas.eq.igas_CH4)then + + ! first do self-induced absorption + interm = indi(nw,igas,igas) + call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indi(nw,igas,igas) = interm + + ! then cross-interactions with other gases + do jgas=1,ngasmx + p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) + dtempc = 0.0d0 + if(jgas.eq.igas_H2)then + interm = indi(nw,igas,jgas) + call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indi(nw,igas,jgas) = interm + elseif(jgas.eq.igas_He)then + interm = indi(nw,igas,jgas) + call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indi(nw,igas,jgas) = interm + endif + dtemp = dtemp + dtempc + enddo + elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then ! Compute self and foreign (with air) continuum of H2O Index: /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 (revision 2654) +++ /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 (revision 2655) @@ -11,5 +11,5 @@ use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig - use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2 + use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4 use comcstfi_mod, only: g, r, mugaz use callkeys_mod, only: kastprof,continuum,graybody,callgasvis @@ -215,4 +215,27 @@ dtemp = dtemp + dtempc enddo + + elseif(igas.eq.igas_CH4)then + + ! first do self-induced absorption + interm = indv(nw,igas,igas) + call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indv(nw,igas,igas) = interm + + ! then cross-interactions with other gases + do jgas=1,ngasmx + p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) + dtempc = 0.0d0 + if(jgas.eq.igas_H2)then + interm = indv(nw,igas,jgas) + call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indv(nw,igas,jgas) = interm + elseif(jgas.eq.igas_He)then + interm = indv(nw,igas,jgas) + call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indv(nw,igas,jgas) = interm + endif + dtemp = dtemp + dtempc + enddo elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then Index: /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 (revision 2654) +++ /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 (revision 2655) @@ -30,5 +30,5 @@ use comcstfi_mod, only: mugaz use gases_h, only: gnom, ngasmx, & - igas_H2, igas_H2O, igas_He, igas_N2 + igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4 use ioipsl_getin_p_mod, only: getin_p use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,& @@ -737,4 +737,22 @@ endif enddo + + + elseif (igas .eq. igas_CH4) then + + ! first do self-induced absorption + dummy = -9999 + call interpolateCH4CH4(600.0,66.7,400.0,testcont,.true.,dummy) + ! then cross-interactions with other gases + do jgas=1,ngasmx + if (jgas .eq. igas_H2) then + dummy = -9999 + call interpolateH2CH4(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) + elseif (jgas .eq. igas_He) then + dummy = -9999 + call interpolateHeCH4(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) + endif + enddo + elseif (igas .eq. igas_H2O) then