Index: trunk/LMDZ.VENUS/libf/phyvenus/interpolate_continuum.F90 =================================================================== --- trunk/LMDZ.VENUS/libf/phyvenus/interpolate_continuum.F90 (revision 3794) +++ trunk/LMDZ.VENUS/libf/phyvenus/interpolate_continuum.F90 (revision 3794) @@ -0,0 +1,833 @@ +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/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/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 + + !!!! TEST !!! + !abs_coef=abs_coef*30 + + !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 infrared (IR) 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 Index: trunk/LMDZ.VENUS/libf/phyvenus/optcv.F90 =================================================================== --- trunk/LMDZ.VENUS/libf/phyvenus/optcv.F90 (revision 3775) +++ trunk/LMDZ.VENUS/libf/phyvenus/optcv.F90 (revision 3794) @@ -11,5 +11,6 @@ use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND use radcommon_h, only: gasv, tlimit, tgasref, pfgasref,wnov,scalep,indv - use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2 + use gases_h, only: gfrac, ngasmx, igas_H2O, igas_CO2 + use interpolate_continuum_mod, only: interpolate_continuum implicit none @@ -72,5 +73,4 @@ real*8 DRAYAER double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc - double precision p_cross ! variable species mixing ratio variables @@ -166,5 +166,4 @@ ! if(continuum)then ! include continua if necessary - wn_cont = dble(wnov(nw)) T_cont = dble(TMID(k)) do igas=1,ngasmx @@ -178,17 +177,12 @@ dtemp=0.0 -! For Venus: only H2O, using CKD - if(igas.eq.igas_H2O.and.T_cont.gt.200.0)then - - p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air! -! if(H2Ocont_simple)then -! call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.) -! else - interm = indv(nw,igas,igas) - call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) - indv(nw,igas,igas) = interm -! endif - - endif +! For Venus: only H2O and CO2 + do jgas=1,ngasmx + if ( ((igas .eq. igas_H2O) .and. (jgas .eq. igas_H2O)) .or. & + ((igas .eq. igas_H2O) .and. (jgas .eq. igas_CO2)) .or. & + ((igas .eq. igas_CO2) .and. (jgas .eq. igas_CO2)) ) then + call interpolate_continuum('',igas,jgas,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.) + endif + enddo DCONT = DCONT + dtemp Index: trunk/LMDZ.VENUS/libf/phyvenus/radcommon_h.F90 =================================================================== --- trunk/LMDZ.VENUS/libf/phyvenus/radcommon_h.F90 (revision 3775) +++ trunk/LMDZ.VENUS/libf/phyvenus/radcommon_h.F90 (revision 3794) @@ -1,4 +1,4 @@ module radcommon_h - use radinc_h, only: L_NSPECTV, NTstart, NTstop, & + use radinc_h, only: L_NSPECTI, L_NSPECTV, NTstart, NTstop, & naerkind, nsizemax implicit none @@ -13,4 +13,13 @@ ! some or all of this common data set ! +! WNOI - Array of wavenumbers at the spectral interval +! centers for the infrared. Array is NSPECTI +! elements long. +! DWNI - Array of "delta wavenumber", i.e., the width, +! in wavenumbers (cm^-1) of each IR spectral +! interval. NSPECTI elements long. +! WAVEI - Array (NSPECTI elements long) of the wavelenght +! (in microns) at the center of each IR spectral +! interval. ! WNOV - Array of wavenumbers at the spectral interval ! center for the VISUAL. Array is NSPECTV @@ -29,4 +38,6 @@ ! TAURAYVAR - Array (NSPECTV elements) of the pressure-independent ! part of Rayleigh scattering optical depth for the variable gas. +! FZEROI - Fraction of zeros in the IR CO2 k-coefficients, for +! each temperature, pressure, and spectral interval ! FZEROV - Fraction of zeros in the VISUAL CO2 k-coefficients, for ! each temperature, pressure, and spectral interval @@ -42,30 +53,51 @@ ! gvis : assymetry factor ! +! Longwave +! ~~~~~~~~ +! +! For the "naerkind" kind of aerosol radiative properties : +! QIRsQREF : Qext / Qext("longrefir") +! omegaIR : mean single scattering albedo +! gIR : mean assymetry factor +! +! +! Note - QIRsQREF in the martian model was scaled to longrefvis, +! here it is scaled to longrefir, which is actually a dummy parameter, +! as we do not output scaled aerosol opacity ... +! + REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) !BWNV read by master in setspv REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV) -!$OMP THREADPRIVATE(WNOV,DWNV,WAVEV,& +!$OMP THREADPRIVATE(WNOI,DWNI,WAVEI,& + !$OMP WNOV,DWNV,WAVEV,& !$OMP STELLARF,TAURAY,TAURAYVAR) + REAL*8 blami(L_NSPECTI+1) REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90 -!$OMP THREADPRIVATE(blamv) +!$OMP THREADPRIVATE(blami,blamv) !! AS: introduced to avoid doing same computations again for continuum + INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indi INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indv -!$OMP THREADPRIVATE(indv) +!$OMP THREADPRIVATE(indi,indv) !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 - REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasv + REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, PFGASREF, GWEIGHT + real*8 FZEROI(L_NSPECTI) real*8 FZEROV(L_NSPECTV) real*8 pgasmin, pgasmax real*8 tgasmin, tgasmax -!$OMP THREADPRIVATE(gasv,& !pgasref,tgasref,pfgasref read by master in sugas_corrk - !$OMP FZEROV) !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk +!$OMP THREADPRIVATE(gasi,gasv,& !pgasref,tgasref,pfgasref read by master in sugas_corrk + !$OMP FZEROI,FZEROV) !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk real QVISsQREF(L_NSPECTV,naerkind,nsizemax) real omegavis(L_NSPECTV,naerkind,nsizemax) real gvis(L_NSPECTV,naerkind,nsizemax) -!$OMP THREADPRIVATE(QVISsQREF,omegavis,gvis) + real QIRsQREF(L_NSPECTV,naerkind,nsizemax) + real omegair(L_NSPECTV,naerkind,nsizemax) + real gir(L_NSPECTV,naerkind,nsizemax) +!$OMP THREADPRIVATE(QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir) @@ -73,5 +105,5 @@ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - REAL lamrefvis(naerkind) + REAL lamrefir(naerkind),lamrefvis(naerkind) ! Actual number of grain size classes in each domain for a @@ -87,16 +119,18 @@ ! Extinction coefficient at reference wavelengths; -! These wavelengths are defined in aeroptproperties, and called longrefvis +! These wavelengths are defined in aeroptproperties, and called longrefvis and longrefir REAL :: QREFvis(naerkind,nsizemax) + REAL :: QREFir(naerkind,nsizemax) + REAL :: omegaREFir(naerkind,nsizemax) REAL,SAVE :: tstellar ! Stellar brightness temperature (SW) - + REAL*8,SAVE :: planckir(naerkind,nsizemax) real*8,save :: PTOP real*8,parameter :: UBARI = 0.5D0 -!$OMP THREADPRIVATE(QREFvis,& ! gweight read by master in sugas_corrk - !$OMP tstellar,PTOP) +!$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFir,& ! gweight read by master in sugas_corrk + !$OMP tstellar,planckir,PTOP) ! If the gas optical depth (top to the surface) is less than Index: trunk/LMDZ.VENUS/libf/phyvenus/radinc_h.F90 =================================================================== --- trunk/LMDZ.VENUS/libf/phyvenus/radinc_h.F90 (revision 3775) +++ trunk/LMDZ.VENUS/libf/phyvenus/radinc_h.F90 (revision 3794) @@ -79,4 +79,5 @@ integer, parameter :: L_NSPECTV = NBvisible + integer, parameter :: L_NSPECTI = NBinfrared integer, parameter :: NAERKIND = 5 Index: trunk/LMDZ.VENUS/libf/phyvenus/sfluxv.F =================================================================== --- trunk/LMDZ.VENUS/libf/phyvenus/sfluxv.F (revision 3775) +++ trunk/LMDZ.VENUS/libf/phyvenus/sfluxv.F (revision 3794) @@ -1,5 +1,5 @@ SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV, * UBAR0,STEL,NFLUXTOPV,FLUXTOPVDN, - * NFLUXOUTV_nu,NFLUXGNDV_nu, + * NFLUXV_nu,NFLUXGNDV_nu, * FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf) @@ -19,5 +19,5 @@ real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD) real*8 NFLUXTOPV, FLUXUP, FLUXDN,FLUXTOPVDN - real*8 NFLUXOUTV_nu(L_NSPECTV) + real*8 NFLUXV_nu(L_NLAYRAD,L_NSPECTV) real*8 NFLUXGNDV_nu(L_NSPECTV) @@ -40,5 +40,5 @@ DO NW=1,L_NSPECTV - NFLUXOUTV_nu(NW)=0.0 + NFLUXV_nu(:,NW)=0.0 NFLUXGNDV_nu(NW)=0.0 END DO @@ -115,7 +115,7 @@ END DO -c band-resolved flux leaving TOA (RDW) - NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW) - * +FLUXUP*GWEIGHT(NG)*(1.0-FZEROV(NW)) +c band-resolved net flux + NFLUXV_nu(:,NW) = NFLUXV_nu(:,NW) + * +(FMUPV-FMDV)*GWEIGHT(NG)*(1.0-FZEROV(NW)) c band-resolved flux at ground (RDW) @@ -173,7 +173,7 @@ END DO -c band-resolved flux leaving TOA (RDW) - NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW) - * +FLUXUP*FZERO +c band-resolved net flux + NFLUXV_nu(:,NW) = NFLUXV_nu(:,NW) + * +(FMUPV-FMDV)*FZERO c band-resolved flux at ground (RDW) Index: trunk/LMDZ.VENUS/libf/phyvenus/su_gases.F90 =================================================================== --- trunk/LMDZ.VENUS/libf/phyvenus/su_gases.F90 (revision 3775) +++ trunk/LMDZ.VENUS/libf/phyvenus/su_gases.F90 (revision 3794) @@ -39,4 +39,12 @@ igas_OCS=5 + ! We don't care about the others + ! but initialized in order to avoid problems in interpolate_continuum.F90 + igas_H2=99 + igas_He=99 + igas_N2=99 + igas_CH4=99 + igas_O2=99 + vgas=0 if(.not.allocated(gfrac)) allocate(gfrac(ngasmx)) Index: trunk/LMDZ.VENUS/libf/phyvenus/sugas_corrk.F90 =================================================================== --- trunk/LMDZ.VENUS/libf/phyvenus/sugas_corrk.F90 (revision 3775) +++ trunk/LMDZ.VENUS/libf/phyvenus/sugas_corrk.F90 (revision 3794) @@ -27,5 +27,6 @@ use radcommon_h, only : WNOV use datafile_mod, only: datadir,corrkdir,banddir - use gases_h, only: gnom, ngasmx, igas_H2O + use gases_h, only: gnom, ngasmx, igas_H2O, igas_CO2 + use interpolate_continuum_mod, only: interpolate_continuum implicit none #include "YOMCST.h" @@ -124,4 +125,6 @@ 'match that in gases.def (',trim(gnom(igas)),'). You should compare ', & 'gases.def with Q.dat in your radiative transfer directory.' + print*,'For Venus, no gases.def, but su_gases.F90 hardcoded... Should be identical ', & + 'to Q.dat in your radiative transfer directory.' message='exiting.' call abort_physic(subname,message,1) @@ -452,22 +455,30 @@ !======================================================================= ! Initialise the continuum absorption data -! if(continuum)then - do igas=1,ngasmx - -! For Venus: only H2O, using CKD - if (igas .eq. igas_H2O) then - - ! H2O is special -! if(H2Ocont_simple)then -! call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.) -! else - dummy = -9999 - call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.,dummy) -! endif - - endif - - enddo -! endif +! if(continuum)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/generic/datagcm/continuum_data/ + +! For Venus: only H2O and CO2 + if (igas .eq. igas_H2O) then + file_id='/continuum_data/' // 'H2O-H2O_2022.cia' + ! H2O-H2O_2022.cia = 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_CO2) then + file_id='/continuum_data/' // 'H2O-CO2_2022.cia' + ! H2O-CO2_2022.cia = 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_2025.cia' + ! CO2-CO2_2025.cia = 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.) + endif + enddo ! igas=1,ngasmx +! endif ! print*,'----------------------------------------------------' Index: trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_corrk.F90 =================================================================== --- trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_corrk.F90 (revision 3775) +++ trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_corrk.F90 (revision 3794) @@ -9,5 +9,5 @@ gweight, pfgasref, pgasmax, pgasmin, & pgasref, tgasmax, tgasmin, tgasref, scalep, & - stellarf, dwnv, tauray + stellarf, dwnv, tauray, wavev, wnov use datafile_mod, only: datadir,banddir,corrkdir use ioipsl_getin_p_mod, only: getin_p @@ -99,5 +99,6 @@ REAL*8 tauaero(L_LEVELS,naerkind) REAL*8 nfluxtopv,fluxtopvdn - REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2). + REAL*8 nfluxv_nu(L_NLAYRAD,L_NSPECTV) ! vertical profil, band-resolved VI net flux (W/m2) + REAL*8 heatrate_nu(L_NLAYRAD,L_NSPECTV) ! vertical profil, band-resolved VI heating rate (K/s/micron) REAL*8 fmnetv(L_NLAYRAD) REAL*8 fluxupv(L_NLAYRAD) @@ -504,5 +505,5 @@ call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, & acosz,stel_fract, & - nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, & + nfluxtopv,fluxtopvdn,nfluxv_nu,nfluxgndv_nu, & fmnetv,fluxupv,fluxdnv,fzerov,taugsurf) @@ -510,5 +511,5 @@ nfluxtopv = 0.0d0 fluxtopvdn = 0.0d0 - nfluxoutv_nu(:) = 0.0d0 + nfluxv_nu(:,:) = 0.0d0 nfluxgndv_nu(:) = 0.0d0 do l=1,L_NLAYRAD @@ -564,4 +565,5 @@ *RG/(cpdet(tmid(1))*scalep*(plevrad(3)-plevrad(2))) + !----------------------------------------------------------------------- !----------------------------------------------------------------------- @@ -570,4 +572,21 @@ !----------------------------------------------------------------------- +! Diagnostic in 1D: +! output the vertical profil of band-resolved heating rates, in K/s/microns + if ((1.eq.1).and.is_master.and.(ngrid.eq.1).and.firstcall) then + open(unit=11,file="nfluxnu.dat") + write(11,*) "lambda(microns)",wavev(:) + write(11,*) "dlbda(microns)",1.e4*dwnv(:)/wnov(:)**2 + write(11,*) "pressure(Pa) dT_nu(K/s/micron)" + DO l=2,L_NLAYRAD + heatrate_nu(L_NLAYRAD+1-l,:)=(nfluxv_nu(l,:)-nfluxv_nu(l-1,:)) & + *RG/(cpdet(tmid(l))*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) & + ! dlbda(microns)=-dnu(cm-1)*1E4/nu(cm-1)^2 + / (1.e4*dwnv(:)/wnov(:)**2) + write(11,*) pplay(1,L_NLAYRAD+1-l),heatrate_nu(L_NLAYRAD+1-l,:) + END DO + close(11) + endif + end subroutine sw_venus_corrk