subroutine mugas_prof(ngrid,nlayer,nq,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.
!
!     Inputs
!     ------
!     ngrid   | Number of atmospheric columns.
!     nlayer  | Number of atmospheric layers.
!     nq      | Number of tracers.
!     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 (2025)
!==================================================================
      
      ! Arguments:
      !~~~~~~~~~~~
      integer,intent(in) :: ngrid
      integer,intent(in) :: nlayer
      integer,intent(in) :: nq
      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, ifine

      real :: dz(ngrid,nlayer)
      real :: zdq_HCN(ngrid,nlayer)
      real :: zdq_C4H2(ngrid,nlayer)
      real :: zdq_C6H6(ngrid,nlayer)
      real :: zdq_C2H2(ngrid,nlayer)
      real :: zdq_C2H6(ngrid,nlayer)
      
      integer Nfine
      parameter(Nfine=999)
      character(len=100) :: file_path
      character(len=100) :: file_name
      real,save,allocatable :: Altdata(:)
      real,save,allocatable :: Pdata(:)
      real,save,allocatable :: Tdata(:)
      real,save,allocatable :: C2H2_fluxdata(:)
      real,save,allocatable :: C2H6_fluxdata(:)
      real,save,allocatable :: C4H2_fluxdata(:)
      real,save,allocatable :: C6H6_fluxdata(:)
      real,save,allocatable :: HCN_fluxdata(:)

      LOGICAL,SAVE :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)

      ! 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

      ! 1. Read data:
      !~~~~~~~~~~~~~~
      if (firstcall) then
        firstcall=.false.

        if (mugasflux_file/='None') then
            file_name = mugasflux_file
        else
            STOP "No filename given for condensable gas fluxes profile. Need to set mugasflux_file..."
        endif
!$OMP MASTER

        if(.not.allocated(Altdata))       allocate(Altdata(Nfine))
        if(.not.allocated(Pdata))         allocate(Pdata(Nfine))
        if(.not.allocated(Tdata))         allocate(Tdata(Nfine))
        if(.not.allocated(C2H2_fluxdata)) allocate(C2H2_fluxdata(Nfine))
        if(.not.allocated(C2H6_fluxdata)) allocate(C2H6_fluxdata(Nfine))
        if(.not.allocated(C4H2_fluxdata)) allocate(C4H2_fluxdata(Nfine))
        if(.not.allocated(C6H6_fluxdata)) allocate(C6H6_fluxdata(Nfine))
        if(.not.allocated(HCN_fluxdata))  allocate(HCN_fluxdata(Nfine))
        
        file_path=file_name
        open(224,file=file_path,form='formatted')
        do ifine = 1, Nfine
            read(224,*) Altdata(ifine), Pdata(ifine), Tdata(ifine), HCN_fluxdata(ifine), &
                        C4H2_fluxdata(ifine), C6H6_fluxdata(ifine), C2H2_fluxdata(ifine), C2H6_fluxdata(ifine)
        enddo
        close(224)
!$OMP END MASTER
!$OMP BARRIER
      endif

      ! 2. Interpolate on the model vertical grid:
      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      do iq = 1, nq
        
        ! C2H2
        if (noms(iq).eq."C2H2_mugas") then
          do ig = 1, ngrid
            call interp_line(Altdata*1000.,C2H2_fluxdata,Nfine,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(Altdata*1000.,C2H6_fluxdata,Nfine,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(Altdata*1000.,C4H2_fluxdata,Nfine,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(Altdata*1000.,C6H6_fluxdata,Nfine,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(Altdata*1000.,HCN_fluxdata,Nfine,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 mugas_prof