Index: trunk/UTIL/corrk_exo_k/bin_corrk.py =================================================================== --- trunk/UTIL/corrk_exo_k/bin_corrk.py (revision 3632) +++ trunk/UTIL/corrk_exo_k/bin_corrk.py (revision 3633) @@ -7,5 +7,5 @@ R=100 wnedges=xk.wavenumber_grid_R(10000/100, 10000/.3, R) -print("wnedges range: ", wnedges.min(), wnedges.max()) +print("Spectral resolution asked: wn=", wnedges.min(), wnedges.max(), ", R =", R) ### folder where corrk will be saved @@ -66,13 +66,22 @@ tmp_ktab=xk.Ktable(mol=mol, remove_zeros=True) + print(f"Initial spectral resolution = {tmp_ktab.Nw}") - if logpgrid is None: # use first molecule to set p and t grid + if logpgrid is None or tgrid is None: # use first molecule to set p and t grid logpgrid = tmp_ktab.logpgrid tgrid = tmp_ktab.tgrid - print(f"Remapping: {mol}") - tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid) print(f"Binning: {mol}") tmp_ktab.bin_down(wnedges) + + if (tmp_ktab.logpgrid != logpgrid).any() or (tmp_ktab.tgrid != tgrid).any(): + print("Initial P grid:", tmp_ktab.pgrid) + print("Initial T grid:", tmp_ktab.tgrid) + print(f"Remapping {mol}") + print("New P grid:", np.power(logpgrid,10)) + print("New T grid:", tgrid) + tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid) + + print(tmp_ktab) print(f"Writing: {mol}") tmp_ktab.write_hdf5(dir_out+mol+f"_R{R}.corrk.h5") Index: trunk/UTIL/corrk_exo_k/dace_to_corrk.py =================================================================== --- trunk/UTIL/corrk_exo_k/dace_to_corrk.py (revision 3632) +++ trunk/UTIL/corrk_exo_k/dace_to_corrk.py (revision 3633) @@ -7,7 +7,7 @@ ### choose spectral resolution for corrk -R=500 +R=20000 wnedges=xk.wavenumber_grid_R(10000/100, 10000/.3, R) -print(wnedges.min(), wnedges.max()) +print("Spectral resolution asked: wn=", wnedges.min(), wnedges.max(), ", R =", R) ### set data paths @@ -15,19 +15,31 @@ corrk_dir = "datadir/corrk_data/" dir_out = corrk_dir+f"/R{R}_from_dace/" -ref=xk.Ktable(filename=corrk_dir + '/R500_from_R15000/CH4_R500.corrk.h5', mol="CH4") + +### you can set your desired p and t grid here (otherwise, first molecule is used) +logpgrid = None +tgrid = None +# ref=xk.Ktable(filename=corrk_dir + '/R500_from_R15000/CH4_R500.corrk.h5', mol="CH4") +# logpgrid = ref.logpgrid +# tgrid = ref.tgrid + + ### select (output) molecules and their corresponding (input) isotopologues molecules = { - "CO": "12C-16O", - "C2H2":"12C2-1H2", - "C2H4":"12C2-1H4", + # "CO": "12C-16O", + # "CO2": "12C-16O2", + # "H2O": "1H2-16O", + # # "C2H2":"12C2-1H2", + # # "C2H4":"12C2-1H4", "CH4": "12C-1H4", + # "TiO": "48Ti-16O", + # "VO": "51V-16O", } -xk.Settings().set_mks(True) +# xk.Settings().set_mks(True) xk.Settings().set_log_interp(False) xk.Settings().set_case_sensitive(True) -def dace_to_corrk(molecule=None, isotopologue=None, order=17): +def dace_to_corrk(molecule=None, isotopologue=None, order=17, logpgrid=logpgrid, tgrid=tgrid): if molecule is None: molecule=isotopologue_to_species(isotopologue) @@ -43,14 +55,31 @@ return - tmp=xk.hires_to_ktable(path=mol_corrk, helios=True, wnedges=wnedges, mol=molecule, order=order) # test conversion from a whole folder to tmp + print(f"Reading: {mol}") + hr = xk.Hires_spectrum(glob.glob(mol_corrk+"/*.bin")[0], helios=True) + print(f"Initial spectral resolution = {hr.Nw}") - tmp.kdata = tmp.kdata * 10 * xk.Molar_mass().fetch(molecule) / xk.N_A + tmp_ktab=xk.hires_to_ktable(path=mol_corrk, helios=True, wnedges=wnedges, mol=molecule, order=order) # test conversion from a whole folder to tmp + + tmp_ktab.kdata = tmp_ktab.kdata * 10 * xk.Molar_mass().fetch(molecule) / xk.N_A # cm^2/g * 10000/1000 * kg/mol / (molecule/mol) -> m^2/molecule - tmp.kdata_unit = "m^2/molecule" - tmp.remap_logPT(ref.logpgrid, ref.tgrid) - tmp.remove_zeros() + tmp_ktab.kdata_unit = "m^2/molecule" + tmp_ktab.remove_zeros() - print(tmp) - tmp.write_hdf5(dir_out+f"{molecule}_corrk_dace_{R}.h5") + + if logpgrid is None or tgrid is None: # use first molecule to set p and t grid + logpgrid = tmp_ktab.logpgrid + tgrid = tmp_ktab.tgrid + + if (tmp_ktab.logpgrid != logpgrid).any() or (tmp_ktab.tgrid != tgrid).any(): + print("Initial P grid:", tmp_ktab.pgrid) + print("Initial T grid:", tmp_ktab.tgrid) + print(f"Remapping {mol}") + print("New P grid:", np.power(logpgrid,10)) + print("New T grid:", tgrid) + tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid) + + print(tmp_ktab) + print(f"Writing: {mol}") + tmp_ktab.write_hdf5(dir_out+f"{molecule}_corrk_dace_{R}.h5") # %matplotlib notebook @@ -58,5 +87,5 @@ t_plot=300 fig,ax=plt.subplots(figsize=(9,6.5)) - tmp.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',label='g=1') + tmp_ktab.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',label='g=1') fig.tight_layout() plt.savefig(dir_out+f"corrk_dace_R{R}_{molecule}_p{p_plot}_t{t_plot}.pdf") Index: trunk/UTIL/corrk_exo_k/lmdz_corrk_from_exomol.py =================================================================== --- trunk/UTIL/corrk_exo_k/lmdz_corrk_from_exomol.py (revision 3632) +++ (revision ) @@ -1,86 +1,0 @@ -import exo_k as xk -import numpy as np - -### folder with corrk data for each molecule (see 'xsec_to_corrk.py' to convert cross-sections to correlated-k tables) -xk.Settings().set_search_path("datadir/exomol/taurex_R15000/") - -corrk_dir = "datadir/corrk_data/" -dir_out = corrk_dir+"pluto_ch4_var_co_c2h2_c2h4/" - -### molecules for which we have corrk data -gases = [ - "CO", - "C2H2", - "C2H4", - "CH4", -] - -### set variable gas and its variance -vgas={"CH4": "background"} -x_array=[ - 1.00e-4, - 1.00e-3, - 5.00e-3, - 1.00e-2, - 5.00e-2, - 1.00e-1, - 5.00e-1 -] - -### set constant gases and their values -bg_gas = { - "CO": 5e-4, - "C2H2": 5e-6, # at 100km (Young et al, 2018) - "C2H4": 6e-7, # at 100km (Young et al, 2018) - "C2H6": 2e-5, # at 100km (Young et al, 2018) - "N2":"background" -} - -### select your spectral resolution here -IR_grid = np.array([ 10., 100., 160., 1180., 1273., 1280., 1301., 1308., 1320., 1450., 1520., 1550., 1670., 2080., 2220., 2450., 2750., 2920.]) - -VI_grid = np.array([ 1670., 2080., 2220., 2450., 2750., 2920., 3200., 3750., 4000., 4180., 4260., 4460., 4700., 5000., 5400., 5730., 5930., 6300., 6700., 7700., 8200., 8750., 9250., 35000.]) - -wnedges = np.unique(np.concatenate([IR_grid,VI_grid])) - -### if you'd rather use a resolution R -# R=20 -# IR_grid=xk.wavenumber_grid_R(10.,2920.,R) -# VI_grid=xk.wavenumber_grid_R(1670.,35000.,R) - -### if you want to use an existing lmdz corrk grid -# corrk=corrk_dir+'pluton_ch4var_COfix/17x23' -# wn_VI = np.loadtxt(corrk+'/narrowbands_VI.in',skiprows=1) -# wn_IR = np.loadtxt(corrk+'/narrowbands_IR.in',skiprows=1) -# IR_grid = np.concatenate([wn_IR[:,0],wn_IR[-1:,1]]) -# VI_grid = np.concatenate([wn_VI[:,0],wn_VI[-1:,1]]) - -### Create new corrk to resolution IR_grid, VI_grid -xk.Settings().set_mks(True) -xk.Settings().set_log_interp(False) -xk.Settings().set_case_sensitive(True) - -kdata = xk.Kdatabase(gases) -for gas in kdata.ktables: - # convert Xtables to Ktables - if type(kdata.ktables[gas]) is xk.Xtable: - kdata.ktables[gas] = xk.Ktable(wnedges=wnedges, xtable=kdata.ktables[gas]) - -mix_var_gas=kdata.create_mix_ktable5d(vgas_comp=vgas,bg_comp=bg_gas,x_array=x_array) - -IR_mix_var_gas_ktab=mix_var_gas.copy() -IR_mix_var_gas_ktab.bin_down(IR_grid) - -VI_mix_var_gas_ktab=mix_var_gas.copy() -VI_mix_var_gas_ktab.bin_down(VI_grid) - -IR_mix_var_gas_ktab.write_LMDZ(dir_out,band='IR') -VI_mix_var_gas_ktab.write_LMDZ(dir_out,band='VI') -xk.finalize_LMDZ_dir(dir_out,IR_mix_var_gas_ktab.Nw,VI_mix_var_gas_ktab.Nw) - -### write Q.dat -with open(dir_out+'Q.dat', 'w') as f: - f.write(f"{len(gases)}\n") - f.write("\n".join(gases) + "\n") - f.write(f"{len(x_array)}\n") - f.write("\n".join(map(str, x_array)) + "\n") Index: trunk/UTIL/corrk_exo_k/lmdz_corrk_using_exo_k.py =================================================================== --- trunk/UTIL/corrk_exo_k/lmdz_corrk_using_exo_k.py (revision 3633) +++ trunk/UTIL/corrk_exo_k/lmdz_corrk_using_exo_k.py (revision 3633) @@ -0,0 +1,87 @@ +import exo_k as xk +import numpy as np + +### folder with corrk data for each molecule (see 'xsec_to_corrk.py' to convert cross-sections to correlated-k tables) +xk.Settings().set_search_path("datadir/exomol/taurex_R15000/") +# xk.Settings().set_search_path("datadir/corrk_data/R20000_from_dace/") + +corrk_dir = "datadir/corrk_data/" +dir_out = corrk_dir+"pluto_ch4_var_co_c2h2_c2h4/" + +### molecules for which we have corrk data +gases = [ + "CO", + "C2H2", + "C2H4", + "CH4", +] + +### set variable gas and its variance +vgas={"CH4": "background"} +x_array=[ + 1.00e-4, + 1.00e-3, + 5.00e-3, + 1.00e-2, + 5.00e-2, + 1.00e-1, + 5.00e-1 +] + +### set constant gases and their values +bg_gas = { + "CO": 5e-4, + "C2H2": 5e-6, # at 100km (Young et al, 2018) + "C2H4": 6e-7, # at 100km (Young et al, 2018) + "C2H6": 2e-5, # at 100km (Young et al, 2018) + "N2":"background" +} + +### select your spectral resolution here +IR_grid = np.array([ 10., 100., 160., 1180., 1273., 1280., 1301., 1308., 1320., 1450., 1520., 1550., 1670., 2080., 2220., 2450., 2750., 2920.]) + +VI_grid = np.array([ 1670., 2080., 2220., 2450., 2750., 2920., 3200., 3750., 4000., 4180., 4260., 4460., 4700., 5000., 5400., 5730., 5930., 6300., 6700., 7700., 8200., 8750., 9250., 35000.]) + +wnedges = np.unique(np.concatenate([IR_grid,VI_grid])) + +### if you'd rather use a resolution R +# R=20 +# IR_grid=xk.wavenumber_grid_R(10.,2920.,R) +# VI_grid=xk.wavenumber_grid_R(1670.,35000.,R) + +### if you want to use an existing lmdz corrk grid +# corrk=corrk_dir+'pluton_ch4var_COfix/17x23' +# wn_VI = np.loadtxt(corrk+'/narrowbands_VI.in',skiprows=1) +# wn_IR = np.loadtxt(corrk+'/narrowbands_IR.in',skiprows=1) +# IR_grid = np.concatenate([wn_IR[:,0],wn_IR[-1:,1]]) +# VI_grid = np.concatenate([wn_VI[:,0],wn_VI[-1:,1]]) + +### Create new corrk to resolution IR_grid, VI_grid +xk.Settings().set_mks(True) +xk.Settings().set_log_interp(False) +xk.Settings().set_case_sensitive(True) + +kdata = xk.Kdatabase(gases) +for gas in kdata.ktables: + # convert Xtables to Ktables + if type(kdata.ktables[gas]) is xk.Xtable: + kdata.ktables[gas] = xk.Ktable(wnedges=wnedges, xtable=kdata.ktables[gas]) + +mix_var_gas=kdata.create_mix_ktable5d(vgas_comp=vgas,bg_comp=bg_gas,x_array=x_array) + +IR_mix_var_gas_ktab=mix_var_gas.copy() +IR_mix_var_gas_ktab.bin_down(IR_grid) + +VI_mix_var_gas_ktab=mix_var_gas.copy() +VI_mix_var_gas_ktab.bin_down(VI_grid) + +IR_mix_var_gas_ktab.write_LMDZ(dir_out,band='IR') +VI_mix_var_gas_ktab.write_LMDZ(dir_out,band='VI') +xk.finalize_LMDZ_dir(dir_out,IR_mix_var_gas_ktab.Nw,VI_mix_var_gas_ktab.Nw) + +### write Q.dat +with open(dir_out+'Q.dat', 'w') as f: + f.write(f"{len(gases)}\n") + f.write("\n".join(gases) + "\n") + f.write(f"{len(x_array)}\n") + f.write("\n".join(map(str, x_array)) + "\n") Index: trunk/UTIL/corrk_exo_k/xsec_to_corrk.py =================================================================== --- trunk/UTIL/corrk_exo_k/xsec_to_corrk.py (revision 3632) +++ trunk/UTIL/corrk_exo_k/xsec_to_corrk.py (revision 3633) @@ -7,5 +7,5 @@ R=500 wnedges=xk.wavenumber_grid_R(10000/100, 10000/.3, R) -print("wnedges range: ", wnedges.min(), wnedges.max()) +print("Spectral resolution asked: wn=", wnedges.min(), wnedges.max(), ", R =", R) ### folder where corrk will be saved @@ -63,12 +63,26 @@ print(f"Reading: {mol}") - tmp_ktab=xk.Ktable(wnedges=wnedges, xtable=xk.Xtable(pattern, mol=mol, remove_zeros=True)) + xtable=xk.Xtable(pattern, mol=mol, remove_zeros=True) + print(f"Initial spectral resolution = {xtable.Nw}") - if logpgrid is None: # use first molecule to set p and t grid + tmp_ktab=xk.Ktable(wnedges=wnedges, xtable=xtable) + + # print(f"Binning: {mol}") + # tmp_ktab.bin_down(wnedges) + + + if logpgrid is None or tgrid is None: # use first molecule to set p and t grid logpgrid = tmp_ktab.logpgrid tgrid = tmp_ktab.tgrid - print(f"Remapping: {mol}") - tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid) + if (tmp_ktab.logpgrid != logpgrid).any() or (tmp_ktab.tgrid != tgrid).any(): + print("Initial P grid:", tmp_ktab.pgrid) + print("Initial T grid:", tmp_ktab.tgrid) + print(f"Remapping {mol}") + print("New P grid:", np.power(logpgrid,10)) + print("New T grid:", tgrid) + tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid) + + print(tmp_ktab) print(f"Writing: {mol}") tmp_ktab.write_hdf5(dir_out+mol+f"_R{R}.corrk.h5")