module interpolate_continuum_mod implicit none contains 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,intent(in) :: ind_WN ! wavenumber index integer,intent(in) :: igas_X ! index of molecule X integer,intent(in) :: igas_Y ! index of molecule Y double precision,intent(in) :: temp ! temperature (Kelvin) double precision,intent(in) :: pres_X ! partial pressure of molecule X (Pascals) double precision,intent(in) :: pres_Y ! partial pressure of molecule Y (Pascals) character(len=*),intent(in) :: filename ! name of the lookup table character(len=2),intent(in) :: c_WN ! wavelength chanel: infrared (IR) or visible (VI) logical,intent(in) :: firstcall ! output double precision,intent(out) :: 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,save :: num_T_N2N2 double precision,save,dimension(:),allocatable :: temp_arr_N2N2 double precision,save,dimension(:,:),allocatable :: abs_arr_N2N2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_N2N2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair O2-O2 integer,save :: num_T_O2O2 double precision,save,dimension(:),allocatable :: temp_arr_O2O2 double precision,save,dimension(:,:),allocatable :: abs_arr_O2O2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_O2O2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair H2-H2 integer,save :: num_T_H2H2 double precision,save,dimension(:),allocatable :: temp_arr_H2H2 double precision,save,dimension(:,:),allocatable :: abs_arr_H2H2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_H2H2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair CO2-CO2 integer,save :: num_T_CO2CO2 double precision,save,dimension(:),allocatable :: temp_arr_CO2CO2 double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CO2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CO2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair CH4-CH4 integer,save :: num_T_CH4CH4 double precision,save,dimension(:),allocatable :: temp_arr_CH4CH4 double precision,save,dimension(:,:),allocatable :: abs_arr_CH4CH4_IR double precision,save,dimension(:,:),allocatable :: abs_arr_CH4CH4_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair H2O-H2O integer,save :: num_T_H2OH2O double precision,save,dimension(:),allocatable :: temp_arr_H2OH2O double precision,save,dimension(:,:),allocatable :: abs_arr_H2OH2O_IR double precision,save,dimension(:,:),allocatable :: abs_arr_H2OH2O_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair H2-He integer,save :: num_T_H2He double precision,save,dimension(:),allocatable :: temp_arr_H2He double precision,save,dimension(:,:),allocatable :: abs_arr_H2He_IR double precision,save,dimension(:,:),allocatable :: abs_arr_H2He_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair H2-CH4 integer,save :: num_T_H2CH4 double precision,save,dimension(:),allocatable :: temp_arr_H2CH4 double precision,save,dimension(:,:),allocatable :: abs_arr_H2CH4_IR double precision,save,dimension(:,:),allocatable :: abs_arr_H2CH4_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair CO2-H2 integer,save :: num_T_CO2H2 double precision,save,dimension(:),allocatable :: temp_arr_CO2H2 double precision,save,dimension(:,:),allocatable :: abs_arr_CO2H2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_CO2H2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair CO2-CH4 integer,save :: num_T_CO2CH4 double precision,save,dimension(:),allocatable :: temp_arr_CO2CH4 double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CH4_IR double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CH4_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair N2-H2 integer,save :: num_T_N2H2 double precision,save,dimension(:),allocatable :: temp_arr_N2H2 double precision,save,dimension(:,:),allocatable :: abs_arr_N2H2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_N2H2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair N2-CH4 integer,save :: num_T_N2CH4 double precision,save,dimension(:),allocatable :: temp_arr_N2CH4 double precision,save,dimension(:,:),allocatable :: abs_arr_N2CH4_IR double precision,save,dimension(:,:),allocatable :: abs_arr_N2CH4_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair CO2-O2 integer,save :: num_T_CO2O2 double precision,save,dimension(:),allocatable :: temp_arr_CO2O2 double precision,save,dimension(:,:),allocatable :: abs_arr_CO2O2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_CO2O2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair N2-O2 integer,save :: num_T_N2O2 double precision,save,dimension(:), allocatable :: temp_arr_N2O2 double precision,save,dimension(:,:), allocatable :: abs_arr_N2O2_IR double precision,save,dimension(:,:), allocatable :: abs_arr_N2O2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair H2O-N2 integer,save :: num_T_H2ON2 double precision,save,dimension(:),allocatable :: temp_arr_H2ON2 double precision,save,dimension(:,:),allocatable :: abs_arr_H2ON2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_H2ON2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair H2O-O2 integer,save :: num_T_H2OO2 double precision,save,dimension(:),allocatable :: temp_arr_H2OO2 double precision,save,dimension(:,:),allocatable :: abs_arr_H2OO2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_H2OO2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared ! Temperature array, continuum absorption grid for the pair H2O-CO2 integer,save :: num_T_H2OCO2 double precision,save,dimension(:),allocatable :: temp_arr_H2OCO2 double precision,save,dimension(:,:),allocatable :: abs_arr_H2OCO2_IR double precision,save,dimension(:,:),allocatable :: abs_arr_H2OCO2_VI ! None of these saved variables are THREADPRIVATE because read by master ! and then only accessed but never modified and thus can be shared 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_physic("interpolate_continuum","missing input file",1) 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) allocate(temp_arr(num_T)) allocate(wn_arr(num_wn)) 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 allocate(temp_arr_CO2CO2(num_T_CO2CO2)) allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2)) 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 allocate(temp_arr_N2N2(num_T_N2N2)) allocate(abs_arr_N2N2_VI(L_NSPECTV,num_T_N2N2)) 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 allocate(temp_arr_O2O2(num_T_O2O2)) allocate(abs_arr_O2O2_VI(L_NSPECTV,num_T_O2O2)) 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 allocate(temp_arr_CH4CH4(num_T_CH4CH4)) allocate(abs_arr_CH4CH4_VI(L_NSPECTV,num_T_CH4CH4)) 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 allocate(temp_arr_H2H2(num_T_H2H2)) allocate(abs_arr_H2H2_VI(L_NSPECTV,num_T_H2H2)) 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 allocate(temp_arr_H2OH2O(num_T_H2OH2O)) allocate(abs_arr_H2OH2O_VI(L_NSPECTV,num_T_H2OH2O)) 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 allocate(temp_arr_N2H2(num_T_N2H2)) allocate(abs_arr_N2H2_VI(L_NSPECTV,num_T_N2H2)) 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 allocate(temp_arr_N2O2(num_T_N2O2)) allocate(abs_arr_N2O2_VI(L_NSPECTV,num_T_N2O2)) 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 allocate(temp_arr_N2CH4(num_T_N2CH4)) allocate(abs_arr_N2CH4_VI(L_NSPECTV,num_T_N2CH4)) 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 allocate(temp_arr_CO2O2(num_T_CO2O2)) allocate(abs_arr_CO2O2_VI(L_NSPECTV,num_T_CO2O2)) 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 allocate(temp_arr_H2CH4(num_T_H2CH4)) allocate(abs_arr_H2CH4_VI(L_NSPECTV,num_T_H2CH4)) 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 allocate(temp_arr_H2He(num_T_H2He)) allocate(abs_arr_H2He_VI(L_NSPECTV,num_T_H2He)) 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 allocate(temp_arr_H2ON2(num_T_H2ON2)) allocate(abs_arr_H2ON2_VI(L_NSPECTV,num_T_H2ON2)) 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 allocate(temp_arr_H2OO2(num_T_H2OO2)) allocate(abs_arr_H2OO2_VI(L_NSPECTV,num_T_H2OO2)) 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 allocate(temp_arr_H2OCO2(num_T_H2OCO2)) allocate(abs_arr_H2OCO2_VI(L_NSPECTV,num_T_H2OCO2)) 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 allocate(temp_arr_CO2CO2(num_T_CO2CO2)) allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2)) 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 allocate(temp_arr_CO2H2(num_T_CO2H2)) allocate(abs_arr_CO2H2_VI(L_NSPECTV,num_T_CO2H2)) 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 allocate(temp_arr_CO2CH4(num_T_CO2CH4)) allocate(abs_arr_CO2CH4_VI(L_NSPECTV,num_T_CO2CH4)) 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' 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,intent(in) :: num_T integer,intent(in) :: num_wn double precision,intent(in) :: wn_arr(num_wn) double precision,intent(in) :: abs_arr_in(num_wn,num_T) double precision,intent(out) :: abs_arr_out_IR(L_NSPECTI,num_T) double precision,intent(out) :: 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 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,intent(out) :: z_temp double precision,intent(in) :: temp integer,intent(in) :: num_T double precision,intent(in) :: temp_arr(num_T) character(len=22) :: rname = "T_boundaries_continuum" 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 call abort_physic(rname,"temperature too low",1) 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 call abort_physic(rname,"temperature too low",1) else z_temp=maxval(temp_arr) endif endif 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' end subroutine interpolate_T_abs_coeff end module interpolate_continuum_mod