subroutine evolmugas(ngrid,nlayer,nq,zls,zzlay,zzlev,pplay,pt,pdqchem) use datafile_mod use comcstfi_mod, only: r use tracer_h, only: noms,igcm_C2H2_mugas,igcm_C2H6_mugas,& igcm_C4H2_mugas,igcm_C6H6_mugas,igcm_HCN_mugas use mod_phys_lmdz_para, only : is_master implicit none !================================================================== ! Purpose ! ------- ! Get fixed fluxes of microphysical condensable gas ! from txt file according to the solar longitude ! (for 1D purpose for now). ! ! Inputs ! ------ ! ngrid | Number of atmospheric columns. ! nlayer | Number of atmospheric layers. ! nq | Number of tracers. ! pls | Solar longitude. ! zzlay | Altitude at the middle of the atmospheric layers. ! zzlev | Altitude at the atmospheric layer boundaries. ! pplay | Mid-layer pressure (Pa) ! pt | Temperature (K). ! ! Outputs ! ------- ! ! Both ! ---- ! pdqchem ! Tracer tendencies for condensable gases through muphi (kg/kg_of_air/s). ! ! Authors ! ------- ! Bruno de Batz de Trenquelléon (2026) !================================================================== ! Arguments: !~~~~~~~~~~~ integer,intent(in) :: ngrid integer,intent(in) :: nlayer integer,intent(in) :: nq real,intent(in) :: zls real,intent(in) :: zzlay(ngrid,nlayer) real,intent(in) :: zzlev(ngrid,nlayer+1) real,intent(in) :: pplay(ngrid,nlayer) real,intent(in) :: pt(ngrid,nlayer) real,intent(inout) :: pdqchem(ngrid,nlayer,nq) ! Local variables: !~~~~~~~~~~~~~~~~~ integer :: ig, iq, l integer :: ils, nls, nalt parameter(nls=300) parameter(nalt=600) character(len=1000) :: file_path character(len=100000) :: line real :: diff_min real :: dz(ngrid,nlayer) real :: zdq_C2H2(ngrid,nlayer) real :: zdq_C2H6(ngrid,nlayer) real :: zdq_C4H2(ngrid,nlayer) real :: zdq_C6H6(ngrid,nlayer) real :: zdq_HCN(ngrid,nlayer) real,save,allocatable :: Ls_data(:) real,save,allocatable :: Alt_data(:) real,save,allocatable :: C2H2_fluxls(:) real,save,allocatable :: C2H6_fluxls(:) real,save,allocatable :: C4H2_fluxls(:) real,save,allocatable :: C6H6_fluxls(:) real,save,allocatable :: HCN_fluxls(:) real,save,allocatable :: Flux_C2H2(:,:) real,save,allocatable :: Flux_C2H6(:,:) real,save,allocatable :: Flux_C4H2(:,:) real,save,allocatable :: Flux_C6H6(:,:) real,save,allocatable :: Flux_HCN(:,:) logical,save :: firstcall=.true. ! 0. Initialization: !~~~~~~~~~~~~~~~~~~~ dz(:,:) = zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer) zdq_C2H2(:,:) = 0.d0 zdq_C2H6(:,:) = 0.d0 zdq_C4H2(:,:) = 0.d0 zdq_C6H6(:,:) = 0.d0 zdq_HCN(:,:) = 0.d0 if(.not.allocated(Ls_data)) then allocate(Ls_data(nls)) Ls_data(:) = 0.d0 endif if(.not.allocated(Alt_data)) then allocate(Alt_data(nalt)) Alt_data(:) = 0.d0 endif if(.not.allocated(C2H2_fluxls)) then allocate(C2H2_fluxls(nalt)) C2H2_fluxls(:) = 0.d0 endif if(.not.allocated(Flux_C2H2)) then allocate(Flux_C2H2(nls,nalt)) Flux_C2H2(:,:) = 0.d0 endif if(.not.allocated(C2H6_fluxls)) then allocate(C2H6_fluxls(nalt)) C2H6_fluxls(:) = 0.d0 endif if(.not.allocated(Flux_C2H6)) then allocate(Flux_C2H6(nls,nalt)) Flux_C2H6(:,:) = 0.d0 endif if(.not.allocated(C4H2_fluxls)) then allocate(C4H2_fluxls(nalt)) C4H2_fluxls(:) = 0.d0 endif if(.not.allocated(Flux_C4H2)) then allocate(Flux_C4H2(nls,nalt)) Flux_C4H2(:,:) = 0.d0 endif if(.not.allocated(C6H6_fluxls)) then allocate(C6H6_fluxls(nalt)) C6H6_fluxls(:) = 0.d0 endif if(.not.allocated(Flux_C6H6)) then allocate(Flux_C6H6(nls,nalt)) Flux_C6H6(:,:) = 0.d0 endif if(.not.allocated(HCN_fluxls)) then allocate(HCN_fluxls(nalt)) HCN_fluxls(:) = 0.d0 endif if(.not.allocated(Flux_HCN)) then allocate(Flux_HCN(nls,nalt)) Flux_HCN(:,:) = 0.d0 endif ! 1. Read data: !~~~~~~~~~~~~~~ IF (firstcall) then firstcall=.false. ! C2H2 !----- file_path=trim(datadir)//'/gas_prop/c2h2_vs_ls.txt' open(unit=288, file=file_path, status="old") read(288,*) Ls_data read(288,*) Alt_data do ils = 1, nls read(288,*) Flux_C2H2(ils,:) enddo close(288) ! C2H6 !----- file_path=trim(datadir)//'/gas_prop/c2h6_vs_ls.txt' open(unit=289, file=file_path, status="old") read(289,*) read(289,*) do ils = 1, nls read(289,*) Flux_C2H6(ils,:) enddo close(289) ! C4H2 !----- file_path=trim(datadir)//'/gas_prop/c4h2_vs_ls.txt' open(unit=290, file=file_path, status="old") read(290,*) read(290,*) do ils = 1, nls read(290,*) Flux_C4H2(ils,:) enddo close(290) ! C6H6 !----- file_path=trim(datadir)//'/gas_prop/c6h6_vs_ls.txt' open(unit=291, file=file_path, status="old") read(291,*) read(291,*) do ils = 1, nls read(291,*) Flux_C6H6(ils,:) enddo close(291) ! HCN !----- file_path=trim(datadir)//'/gas_prop/hcn_vs_ls.txt' open(unit=292, file=file_path, status="old") read(292,*) read(292,*) do ils = 1, nls read(292,*) Flux_HCN(ils,:) enddo close(292) ENDIF ! 2. Interpolate on the model vertical grid: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff_min = abs(Ls_data(1) - zls) C2H2_fluxls(:) = Flux_C2H2(1,:) C2H6_fluxls(:) = Flux_C2H6(1,:) C4H2_fluxls(:) = Flux_C4H2(1,:) C6H6_fluxls(:) = Flux_C6H6(1,:) HCN_fluxls(:) = Flux_HCN(1,:) do ils = 2, nls if (abs(Ls_data(ils) - zls) < diff_min) then diff_min = abs(Ls_data(ils) - zls) C2H2_fluxls(:) = Flux_C2H2(ils,:) C2H6_fluxls(:) = Flux_C2H6(ils,:) C4H2_fluxls(:) = Flux_C4H2(ils,:) C6H6_fluxls(:) = Flux_C6H6(ils,:) HCN_fluxls(:) = Flux_HCN(ils,:) end if end do do iq = 1, nq ! C2H2 !----- if (noms(iq).eq."C2H2_mugas") then do ig = 1, ngrid call interp_line(Alt_data*1000.,C2H2_fluxls,nalt,zzlay(ig,:),zdq_C2H2(ig,:),nlayer) do l = 1, nlayer pdqchem(ig,l,igcm_C2H2_mugas) = zdq_C2H2(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) enddo enddo endif ! C2H6 !----- if (noms(iq).eq."C2H6_mugas") then do ig = 1, ngrid call interp_line(Alt_data*1000.,C2H6_fluxls,nalt,zzlay(ig,:),zdq_C2H6(ig,:),nlayer) do l = 1, nlayer pdqchem(ig,l,igcm_C2H6_mugas) = zdq_C2H6(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) enddo enddo endif ! C4H2 !----- if (noms(iq).eq."C4H2_mugas") then do ig = 1, ngrid call interp_line(Alt_data*1000.,C4H2_fluxls,nalt,zzlay(ig,:),zdq_C4H2(ig,:),nlayer) do l = 1, nlayer pdqchem(ig,l,igcm_C4H2_mugas) = zdq_C4H2(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) enddo enddo endif ! C6H6 !----- if (noms(iq).eq."C6H6_mugas") then do ig = 1, ngrid call interp_line(Alt_data*1000.,C6H6_fluxls,nalt,zzlay(ig,:),zdq_C6H6(ig,:),nlayer) do l = 1, nlayer pdqchem(ig,l,igcm_C6H6_mugas) = zdq_C6H6(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) enddo enddo endif ! HCN !---- if (noms(iq).eq."HCN_mugas") then do ig = 1, ngrid call interp_line(Alt_data*1000.,HCN_fluxls,nalt,zzlay(ig,:),zdq_HCN(ig,:),nlayer) do l = 1, nlayer pdqchem(ig,l,igcm_HCN_mugas) = zdq_HCN(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) enddo enddo endif enddo ! end of nq end subroutine evolmugas