Index: trunk/LMDZ.GENERIC/changelog.txt =================================================================== --- trunk/LMDZ.GENERIC/changelog.txt (revision 3640) +++ trunk/LMDZ.GENERIC/changelog.txt (revision 3641) @@ -1994,2 +1994,6 @@ Follow-up of recent changes in iostart (making 1D restarts) from r3562, Correspondingly adapt newstart. + +== 23/02/2025 == MT +Implement the new generic continuum database +More details here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Continuum_Database Index: trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90 (revision 3640) +++ trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90 (revision 3641) @@ -45,4 +45,6 @@ elseif(igas.eq.igas_N2)then mugaz_c = mugaz_c + 28.01*gfrac(igas) + elseif(igas.eq.igas_O2)then + mugaz_c = mugaz_c + 31.999*gfrac(igas) elseif(igas.eq.igas_H2)then mugaz_c = mugaz_c + 2.01*gfrac(igas) @@ -98,4 +100,6 @@ elseif(igas.eq.igas_N2)then cpp_c = cpp_c + 1.040*gfrac(igas)*28.01/mugaz_c + elseif(igas.eq.igas_O2)then + cpp_c = cpp_c + 0.918*gfrac(igas)*31.999/mugaz_c elseif(igas.eq.igas_H2)then cpp_c = cpp_c + 14.31*gfrac(igas)*2.01/mugaz_c Index: trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90 (revision 3640) +++ trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90 (revision 3641) @@ -29,5 +29,6 @@ use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep use gases_h, only: ngasmx, vgas, gnom, gfrac, igas_CO2, igas_H2, & - igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO + igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO, & + igas_O2 use comcstfi_mod, only: g, mugaz, pi Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 (revision 3640) +++ trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 (revision 3641) @@ -13,6 +13,6 @@ logical,save :: season,diurnal,tlocked,rings_shadow,lwrite !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite) - logical,save :: callgasvis,continuum,graybody -!$OMP THREADPRIVATE(callgasvis,continuum,graybody) + logical,save :: callgasvis,continuum,generic_continuum_database,graybody +!$OMP THREADPRIVATE(callgasvis,continuum,generic_continuum_database,graybody) logical,save :: strictboundcorrk !$OMP THREADPRIVATE(strictboundcorrk) Index: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 (revision 3640) +++ trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 (revision 3641) @@ -341,4 +341,12 @@ call getin_p("continuum",continuum) if (is_master) write(*,*) trim(rname)//": continuum = ",continuum + + if (is_master) write(*,*) trim(rname)//& + ": use generic continuum opacity database ?"//& + " (matters only if callrad=T)" + generic_continuum_database=.true. ! default value + call getin_p("generic_continuum_database",generic_continuum_database) + if (is_master) write(*,*) trim(rname)//": generic_continuum_database = ", & + generic_continuum_database if (is_master) write(*,*) trim(rname)//": version for H2H2 CIA file ?" Index: trunk/LMDZ.GENERIC/libf/phystd/interpolate_continuum.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/interpolate_continuum.F90 (revision 3641) +++ trunk/LMDZ.GENERIC/libf/phystd/interpolate_continuum.F90 (revision 3641) @@ -0,0 +1,809 @@ + + subroutine interpolate_continuum(filename,igas_X,igas_Y,c_WN,ind_WN,temp,pres_X,pres_Y,abs_coef,firstcall) + +!================================================================== +! +! Purpose +! ------- +! Generic routine to calculate continuum opacities, using lookup tables provided here: https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/ +! More information on the data here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Continuum_Database +! +! Author +! ------- +! M. Turbet (2025) +! +!================================================================== + + use datafile_mod, only: datadir + use mod_phys_lmdz_para, only : is_master + + use gases_h, only: ngasmx, gnom, & + igas_H2, igas_H2O, igas_He, igas_N2, & + igas_CH4, igas_CO2, igas_O2 + + use radinc_h, only: L_NSPECTI, L_NSPECTV + + use radcommon_h, only : BWNV,BWNI,WNOI,WNOV + + + implicit none + + ! input + integer ind_WN ! wavenumber index + integer igas_X ! index of molecule X + integer igas_Y ! index of molecule Y + double precision temp ! temperature (Kelvin) + double precision pres_X ! partial pressure of molecule X (Pascals) + double precision pres_Y ! partial pressure of molecule Y (Pascals) + character*200 filename ! name of the lookup table + character*2 c_WN ! wavelength chanel: infrared (IR) or visible (VI) + logical firstcall + + ! output + double precision abs_coef ! absorption coefficient (m^-1) + + ! intermediate variables + double precision amagat_X ! density of molecule X (in amagat units) + double precision amagat_Y ! density of molecule Y (in amagat units) + + character(len=512) :: line + + integer i, pos, iT, iW, iB, count_norm, igas + + double precision temp_value, temp_abs, temp_wn + + double precision z_temp + + integer num_wn, num_T + + double precision, dimension(:), allocatable :: temp_arr + double precision, dimension(:), allocatable :: wn_arr + double precision, dimension(:,:), allocatable :: abs_arr + + integer ios + + ! Temperature array, continuum absorption grid for the pair N2-N2 + integer num_T_N2N2 + double precision, dimension(:), allocatable :: temp_arr_N2N2 + double precision, dimension(:,:), allocatable :: abs_arr_N2N2_IR + double precision, dimension(:,:), allocatable :: abs_arr_N2N2_VI + save num_T_N2N2,temp_arr_N2N2,abs_arr_N2N2_IR,abs_arr_N2N2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair O2-O2 + integer num_T_O2O2 + double precision, dimension(:), allocatable :: temp_arr_O2O2 + double precision, dimension(:,:), allocatable :: abs_arr_O2O2_IR + double precision, dimension(:,:), allocatable :: abs_arr_O2O2_VI + save num_T_O2O2,temp_arr_O2O2,abs_arr_O2O2_IR,abs_arr_O2O2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair H2-H2 + integer num_T_H2H2 + double precision, dimension(:), allocatable :: temp_arr_H2H2 + double precision, dimension(:,:), allocatable :: abs_arr_H2H2_IR + double precision, dimension(:,:), allocatable :: abs_arr_H2H2_VI + save num_T_H2H2,temp_arr_H2H2,abs_arr_H2H2_IR,abs_arr_H2H2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair CO2-CO2 + integer num_T_CO2CO2 + double precision, dimension(:), allocatable :: temp_arr_CO2CO2 + double precision, dimension(:,:), allocatable :: abs_arr_CO2CO2_IR + double precision, dimension(:,:), allocatable :: abs_arr_CO2CO2_VI + save num_T_CO2CO2,temp_arr_CO2CO2,abs_arr_CO2CO2_IR,abs_arr_CO2CO2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair CH4-CH4 + integer num_T_CH4CH4 + double precision, dimension(:), allocatable :: temp_arr_CH4CH4 + double precision, dimension(:,:), allocatable :: abs_arr_CH4CH4_IR + double precision, dimension(:,:), allocatable :: abs_arr_CH4CH4_VI + save num_T_CH4CH4,temp_arr_CH4CH4,abs_arr_CH4CH4_IR,abs_arr_CH4CH4_VI !read by master + + ! Temperature array, continuum absorption grid for the pair H2O-H2O + integer num_T_H2OH2O + double precision, dimension(:), allocatable :: temp_arr_H2OH2O + double precision, dimension(:,:), allocatable :: abs_arr_H2OH2O_IR + double precision, dimension(:,:), allocatable :: abs_arr_H2OH2O_VI + save num_T_H2OH2O,temp_arr_H2OH2O,abs_arr_H2OH2O_IR,abs_arr_H2OH2O_VI !read by master + + ! Temperature array, continuum absorption grid for the pair H2-He + integer num_T_H2He + double precision, dimension(:), allocatable :: temp_arr_H2He + double precision, dimension(:,:), allocatable :: abs_arr_H2He_IR + double precision, dimension(:,:), allocatable :: abs_arr_H2He_VI + save num_T_H2He,temp_arr_H2He,abs_arr_H2He_IR,abs_arr_H2He_VI !read by master + + ! Temperature array, continuum absorption grid for the pair H2-CH4 + integer num_T_H2CH4 + double precision, dimension(:), allocatable :: temp_arr_H2CH4 + double precision, dimension(:,:), allocatable :: abs_arr_H2CH4_IR + double precision, dimension(:,:), allocatable :: abs_arr_H2CH4_VI + save num_T_H2CH4,temp_arr_H2CH4,abs_arr_H2CH4_IR,abs_arr_H2CH4_VI !read by master + + ! Temperature array, continuum absorption grid for the pair CO2-H2 + integer num_T_CO2H2 + double precision, dimension(:), allocatable :: temp_arr_CO2H2 + double precision, dimension(:,:), allocatable :: abs_arr_CO2H2_IR + double precision, dimension(:,:), allocatable :: abs_arr_CO2H2_VI + save num_T_CO2H2,temp_arr_CO2H2,abs_arr_CO2H2_IR,abs_arr_CO2H2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair CO2-CH4 + integer num_T_CO2CH4 + double precision, dimension(:), allocatable :: temp_arr_CO2CH4 + double precision, dimension(:,:), allocatable :: abs_arr_CO2CH4_IR + double precision, dimension(:,:), allocatable :: abs_arr_CO2CH4_VI + save num_T_CO2CH4,temp_arr_CO2CH4,abs_arr_CO2CH4_IR,abs_arr_CO2CH4_VI !read by master + + ! Temperature array, continuum absorption grid for the pair N2-H2 + integer num_T_N2H2 + double precision, dimension(:), allocatable :: temp_arr_N2H2 + double precision, dimension(:,:), allocatable :: abs_arr_N2H2_IR + double precision, dimension(:,:), allocatable :: abs_arr_N2H2_VI + save num_T_N2H2,temp_arr_N2H2,abs_arr_N2H2_IR,abs_arr_N2H2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair N2-CH4 + integer num_T_N2CH4 + double precision, dimension(:), allocatable :: temp_arr_N2CH4 + double precision, dimension(:,:), allocatable :: abs_arr_N2CH4_IR + double precision, dimension(:,:), allocatable :: abs_arr_N2CH4_VI + save num_T_N2CH4,temp_arr_N2CH4,abs_arr_N2CH4_IR,abs_arr_N2CH4_VI !read by master + + ! Temperature array, continuum absorption grid for the pair CO2-O2 + integer num_T_CO2O2 + double precision, dimension(:), allocatable :: temp_arr_CO2O2 + double precision, dimension(:,:), allocatable :: abs_arr_CO2O2_IR + double precision, dimension(:,:), allocatable :: abs_arr_CO2O2_VI + save num_T_CO2O2,temp_arr_CO2O2,abs_arr_CO2O2_IR,abs_arr_CO2O2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair N2-O2 + integer num_T_N2O2 + double precision, dimension(:), allocatable :: temp_arr_N2O2 + double precision, dimension(:,:), allocatable :: abs_arr_N2O2_IR + double precision, dimension(:,:), allocatable :: abs_arr_N2O2_VI + save num_T_N2O2,temp_arr_N2O2,abs_arr_N2O2_IR,abs_arr_N2O2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair H2O-N2 + integer num_T_H2ON2 + double precision, dimension(:), allocatable :: temp_arr_H2ON2 + double precision, dimension(:,:), allocatable :: abs_arr_H2ON2_IR + double precision, dimension(:,:), allocatable :: abs_arr_H2ON2_VI + save num_T_H2ON2,temp_arr_H2ON2,abs_arr_H2ON2_IR,abs_arr_H2ON2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair H2O-O2 + integer num_T_H2OO2 + double precision, dimension(:), allocatable :: temp_arr_H2OO2 + double precision, dimension(:,:), allocatable :: abs_arr_H2OO2_IR + double precision, dimension(:,:), allocatable :: abs_arr_H2OO2_VI + save num_T_H2OO2,temp_arr_H2OO2,abs_arr_H2OO2_IR,abs_arr_H2OO2_VI !read by master + + ! Temperature array, continuum absorption grid for the pair H2O-CO2 + integer num_T_H2OCO2 + double precision, dimension(:), allocatable :: temp_arr_H2OCO2 + double precision, dimension(:,:), allocatable :: abs_arr_H2OCO2_IR + double precision, dimension(:,:), allocatable :: abs_arr_H2OCO2_VI + save num_T_H2OCO2,temp_arr_H2OCO2,abs_arr_H2OCO2_IR,abs_arr_H2OCO2_VI !read by master + + + + if(firstcall)then ! called by sugas_corrk only + if (is_master) print*,'----------------------------------------------------' + if (is_master) print*,'Initialising continuum (interpolate_continuum routine) from ', trim(filename) + +!$OMP MASTER + + open(unit=33, file=trim(filename), status="old", action="read",iostat=ios) + + if (ios.ne.0) then ! file not found + if (is_master) then + write(*,*) 'Error from interpolate_continuum routine' + write(*,*) 'Data file ',trim(filename),' 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 is there.' + write(*,*) 'Latest continuum data can be downloaded here:' + write(*,*) 'https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/' + endif + call abort + endif + + ! We read the first line of the file to get the number of temperatures provided in the data file + read(33, '(A)') line + + i = 1 + iT = 0 + + do while (i .lt. len_trim(line)) + pos = index(line(i:), 'T=') + if (pos == 0) exit + i = i + pos + iT = iT + 1 + read(line(i+2:i+10), '(E9.2)') temp_value + end do + + num_T=iT ! num_T is the number of temperatures provided in the data file + + ! We read all the remaining lines of the file to get the number of wavenumbers provided in the data file + iW = 0 + do + read(33,*, end=501) line + iW = iW + 1 + end do + +501 continue + + num_wn=iW ! num_wn is the number of wavenumbers provided in the data file + + close(33) + + if(.not.allocated(temp_arr)) allocate(temp_arr(num_T)) + if(.not.allocated(wn_arr)) allocate(wn_arr(num_wn)) + if(.not.allocated(abs_arr)) allocate(abs_arr(num_wn,num_T)) + + ! We now open and read the file a second time to extract the temperature array, wavenumber array and continuum absorption data + + open(unit=33, file=trim(filename), status="old", action="read") + + ! We extract the temperature array (temp_arr) + + read(33, '(A)') line + + i = 1 + iT = 0 + + do while (i .lt. len_trim(line)) + pos = index(line(i:), 'T=') + if (pos == 0) exit + i = i + pos + iT = iT + 1 + read(line(i+2:i+10), '(E9.2)') temp_arr(iT) + end do + + ! We extract the wavenumber array (wn_arr) and continuum absorption (abs_arr) + + do iW=1,num_wn + read(33,*) wn_arr(iW), (abs_arr(iW, iT), iT=1,num_T) + end do + + close(33) + + print*,'We read continuum absorption data for the pair ', trim(gnom(igas_X)),'-',trim(gnom(igas_Y)) + print*,'Temperature grid of the dataset: ', temp_arr(:) + + ! We loop on all molecular pairs with available continuum data and fill the corresponding array + + if ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then + num_T_CO2CO2=num_T + if(.not.allocated(temp_arr_CO2CO2)) allocate(temp_arr_CO2CO2(num_T_CO2CO2)) + if(.not.allocated(abs_arr_CO2CO2_VI)) allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2)) + if(.not.allocated(abs_arr_CO2CO2_IR)) allocate(abs_arr_CO2CO2_IR(L_NSPECTI,num_T_CO2CO2)) + temp_arr_CO2CO2(:)=temp_arr(:) + abs_arr_CO2CO2_VI(:,:)=0. + abs_arr_CO2CO2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CO2_VI,abs_arr_CO2CO2_IR,num_T_CO2CO2) + elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_N2)) then + num_T_N2N2=num_T + if(.not.allocated(temp_arr_N2N2)) allocate(temp_arr_N2N2(num_T_N2N2)) + if(.not.allocated(abs_arr_N2N2_VI)) allocate(abs_arr_N2N2_VI(L_NSPECTV,num_T_N2N2)) + if(.not.allocated(abs_arr_N2N2_IR)) allocate(abs_arr_N2N2_IR(L_NSPECTI,num_T_N2N2)) + temp_arr_N2N2(:)=temp_arr(:) + abs_arr_N2N2_VI(:,:)=0. + abs_arr_N2N2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2N2_VI,abs_arr_N2N2_IR,num_T_N2N2) + elseif ((igas_X .eq. igas_O2) .and. (igas_Y .eq. igas_O2)) then + num_T_O2O2=num_T + if(.not.allocated(temp_arr_O2O2)) allocate(temp_arr_O2O2(num_T_O2O2)) + if(.not.allocated(abs_arr_O2O2_VI)) allocate(abs_arr_O2O2_VI(L_NSPECTV,num_T_O2O2)) + if(.not.allocated(abs_arr_O2O2_IR)) allocate(abs_arr_O2O2_IR(L_NSPECTI,num_T_O2O2)) + temp_arr_O2O2(:)=temp_arr(:) + abs_arr_O2O2_VI(:,:)=0. + abs_arr_O2O2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_O2O2_VI,abs_arr_O2O2_IR,num_T_O2O2) + elseif ((igas_X .eq. igas_CH4) .and. (igas_Y .eq. igas_CH4)) then + num_T_CH4CH4=num_T + if(.not.allocated(temp_arr_CH4CH4)) allocate(temp_arr_CH4CH4(num_T_CH4CH4)) + if(.not.allocated(abs_arr_CH4CH4_VI)) allocate(abs_arr_CH4CH4_VI(L_NSPECTV,num_T_CH4CH4)) + if(.not.allocated(abs_arr_CH4CH4_IR)) allocate(abs_arr_CH4CH4_IR(L_NSPECTI,num_T_CH4CH4)) + temp_arr_CH4CH4(:)=temp_arr(:) + abs_arr_CH4CH4_VI(:,:)=0. + abs_arr_CH4CH4_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CH4CH4_VI,abs_arr_CH4CH4_IR,num_T_CH4CH4) + elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_H2)) then + num_T_H2H2=num_T + if(.not.allocated(temp_arr_H2H2)) allocate(temp_arr_H2H2(num_T_H2H2)) + if(.not.allocated(abs_arr_H2H2_VI)) allocate(abs_arr_H2H2_VI(L_NSPECTV,num_T_H2H2)) + if(.not.allocated(abs_arr_H2H2_IR)) allocate(abs_arr_H2H2_IR(L_NSPECTI,num_T_H2H2)) + temp_arr_H2H2(:)=temp_arr(:) + abs_arr_H2H2_VI(:,:)=0. + abs_arr_H2H2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2H2_VI,abs_arr_H2H2_IR,num_T_H2H2) + elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_H2O)) then + num_T_H2OH2O=num_T + if(.not.allocated(temp_arr_H2OH2O)) allocate(temp_arr_H2OH2O(num_T_H2OH2O)) + if(.not.allocated(abs_arr_H2OH2O_VI)) allocate(abs_arr_H2OH2O_VI(L_NSPECTV,num_T_H2OH2O)) + if(.not.allocated(abs_arr_H2OH2O_IR)) allocate(abs_arr_H2OH2O_IR(L_NSPECTI,num_T_H2OH2O)) + temp_arr_H2OH2O(:)=temp_arr(:) + abs_arr_H2OH2O_VI(:,:)=0. + abs_arr_H2OH2O_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OH2O_VI,abs_arr_H2OH2O_IR,num_T_H2OH2O) + elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_H2)) then + num_T_N2H2=num_T + if(.not.allocated(temp_arr_N2H2)) allocate(temp_arr_N2H2(num_T_N2H2)) + if(.not.allocated(abs_arr_N2H2_VI)) allocate(abs_arr_N2H2_VI(L_NSPECTV,num_T_N2H2)) + if(.not.allocated(abs_arr_N2H2_IR)) allocate(abs_arr_N2H2_IR(L_NSPECTI,num_T_N2H2)) + temp_arr_N2H2(:)=temp_arr(:) + abs_arr_N2H2_VI(:,:)=0. + abs_arr_N2H2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2H2_VI,abs_arr_N2H2_IR,num_T_N2H2) + elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_O2)) then + num_T_N2O2=num_T + if(.not.allocated(temp_arr_N2O2)) allocate(temp_arr_N2O2(num_T_N2O2)) + if(.not.allocated(abs_arr_N2O2_VI)) allocate(abs_arr_N2O2_VI(L_NSPECTV,num_T_N2O2)) + if(.not.allocated(abs_arr_N2O2_IR)) allocate(abs_arr_N2O2_IR(L_NSPECTI,num_T_N2O2)) + temp_arr_N2O2(:)=temp_arr(:) + abs_arr_N2O2_VI(:,:)=0. + abs_arr_N2O2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2O2_VI,abs_arr_N2O2_IR,num_T_N2O2) + elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_CH4)) then + num_T_N2CH4=num_T + if(.not.allocated(temp_arr_N2CH4)) allocate(temp_arr_N2CH4(num_T_N2CH4)) + if(.not.allocated(abs_arr_N2CH4_VI)) allocate(abs_arr_N2CH4_VI(L_NSPECTV,num_T_N2CH4)) + if(.not.allocated(abs_arr_N2CH4_IR)) allocate(abs_arr_N2CH4_IR(L_NSPECTI,num_T_N2CH4)) + temp_arr_N2CH4(:)=temp_arr(:) + abs_arr_N2CH4_VI(:,:)=0. + abs_arr_N2CH4_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2CH4_VI,abs_arr_N2CH4_IR,num_T_N2CH4) + elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_O2)) then + num_T_CO2O2=num_T + if(.not.allocated(temp_arr_CO2O2)) allocate(temp_arr_CO2O2(num_T_CO2O2)) + if(.not.allocated(abs_arr_CO2O2_VI)) allocate(abs_arr_CO2O2_VI(L_NSPECTV,num_T_CO2O2)) + if(.not.allocated(abs_arr_CO2O2_IR)) allocate(abs_arr_CO2O2_IR(L_NSPECTI,num_T_CO2O2)) + temp_arr_CO2O2(:)=temp_arr(:) + abs_arr_CO2O2_VI(:,:)=0. + abs_arr_CO2O2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2O2_VI,abs_arr_CO2O2_IR,num_T_CO2O2) + elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_CH4)) then + num_T_H2CH4=num_T + if(.not.allocated(temp_arr_H2CH4)) allocate(temp_arr_H2CH4(num_T_H2CH4)) + if(.not.allocated(abs_arr_H2CH4_VI)) allocate(abs_arr_H2CH4_VI(L_NSPECTV,num_T_H2CH4)) + if(.not.allocated(abs_arr_H2CH4_IR)) allocate(abs_arr_H2CH4_IR(L_NSPECTI,num_T_H2CH4)) + temp_arr_H2CH4(:)=temp_arr(:) + abs_arr_H2CH4_VI(:,:)=0. + abs_arr_H2CH4_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2CH4_VI,abs_arr_H2CH4_IR,num_T_H2CH4) + elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_He)) then + num_T_H2He=num_T + if(.not.allocated(temp_arr_H2He)) allocate(temp_arr_H2He(num_T_H2He)) + if(.not.allocated(abs_arr_H2He_VI)) allocate(abs_arr_H2He_VI(L_NSPECTV,num_T_H2He)) + if(.not.allocated(abs_arr_H2He_IR)) allocate(abs_arr_H2He_IR(L_NSPECTI,num_T_H2He)) + temp_arr_H2He(:)=temp_arr(:) + abs_arr_H2He_VI(:,:)=0. + abs_arr_H2He_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2He_VI,abs_arr_H2He_IR,num_T_H2He) + elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_N2)) then + num_T_H2ON2=num_T + if(.not.allocated(temp_arr_H2ON2)) allocate(temp_arr_H2ON2(num_T_H2ON2)) + if(.not.allocated(abs_arr_H2ON2_VI)) allocate(abs_arr_H2ON2_VI(L_NSPECTV,num_T_H2ON2)) + if(.not.allocated(abs_arr_H2ON2_IR)) allocate(abs_arr_H2ON2_IR(L_NSPECTI,num_T_H2ON2)) + temp_arr_H2ON2(:)=temp_arr(:) + abs_arr_H2ON2_VI(:,:)=0. + abs_arr_H2ON2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2ON2_VI,abs_arr_H2ON2_IR,num_T_H2ON2) + elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_O2)) then + num_T_H2OO2=num_T + if(.not.allocated(temp_arr_H2OO2)) allocate(temp_arr_H2OO2(num_T_H2OO2)) + if(.not.allocated(abs_arr_H2OO2_VI)) allocate(abs_arr_H2OO2_VI(L_NSPECTV,num_T_H2OO2)) + if(.not.allocated(abs_arr_H2OO2_IR)) allocate(abs_arr_H2OO2_IR(L_NSPECTI,num_T_H2OO2)) + temp_arr_H2OO2(:)=temp_arr(:) + abs_arr_H2OO2_VI(:,:)=0. + abs_arr_H2OO2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OO2_VI,abs_arr_H2OO2_IR,num_T_H2OO2) + elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_CO2)) then + num_T_H2OCO2=num_T + if(.not.allocated(temp_arr_H2OCO2)) allocate(temp_arr_H2OCO2(num_T_H2OCO2)) + if(.not.allocated(abs_arr_H2OCO2_VI)) allocate(abs_arr_H2OCO2_VI(L_NSPECTV,num_T_H2OCO2)) + if(.not.allocated(abs_arr_H2OCO2_IR)) allocate(abs_arr_H2OCO2_IR(L_NSPECTI,num_T_H2OCO2)) + temp_arr_H2OCO2(:)=temp_arr(:) + abs_arr_H2OCO2_VI(:,:)=0. + abs_arr_H2OCO2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OCO2_VI,abs_arr_H2OCO2_IR,num_T_H2OCO2) + elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then + num_T_CO2CO2=num_T + if(.not.allocated(temp_arr_CO2CO2)) allocate(temp_arr_CO2CO2(num_T_CO2CO2)) + if(.not.allocated(abs_arr_CO2CO2_VI)) allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2)) + if(.not.allocated(abs_arr_CO2CO2_IR)) allocate(abs_arr_CO2CO2_IR(L_NSPECTI,num_T_CO2CO2)) + temp_arr_CO2CO2(:)=temp_arr(:) + abs_arr_CO2CO2_VI(:,:)=0. + abs_arr_CO2CO2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CO2_VI,abs_arr_CO2CO2_IR,num_T_CO2CO2) + elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_H2)) then + num_T_CO2H2=num_T + if(.not.allocated(temp_arr_CO2H2)) allocate(temp_arr_CO2H2(num_T_CO2H2)) + if(.not.allocated(abs_arr_CO2H2_VI)) allocate(abs_arr_CO2H2_VI(L_NSPECTV,num_T_CO2H2)) + if(.not.allocated(abs_arr_CO2H2_IR)) allocate(abs_arr_CO2H2_IR(L_NSPECTI,num_T_CO2H2)) + temp_arr_CO2H2(:)=temp_arr(:) + abs_arr_CO2H2_VI(:,:)=0. + abs_arr_CO2H2_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2H2_VI,abs_arr_CO2H2_IR,num_T_CO2H2) + elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CH4)) then + num_T_CO2CH4=num_T + if(.not.allocated(temp_arr_CO2CH4)) allocate(temp_arr_CO2CH4(num_T_CO2CH4)) + if(.not.allocated(abs_arr_CO2CH4_VI)) allocate(abs_arr_CO2CH4_VI(L_NSPECTV,num_T_CO2CH4)) + if(.not.allocated(abs_arr_CO2CH4_IR)) allocate(abs_arr_CO2CH4_IR(L_NSPECTI,num_T_CO2CH4)) + temp_arr_CO2CH4(:)=temp_arr(:) + abs_arr_CO2CH4_VI(:,:)=0. + abs_arr_CO2CH4_IR(:,:)=0. + call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CH4_VI,abs_arr_CO2CH4_IR,num_T_CO2CH4) + endif ! igas_X / igas_Y condition + + +!$OMP END MASTER +!$OMP BARRIER + + + endif ! firstcall + + ! We loop on all molecular pairs with available continuum data and interpolate in the temperature field + ! Two options: we call visible (VI) or infrared (IR) tables, depending on the value of c_WN + + if ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_CO2CO2,num_T_CO2CO2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CO2,num_T_CO2CO2,abs_coef,abs_arr_CO2CO2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CO2,num_T_CO2CO2,abs_coef,abs_arr_CO2CO2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_N2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_N2N2,num_T_N2N2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_N2N2,num_T_N2N2,abs_coef,abs_arr_N2N2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_N2N2,num_T_N2N2,abs_coef,abs_arr_N2N2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_O2) .and. (igas_Y .eq. igas_O2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_O2O2,num_T_O2O2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_O2O2,num_T_O2O2,abs_coef,abs_arr_O2O2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_O2O2,num_T_O2O2,abs_coef,abs_arr_O2O2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_CH4) .and. (igas_Y .eq. igas_CH4)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_CH4CH4,num_T_CH4CH4) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CH4CH4,num_T_CH4CH4,abs_coef,abs_arr_CH4CH4_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CH4CH4,num_T_CH4CH4,abs_coef,abs_arr_CH4CH4_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_H2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_H2H2,num_T_H2H2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2H2,num_T_H2H2,abs_coef,abs_arr_H2H2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2H2,num_T_H2H2,abs_coef,abs_arr_H2H2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_H2O)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_H2OH2O,num_T_H2OH2O) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2OH2O,num_T_H2OH2O,abs_coef,abs_arr_H2OH2O_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2OH2O,num_T_H2OH2O,abs_coef,abs_arr_H2OH2O_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_H2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_N2H2,num_T_N2H2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_N2H2,num_T_N2H2,abs_coef,abs_arr_N2H2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_N2H2,num_T_N2H2,abs_coef,abs_arr_N2H2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_O2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_N2O2,num_T_N2O2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_N2O2,num_T_N2O2,abs_coef,abs_arr_N2O2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_N2O2,num_T_N2O2,abs_coef,abs_arr_N2O2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_CH4)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_N2CH4,num_T_N2CH4) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_N2CH4,num_T_N2CH4,abs_coef,abs_arr_N2CH4_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_N2CH4,num_T_N2CH4,abs_coef,abs_arr_N2CH4_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_O2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_CO2O2,num_T_CO2O2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CO2O2,num_T_CO2O2,abs_coef,abs_arr_CO2O2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CO2O2,num_T_CO2O2,abs_coef,abs_arr_CO2O2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_CH4)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_H2CH4,num_T_H2CH4) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2CH4,num_T_H2CH4,abs_coef,abs_arr_H2CH4_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2CH4,num_T_H2CH4,abs_coef,abs_arr_H2CH4_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_He)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_H2He,num_T_H2He) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2He,num_T_H2He,abs_coef,abs_arr_H2He_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2He,num_T_H2He,abs_coef,abs_arr_H2He_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_N2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_H2ON2,num_T_H2ON2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2ON2,num_T_H2ON2,abs_coef,abs_arr_H2ON2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2ON2,num_T_H2ON2,abs_coef,abs_arr_H2ON2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_O2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_H2OO2,num_T_H2OO2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2OO2,num_T_H2OO2,abs_coef,abs_arr_H2OO2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2OO2,num_T_H2OO2,abs_coef,abs_arr_H2OO2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_CO2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_H2OCO2,num_T_H2OCO2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2OCO2,num_T_H2OCO2,abs_coef,abs_arr_H2OCO2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_H2OCO2,num_T_H2OCO2,abs_coef,abs_arr_H2OCO2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_H2)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_CO2H2,num_T_CO2H2) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CO2H2,num_T_CO2H2,abs_coef,abs_arr_CO2H2_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CO2H2,num_T_CO2H2,abs_coef,abs_arr_CO2H2_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CH4)) then + call T_boundaries_continuum(z_temp,temp,temp_arr_CO2CH4,num_T_CO2CH4) + if(c_WN .eq. 'IR') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CH4,num_T_CO2CH4,abs_coef,abs_arr_CO2CH4_IR(ind_WN,:)) + elseif(c_WN .eq. 'VI') then + call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CH4,num_T_CO2CH4,abs_coef,abs_arr_CO2CH4_VI(ind_WN,:)) + else + print*,'You must select visible (VI) or infrared (IR) canal.' + stop + endif + endif ! igas_X / igas_Y condition + + ! We compute the values of amagat for molecules X and Y + amagat_X = (273.15/temp)*(pres_X/101325.0) + amagat_Y = (273.15/temp)*(pres_Y/101325.0) + + ! We convert the absorption coefficient from cm^-1 amagat^-2 into m^-1 + abs_coef=abs_coef*100.0*amagat_X*amagat_Y + + !print*,'We have ',amagat_X,' amagats of molecule ', trim(gnom(igas_X)) + !print*,'We have ',amagat_X,' amagats of molecule ', trim(gnom(igas_Y)) + !print*,'So the absorption is ',abs_coef,' m^-1' + + return + end subroutine interpolate_continuum + + + subroutine interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr_in,abs_arr_out_VI,abs_arr_out_IR,num_T) + +!================================================================== +! +! Purpose +! ------- +! Interpolate the continuum data into the visible (VI) and infrared (IR) spectral chanels. +! +! Author +! ------- +! M. Turbet (2025) +! +!================================================================== + + use radcommon_h, only : BWNV,BWNI,WNOI,WNOV + use radinc_h, only: L_NSPECTI, L_NSPECTV + use mod_phys_lmdz_para, only : is_master + + implicit none + + integer iW, iB, count_norm + integer num_T + integer num_wn + double precision wn_arr(num_wn) + double precision abs_arr_in(num_wn,num_T) + double precision abs_arr_out_IR(L_NSPECTI,num_T) + double precision abs_arr_out_VI(L_NSPECTV,num_T) + + ! First visible (VI) chanel + + ! We get read of all the wavenumbers lower than the minimum wavenumber in the visible wavenumber grid + iW=1 + do while((wn_arr(iW) .lt. BWNV(1)) .and. (iW .lt. num_wn)) + iW=iW+1 + enddo + + ! We compute the mean of the continuum absorption inside each wavenumber visible (VI) chanel + do iB = 1, L_NSPECTV + count_norm=0 + do while((wn_arr(iW) .lt. BWNV(iB+1)) .and. (iW .lt. num_wn)) + abs_arr_out_VI(iB,:)=abs_arr_out_VI(iB,:)+abs_arr_in(iW,:) + count_norm=count_norm+1 + iW=iW+1 + enddo + if(count_norm .ge. 1) abs_arr_out_VI(iB,:)=abs_arr_out_VI(iB,:)/count_norm + end do + + ! Then infrared (IR) chanel + + ! We get read of all the wavenumbers lower than the minimum wavenumber in the infrared wavenumber grid + iW=1 + do while((wn_arr(iW) .lt. BWNI(1)) .and. (iW .lt. num_wn)) + iW=iW+1 + enddo + + ! We compute the mean of the continuum absorption inside each wavenumber visible (VI) chanel + do iB = 1, L_NSPECTI + count_norm=0 + do while((wn_arr(iW) .lt. BWNI(iB+1)) .and. (iW .lt. num_wn)) + abs_arr_out_IR(iB,:)=abs_arr_out_IR(iB,:)+abs_arr_in(iW,:) + count_norm=count_norm+1 + iW=iW+1 + enddo + if(count_norm .ge. 1) abs_arr_out_IR(iB,:)=abs_arr_out_IR(iB,:)/count_norm + end do + + if (is_master) then + print*, 'Continuum absorption, first temperature, visible (VI):' + do iB = 1, L_NSPECTV + print*,WNOV(iB),' cm-1',abs_arr_out_VI(iB,1), ' cm-1 amagat-2' + end do + + print*, 'Continuum absorption, first temperature, infrared (IR):' + do iB = 1, L_NSPECTI + print*,WNOI(iB),' cm-1',abs_arr_out_IR(iB,1), ' cm-1 amagat-2' + end do + endif + + return + end subroutine interpolate_wn_abs_coeff + + + subroutine T_boundaries_continuum(z_temp,temp,temp_arr,num_T) + +!================================================================== +! +! Purpose +! ------- +! Check if the temperature is outside the boundaries of the continuum data temperatures. +! +! Author +! ------- +! M. Turbet (2025) +! +!================================================================== + + use callkeys_mod, only: strictboundcia + use mod_phys_lmdz_para, only : is_master + + implicit none + + double precision z_temp + double precision temp + integer num_T + double precision temp_arr(num_T) + + z_temp=temp + + if(z_temp .lt. minval(temp_arr)) then + if (strictboundcia) then + if (is_master) then + print*,'Your temperatures are too low for this continuum dataset' + print*, 'Minimum temperature is ', minval(temp_arr), ' K' + endif + stop + else + z_temp=minval(temp_arr) + endif + elseif(z_temp .gt. maxval(temp_arr)) then + if (strictboundcia) then + if (is_master) then + print*,'Your temperatures are too high for this continuum dataset' + print*, 'Maximum temperature is ', maxval(temp_arr), ' K' + endif + stop + else + z_temp=maxval(temp_arr) + endif + endif + + return + end subroutine T_boundaries_continuum + + + subroutine interpolate_T_abs_coeff(z_temp,temp_arr,num_T,abs_coef,abs_arr) + +!================================================================== +! +! Purpose +! ------- +! Interpolate in the continuum data using the temperature field +! +! Author +! ------- +! M. Turbet (2025) +! +!================================================================== + + implicit none + + integer iT + double precision z_temp + integer num_T + double precision temp_arr(num_T) + + double precision abs_coef + double precision abs_arr(num_T) + + ! Check where to interpolate + iT=1 + do while ( z_temp .gt. temp_arr(iT) ) + iT=iT+1 + end do + + ! We proceed to a simple linear interpolation using the two most nearby temperatures + if(iT .lt. num_T) then + abs_coef=abs_arr(iT-1)+(abs_arr(iT)-abs_arr(iT-1))*(z_temp-temp_arr(iT-1))/(temp_arr(iT)-temp_arr(iT-1)) + else + abs_coef=abs_arr(iT) + endif + + !print*,'the absorption is ',abs_coef,' cm^-1 amagat^-2' + + + return + end subroutine interpolate_T_abs_coeff Index: trunk/LMDZ.GENERIC/libf/phystd/optci.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/optci.F90 (revision 3640) +++ trunk/LMDZ.GENERIC/libf/phystd/optci.F90 (revision 3641) @@ -13,7 +13,8 @@ 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, & - igas_CH4, igas_CO2 + igas_CH4, igas_CO2, igas_O2 use comcstfi_mod, only: g, r, mugaz - use callkeys_mod, only: kastprof,continuum,graybody,varspec + use callkeys_mod, only: kastprof,continuum,graybody,varspec, & + generic_continuum_database use recombin_corrk_mod, only: corrk_recombin, gasi_recomb use tpindex_mod, only: tpindex @@ -188,7 +189,77 @@ if(continuum.and.(.not.graybody))then ! include continua if necessary - wn_cont = dble(wnoi(nw)) - T_cont = dble(TMID(k)) - do igas=1,ngasmx + + if(generic_continuum_database)then + T_cont = dble(TMID(k)) + do igas=1,ngasmx + + if(gfrac(igas).eq.-1)then ! variable + p_cont = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol + elseif(varspec) then + p_cont = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k))) + else + p_cont = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k))) + endif + + dtemp=0.0 + + if (igas .eq. igas_N2) then + call interpolate_continuum('',igas_N2,igas_N2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_H2) then + call interpolate_continuum('',igas_N2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_O2) then + call interpolate_continuum('',igas_N2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_CH4) then + call interpolate_continuum('',igas_N2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + elseif (igas .eq. igas_O2) then + call interpolate_continuum('',igas_O2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_CO2) then + call interpolate_continuum('',igas_CO2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + elseif (igas .eq. igas_H2) then + call interpolate_continuum('',igas_H2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_CH4) then + call interpolate_continuum('',igas_H2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_He) then + call interpolate_continuum('',igas_H2,igas_He,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + elseif (igas .eq. igas_CH4) then + call interpolate_continuum('',igas_CH4,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (igas .eq. igas_H2O) then + call interpolate_continuum('',igas_H2O,igas_H2O,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_N2) then + call interpolate_continuum('',igas_H2O,igas_N2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_O2) then + call interpolate_continuum('',igas_H2O,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_CO2) then + call interpolate_continuum('',igas_H2O,igas_CO2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + elseif (igas .eq. igas_CO2) then + call interpolate_continuum('',igas_CO2,igas_CO2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_H2) then + call interpolate_continuum('',igas_CO2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_CH4) then + call interpolate_continuum('',igas_CO2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + endif + + DCONT = DCONT + dtemp + + enddo ! igas=1,ngasmx + else ! generic_continuum_database + wn_cont = dble(wnoi(nw)) + T_cont = dble(TMID(k)) + do igas=1,ngasmx if(gfrac(igas).eq.-1)then ! variable @@ -295,5 +366,7 @@ DCONT = DCONT*dz(k) - endif + endif ! generic_continuum_database + + endif ! continuum do ng=1,L_NGAUSS-1 Index: trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 (revision 3640) +++ trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 (revision 3641) @@ -12,7 +12,8 @@ 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, & - igas_CH4, igas_CO2 + igas_CH4, igas_CO2, igas_O2 use comcstfi_mod, only: g, r, mugaz - use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec + use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec, & + generic_continuum_database use recombin_corrk_mod, only: corrk_recombin, gasv_recomb use tpindex_mod, only: tpindex @@ -193,4 +194,77 @@ if(continuum.and.(.not.graybody).and.callgasvis)then ! include continua if necessary + + if(generic_continuum_database)then + T_cont = dble(TMID(k)) + do igas=1,ngasmx + + if(gfrac(igas).eq.-1)then ! variable + p_cont = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol + elseif(varspec) then + p_cont = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k))) + else + p_cont = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k))) + endif + + dtemp=0.0 + + if (igas .eq. igas_N2) then + call interpolate_continuum('',igas_N2,igas_N2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_H2) then + call interpolate_continuum('',igas_N2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_O2) then + call interpolate_continuum('',igas_N2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_CH4) then + call interpolate_continuum('',igas_N2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + elseif (igas .eq. igas_O2) then + call interpolate_continuum('',igas_O2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_CO2) then + call interpolate_continuum('',igas_CO2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + elseif (igas .eq. igas_H2) then + call interpolate_continuum('',igas_H2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_CH4) then + call interpolate_continuum('',igas_H2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_He) then + call interpolate_continuum('',igas_H2,igas_He,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + elseif (igas .eq. igas_CH4) then + call interpolate_continuum('',igas_CH4,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (igas .eq. igas_H2O) then + call interpolate_continuum('',igas_H2O,igas_H2O,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_N2) then + call interpolate_continuum('',igas_H2O,igas_N2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_O2) then + call interpolate_continuum('',igas_H2O,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_CO2) then + call interpolate_continuum('',igas_H2O,igas_CO2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + elseif (igas .eq. igas_CO2) then + call interpolate_continuum('',igas_CO2,igas_CO2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + do jgas=1,ngasmx + if (jgas .eq. igas_H2) then + call interpolate_continuum('',igas_CO2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + elseif (jgas .eq. igas_CH4) then + call interpolate_continuum('',igas_CO2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo + endif + + DCONT = DCONT + dtemp + + enddo ! igas=1,ngasmx + + else ! generic_continuum_database + + wn_cont = dble(wnov(nw)) T_cont = dble(TMID(k)) @@ -293,6 +367,7 @@ DCONT = DCONT*dz(k) - endif - + endif ! generic_continuum_database + endif ! continuum + do ng=1,L_NGAUSS-1 Index: trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 (revision 3640) +++ trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 (revision 3641) @@ -31,8 +31,8 @@ use gases_h, only: gnom, ngasmx, & igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4, & - igas_CO2 + igas_CO2, igas_O2 use ioipsl_getin_p_mod, only: getin_p use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,& - continuum + continuum,generic_continuum_database use tracer_h, only : nqtot, moderntracdef, is_recomb, noms use recombin_corrk_mod, only: su_recombin, & @@ -716,6 +716,95 @@ ! Initialise the continuum absorption data if(continuum)then - do igas=1,ngasmx - + if(generic_continuum_database)then + do igas=1,ngasmx ! we loop on all pairs of molecules that have data available + ! data can be downloaded from https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/ + if (igas .eq. igas_N2) then + file_id='/continuum_data/' // 'N2-N2_continuum_70-500K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_N2,igas_N2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + do jgas=1,ngasmx + if (jgas .eq. igas_H2) then + file_id='/continuum_data/' // 'H2-N2_continuum_40-600K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_N2,igas_H2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + elseif (jgas .eq. igas_O2) then + file_id='/continuum_data/' // 'O2-N2_continuum_100-500K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_N2,igas_O2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + elseif (jgas .eq. igas_CH4) then + file_id='/continuum_data/' // 'CH4-N2_continuum_40-600K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_N2,igas_CH4,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + endif + enddo + elseif (igas .eq. igas_O2) then + file_id='/continuum_data/' // 'O2-O2_continuum_100-400K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_O2,igas_O2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + do jgas=1,ngasmx + if (jgas .eq. igas_CO2) then + file_id='/continuum_data/' // 'CO2-O2_continuum_150-600K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_CO2,igas_O2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + endif + enddo + elseif (igas .eq. igas_H2) then + file_id='/continuum_data/' // 'H2-H2_continuum_40-7000K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_H2,igas_H2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + do jgas=1,ngasmx + if (jgas .eq. igas_CH4) then + file_id='/continuum_data/' // 'H2-CH4_continuum_40-600K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_H2,igas_CH4,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + elseif (jgas .eq. igas_He) then + file_id='/continuum_data/' // 'H2-He_continuum_40-5500K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_H2,igas_He,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + endif + enddo + elseif (igas .eq. igas_CH4) then + file_id='/continuum_data/' // 'CH4-CH4_continuum_40-500K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_CH4,igas_CH4,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + elseif (igas .eq. igas_H2O) then + file_id='/continuum_data/' // 'H2O-H2O_continuum_100-2000K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_H2O,igas_H2O,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + do jgas=1,ngasmx + if (jgas .eq. igas_N2) then + file_id='/continuum_data/' // 'H2O-N2_continuum_100-2000K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_H2O,igas_N2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + elseif (jgas .eq. igas_O2) then + file_id='/continuum_data/' // 'H2O-O2_continuum_100-2000K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_H2O,igas_O2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + elseif (jgas .eq. igas_CO2) then + file_id='/continuum_data/' // 'H2O-CO2_continuum_100-800K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_H2O,igas_CO2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + endif + enddo + elseif (igas .eq. igas_CO2) then + file_id='/continuum_data/' // 'CO2-CO2_continuum_100-800K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_CO2,igas_CO2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + do jgas=1,ngasmx + if (jgas .eq. igas_H2) then + file_id='/continuum_data/' // 'CO2-H2_continuum_100-800K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_CO2,igas_H2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + elseif (jgas .eq. igas_CH4) then + file_id='/continuum_data/' // 'CO2-CH4_continuum_100-800K_2025.dat' + file_path=TRIM(datadir)//TRIM(file_id) + call interpolate_continuum(file_path,igas_CO2,igas_CH4,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.) + endif + enddo + endif + enddo ! igas=1,ngasmx + + else ! generic_continuum_database tag + do igas=1,ngasmx if (igas .eq. igas_N2) then @@ -771,6 +860,7 @@ endif - enddo - endif + enddo ! igas=1,ngasmx + endif ! generic_continuum_database tag + endif ! continuum flag print*,'----------------------------------------------------'