Index: trunk/LMDZ.PLUTO/libf/phypluto/evolmugas.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/evolmugas.F90 (revision 4133) +++ trunk/LMDZ.PLUTO/libf/phypluto/evolmugas.F90 (revision 4133) @@ -0,0 +1,301 @@ +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