SUBROUTINE initracer(ngrid,nq)

      use surfdat_h, ONLY: dryness
      USE tracer_h
      USE callkeys_mod, only: optichaze,nb_monomer,haze,fractal,fasthaze,rad_haze,callmufi
      USE recombin_corrk_mod, ONLY: ini_recombin
      USE mod_phys_lmdz_para, only: is_master, bcast
      use aerosol_mod, only: iaero_haze,i_haze
      IMPLICIT NONE
!=======================================================================
!   subject:
!   --------
!   Initialization related to tracer
!   (transported dust, water, chemical species, ice...)
!
!   Name of the tracer
!
!   Test of dimension :
!   Initialize COMMON tracer in tracer.h, using tracer names provided
!   by the argument nametrac
!
!   author: F.Forget
!   ------
!            Ehouarn Millour (oct. 2008): identify tracers by their names
!            Y. Jaziri & J. Vatant d'Ollone (2020) : modern traceur.def
!            B. de Batz de Trenquelléon (2024): specific microphysical tracers   
!=======================================================================

      integer,intent(in) :: ngrid,nq

      character(len=500) :: tracline   ! to read traceur.def lines
      integer :: blank      !to store the index of 1st blank when reading tracers names
      integer iq,ig,count,ierr
      real r0_lift , reff_lift, rho_haze
      integer nqhaze(nq)               ! to store haze tracers
      integer i, ia, block, j
      character(len=20) :: txt ! to store some text
      character(LEN=20) :: tracername ! to temporarily store text
      character(LEN=20) :: str

!-----------------------------------------------------------------------
!  radius(nq)      ! aerosol particle radius (m)
!  rho_q(nq)       ! tracer densities (kg.m-3)
!  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
!  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
!  alpha_devil(nq) ! lifting coeeficient by dust devil
!  rho_dust          ! Mars dust density
!  rho_ice           ! Water ice density
!  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
!  varian            ! Characteristic variance of log-normal distribution
!-----------------------------------------------------------------------

      if (is_master) then ! only the master proc/thread needs do this

        moderntracdef=.false. ! For modern traceur.def (default false, old type)

        open(407, form = 'formatted', status = 'old', &
             file = 'traceur.def', iostat=ierr)
        if (ierr /=0) then
          ! call abort_physic('initracer',&
          ! 'Problem in opening traceur.def',1)
          print*,'Problem in opening traceur.def'
          stop
        end if
!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
        READ(407,'(A)') tracline
        IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
          READ(tracline,*) nqtot ! Try standard traceur.def
        ELSE
         moderntracdef = .true.
         DO
           READ(407,'(A)',iostat=ierr) tracline
           IF (ierr==0) THEN
             IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
               READ(tracline,*) nqtot
               ! Temporary not implemented solution
               if (nqtot/=nq) then
      !            call abort_physic('initracer','Different number of '// &
      ! 'tracers in dynamics and physics not managed yet',1)
                print*,'!= nbr oftracers in dynamics and physics not managed yet'
                stop
               endif
               EXIT
             ENDIF
           ELSE ! If pb, or if reached EOF without having found nqtot
      !        call abort_physic('initracer','Unable to read numbers '// &
      !  'of tracers in traceur.def',1)
            print*,"unable to read numbers of tracer in tracer.def"
            stop
           ENDIF
         ENDDO
        ENDIF ! if modern or standard traceur.def

       endif ! of if (is_master)

       ! share the information with other procs/threads (if any)
       CALL bcast(nqtot)
       CALL bcast(moderntracdef)

!! -----------------------------------------------------------------------
       !! For the moment number of tracers in dynamics and physics are equal
       nqtot=nq
       !! we allocate once for all arrays in common in tracer_h.F90
       !! (supposedly those are not used before call to initracer)
       IF (.NOT.ALLOCATED(noms))           ALLOCATE(noms(nq))
       IF (.NOT.ALLOCATED(mmol))           ALLOCATE(mmol(nq))
       IF (.NOT.ALLOCATED(aki))            ALLOCATE(aki(nqtot))
       IF (.NOT.ALLOCATED(cpi))            ALLOCATE(cpi(nqtot))
       IF (.NOT.ALLOCATED(radius))         ALLOCATE(radius(nq))
       IF (.NOT.ALLOCATED(rho_q))          ALLOCATE(rho_q(nq))
       IF (.NOT.ALLOCATED(qext))           ALLOCATE(qext(nq))
       IF (.NOT.ALLOCATED(alpha_lift))     ALLOCATE(alpha_lift(nq))
       IF (.NOT.ALLOCATED(alpha_devil))    ALLOCATE(alpha_devil(nq))
       IF (.NOT.ALLOCATED(qextrhor))       ALLOCATE(qextrhor(nq))
      !  IF (.NOT.ALLOCATED(igcm_dustbin))   ALLOCATE(igcm_dustbin(nq))
       IF (.NOT.ALLOCATED(is_chim))        ALLOCATE(is_chim(nqtot))
       IF (.NOT.ALLOCATED(is_rad))         ALLOCATE(is_rad(nqtot))
       IF (.NOT.ALLOCATED(is_recomb))      ALLOCATE(is_recomb(nqtot))
       IF (.NOT.ALLOCATED(is_recomb_qset)) THEN
         ALLOCATE(is_recomb_qset(nqtot))
       ENDIF
       IF (.NOT.ALLOCATED(is_recomb_qotf)) THEN
         ALLOCATE(is_recomb_qotf(nqtot))
       ENDIF

       !! initialization
       alpha_lift(:)     = 0.
       alpha_devil(:)    = 0.
       mmol(:)           = 0.
       aki(:)            = 0.
       cpi(:)            = 0.
       is_chim(:)        = 0
       is_rad(:)         = 0
       is_recomb(:)      = 0
       is_recomb_qset(:) = 0
       is_recomb_qotf(:) = 0

       ! Added by JVO 2017 : these arrays are handled later
       ! -> initialization is the least we can do, please !!!
       radius(:)=0.
       qext(:)=0.

       rho_q(:) = 0. !need to be init here if we want to read it from modern traceur with get_tracdat


! Initialization: Read tracers names from traceur.def
        do iq=1,nq
          if (is_master) read(407,'(A)') tracline
          call bcast(tracline)
          blank = index(tracline,' ') ! Find position of 1st blank in tracline
          noms(iq) = tracline(1:blank) !ensure that in Modern-trac case, noms is actually the name and not all properties
        enddo

! Identify tracers by their names: (and set corresponding values of mmol)
      ! 0. initialize tracer indexes to zero:
      ! NB: igcm_* indexes are commons in 'tracer.h'
      igcm_n2=0
      igcm_ch4_gas=0
      igcm_ch4_ice=0
      igcm_prec_haze=0
      igcm_co_gas=0
      igcm_co_ice=0

      nqhaze(:)=0
      i=0
      DO iq=1,nq
         txt=noms(iq)
         IF (txt(1:4).eq."haze") THEN
            i=i+1
            nqhaze(i)=iq
         ENDIF
      ENDDO
      if ((haze.or.fasthaze).and.i==0) then
         print*, 'Haze active but no haze tracer in traceur.def'
         stop
      endif

      igcm_haze=0
      igcm_haze10=0
      igcm_haze30=0
      igcm_haze50=0
      igcm_haze100=0

!     Eddy diffusion tracers
      igcm_eddy1e6=0
      igcm_eddy1e7=0
      igcm_eddy5e7=0
      igcm_eddy1e8=0
      igcm_eddy5e8=0
      write(*,*) 'initracer: noms() ', noms
      rho_n2=1000        ! n2 ice
      rho_ch4_ice=520.       ! rho ch4 ice
      rho_co_ice=520.       ! rho ch4 ice
      rho_haze=800.     ! haze

      ! 1. find dust tracers
      count=0

      ! 2. find chemistry and water tracers
      do iq=1,nq
        if (noms(iq).eq."n2") then
          igcm_n2=iq
          mmol(igcm_n2)=28.
          count=count+1
          write(*,*) 'Tracer ',count,' = n2'
        endif
        if (noms(iq).eq."ch4_gas") then
          igcm_ch4_gas=iq
          mmol(igcm_ch4_gas)=16.
          count=count+1
          write(*,*) 'Tracer ',count,' = ch4 gas'
        endif
        if (noms(iq).eq."ch4_ice") then
          igcm_ch4_ice=iq
          mmol(igcm_ch4_ice)=16.
          radius(iq)=3.e-6
          rho_q(iq)=rho_ch4_ice
          count=count+1
          write(*,*) 'Tracer ',count,' = ch4 ice'
        endif
        if (noms(iq).eq."co_gas") then
          igcm_co_gas=iq
          mmol(igcm_co_gas)=28.
          count=count+1
          write(*,*) 'Tracer ',count,' = co gas'
        endif
        if (noms(iq).eq."co_ice") then
          igcm_co_ice=iq
          mmol(igcm_co_ice)=28.
          radius(iq)=3.e-6
          rho_q(iq)=rho_co_ice
          count=count+1
          write(*,*) 'Tracer ',count,' = co ice'
        endif
        if (noms(iq).eq."prec_haze") then
          igcm_prec_haze=iq
          count=count+1
          write(*,*) 'Tracer ',count,' = prec haze'
        endif
        if (noms(iq).eq."haze") then
          igcm_haze=iq
          count=count+1
          radius(iq)=rad_haze
          rho_q(iq)=rho_haze
          write(*,*) 'Tracer ',count,' = haze'
        endif
        if (noms(iq).eq."haze10") then
          igcm_haze10=iq
          count=count+1
          radius(iq)=10.e-9
          rho_q(iq)=rho_haze
          write(*,*) 'Tracer ',count,' = haze10'
        endif
        if (noms(iq).eq."haze30") then
          igcm_haze30=iq
          count=count+1
          radius(iq)=30.e-9
          rho_q(iq)=rho_haze
          write(*,*) 'Tracer ',count,' = haze30'
        endif
        if (noms(iq).eq."haze50") then
          igcm_haze50=iq
          count=count+1
          radius(iq)=50.e-9
          rho_q(iq)=rho_haze
          write(*,*) 'Tracer ',count,' = haze50'
        endif
        if (noms(iq).eq."haze100") then
          igcm_haze100=iq
          count=count+1
          radius(iq)=100.e-9
          rho_q(iq)=rho_haze
          write(*,*) 'Tracer ',count,' = haze100'
        endif
!       Eddy diffusion tracers
        if (noms(iq).eq."eddy1e6") then
          igcm_eddy1e6=iq
          count=count+1
          write(*,*) 'Tracer ',count,' = eddy1e6'
        endif
        if (noms(iq).eq."eddy1e7") then
          igcm_eddy1e7=iq
          count=count+1
          write(*,*) 'Tracer ',count,' = eddy1e7'
        endif
        if (noms(iq).eq."eddy5e7") then
          igcm_eddy5e7=iq
          count=count+1
          write(*,*) 'Tracer ',count,' = eddy5e7'
        endif
        if (noms(iq).eq."eddy1e8") then
          igcm_eddy1e8=iq
          count=count+1
          write(*,*) 'Tracer ',count,' = eddy1e8'
        endif
        if (noms(iq).eq."eddy5e8") then
          igcm_eddy5e8=iq
          count=count+1
          write(*,*) 'Tracer ',count,' = eddy5e8'
        endif
      enddo ! of do iq=1,nq

      ! ! 3. find condensable traceurs different from h2o and n2
      ! do iq=1,nq
      !   if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
      !     count=count+1
      !   endif
      !   if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
      !     count=count+1
      !   endif

      ! enddo ! of do iq=1,nq

      ! check that we identified all tracers:
      if (count.ne.nq) then
        write(*,*) "initracer: found only ",count," tracers"
        write(*,*) "               expected ",nq
        do iq=1,count
          write(*,*)'      ',iq,' ',trim(noms(iq))
        enddo
!        stop
      else
        write(*,*) "initracer: found all expected tracers, namely:"
        do iq=1,nq
          write(*,*)'      ',iq,' ',trim(noms(iq))
        enddo
      endif
      
      ! Compute number of microphysics tracers:
      ! By convention they all have the prefix "mu_" (case sensitive !)
      nmicro = 0
      IF (callmufi) THEN
         DO iq=1,nq
            str = noms(iq)
            IF (str(1:3) == "mu_") THEN
               nmicro = nmicro+1
               count = count+1
            ENDIF
         ENDDO

         ! Checking the expected number of tracers:
         ! Microphysics moment model: nmicro = 4
         IF (nmicro < 4) THEN
            WRITE(*,*) "initracer:error:"," Inconsistent number of microphysical tracers"
            WRITE(*,*) "expected at least 4 tracers,", nmicro, " given"
            CALL abort
         ELSE IF (nmicro > 4) THEN
            WRITE(*,*) "!!! WARNING !!! initracer: I was expecting only four tracers, you gave me more."
            CALL abort
         ENDIF

         ! microphysics indexes share the same values than original tracname.
         IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro))
         j = 1
         DO i=1,nq
            str = noms(i)
            IF (str(1:3) == "mu_") THEN
               micro_indx(j) = i 
               j=j+1
            ENDIF
         ENDDO
      
      ELSE
         IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro))
      
      ENDIF ! end of callmufi

      ! Get data of tracers. Need to rewind traceur.def first
      if (is_master) then
       rewind(407)
       do
        read(407,'(A)') tracline
        if (index(tracline,"#") .ne. 1) then
          exit
        endif
       enddo
      endif
      do iq=1,nqtot
        if (is_master) read(407,'(A)') tracline
        call bcast(tracline)
        call get_tracdat(iq, tracline)
      enddo

      if (is_master) close(407)

      ! Calculate number of species in the chemistry
      nesp = sum(is_chim)
      write(*,*) 'Number of species in the chemistry nesp = ',nesp

      ! Calculate number of microphysical tracer
      write(*,*) 'Number of microphysical tracer nmicro = ',nmicro
      IF (callmufi) THEN
         call dumptracers(micro_indx)
      ENDIF

!     Processing modern traceur options
      if(moderntracdef) then
        call ini_recombin
      endif

!------------------------------------------------------------
!     Initialisation tracers ....
!------------------------------------------------------------
      ! rho_q(1:nq)=0

      rho_n2=1000.   ! N2 ice density (kg.m-3)

      lw_co=274000.
      lw_ch4=586700.
      lw_n2=2.5e5
      write(*,*) "lw_n2 = ", lw_n2

      if (callmufi) then
        if (optichaze) then
          iaero_haze = 2
          write(*,*) 'Microphysical moment model'
          write(*,*) '--- number of haze aerosol = ', iaero_haze
        endif ! end optichaze
      
      elseif (haze) then
        ! the sedimentation radius remains radius(igcm_haze)
        if (fractal) then
          nmono=nb_monomer
        else
          nmono=1
        endif ! end fractal
        
        ia = 0
        if (optichaze) then
          ia = ia + 1
          iaero_haze = ia
          write(*,*) '--- number of haze aerosol = ', iaero_haze
          block = 0 ! Only one type of haze is active : the first one set in traceur.def
          do iq=1,nq
            tracername=noms(iq)
            write(*,*) "--> tracername ",iq,'/',nq,' = ',tracername
            if (tracername(1:4).eq."haze".and.block.eq.0) then
              i_haze=iq
              block=1
              write(*,*) "i_haze=",i_haze
              write(*,*) "Careful: if you set many haze traceurs in &
              traceur.def,only ",tracername," will be radiatively active &
              (first one in traceur.def)"
            endif
            enddo
        endif ! end optichaze
      endif ! end callmufi or haze

!     Output for records:
!     ~~~~~~~~~~~~~~~~~~
      write(*,*)
      Write(*,*) '******** initracer : dust transport parameters :'
      write(*,*) 'alpha_lift = ', alpha_lift
      write(*,*) 'alpha_devil = ', alpha_devil
      write(*,*) 'radius  = ', radius
      write(*,*) 'Qext  = ', qext
      write(*,*)

      contains

      subroutine get_tracdat(iq,tracline)
        !-------------------ADDING NEW OPTIONS-------------------
        ! Duplicate if sentence for adding new options
        ! Do not forget to add the new variables in tracer_h.F90
        ! Do not forget to allocate and initialize the new variables above
        ! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern"
        !-------------------------------------------------------
          use tracer_h
          implicit none
          integer,           intent(in) :: iq       ! tracer index
          character(len=500),intent(in) :: tracline ! traceur.def lines with parameters

          read(tracline,*) noms(iq)
          ! JVO 20 : We should add a sanity check aborting when duplicates in names !
          write(*,*)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
          ! option mmol
          if (index(tracline,'mmol='   ) /= 0) then
              read(tracline(index(tracline,'mmol=')+len('mmol='):),*)&
                  mmol(iq)
              write(*,*) ' Parameter value (traceur.def) : mmol=', &
                  mmol(iq)
          else
              write(*,*) ' Parameter value (default)     : mmol=', &
                 mmol(iq)
          end if
          ! option aki
          if (index(tracline,'aki='   ) /= 0) then
              read(tracline(index(tracline,'aki=')+len('aki='):),*) &
                  aki(iq)
              write(*,*) ' Parameter value (traceur.def) : aki=', &
                  aki(iq)
          else
              write(*,*) ' Parameter value (default)     : aki=', &
                  aki(iq)
          end if
          ! option cpi
          if (index(tracline,'cpi='   ) /= 0) then
              read(tracline(index(tracline,'cpi=')+len('cpi='):),*) &
                  cpi(iq)
              write(*,*) ' Parameter value (traceur.def) : cpi=', &
                  cpi(iq)
          else
              write(*,*) ' Parameter value (default)     : cpi=', &
                  cpi(iq)
          end if
          ! option is_chim
          if (index(tracline,'is_chim='   ) /= 0) then
          read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*) &
                  is_chim(iq)
              write(*,*) ' Parameter value (traceur.def) : is_chim=', &
                  is_chim(iq)
          else
              write(*,*) ' Parameter value (default)     : is_chim=', &
                  is_chim(iq)
          end if
          ! option is_rad
          if (index(tracline,'is_rad=') /= 0) then
              read(tracline(index(tracline,'is_rad=') &
              +len('is_rad='):),*) is_rad(iq)
              write(*,*) ' Parameter value (traceur.def) : is_rad=', &
              is_rad(iq)
          else
              write(*,*) ' Parameter value (default)     : is_rad=',  &
              is_rad(iq)
          end if
          ! option is_recomb
          if (index(tracline,'is_recomb=') /= 0) then
              read(tracline(index(tracline,'is_recomb=') &
              +len('is_recomb='):),*) is_recomb(iq)
              write(*,*) ' Parameter value (traceur.def) : is_recomb=', &
              is_recomb(iq)
          else
              write(*,*) ' Parameter value (default)     : is_recomb=', &
              is_recomb(iq)
          end if
          ! option is_recomb_qset
          if (index(tracline,'is_recomb_qset=') /= 0) then
              read(tracline(index(tracline,'is_recomb_qset=') &
              +len('is_recomb_qset='):),*) is_recomb_qset(iq)
              write(*,*) ' Parameter value (traceur.def) :'// &
              ' is_recomb_qset=', &
              is_recomb_qset(iq)
          else
              write(*,*) ' Parameter value (default)     :'// &
              ' is_recomb_qset=', &
              is_recomb_qset(iq)
          endif
          ! option is_recomb_qotf
          if (index(tracline,'is_recomb_qotf=') /= 0) then
              read(tracline(index(tracline,'is_recomb_qotf=') &
              +len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
              write(*,*) ' Parameter value (traceur.def) :'// &
              ' is_recomb_qotf=', &
              is_recomb_qotf(iq)
          else
              write(*,*) ' Parameter value (default)     :'// &
              ' is_recomb_qotf=',  &
              is_recomb_qotf(iq)
          end if
          !option radius
          if (index(tracline,'radius=') .ne. 0) then
            read(tracline(index(tracline,'radius=') &
              +len('radius='):),*) radius(iq)
            write(*,*)'Parameter value (traceur.def) :'// &
              "radius=",radius(iq), "m"
          else
            write(*,*) ' Parameter value (default)     :'// &
            ' radius=', radius(iq)," m"
          endif
          !option rho
          if (index(tracline,'rho=') .ne. 0) then
            read(tracline(index(tracline,'rho=') &
              +len('rho='):),*) rho_q(iq)
            write(*,*)'Parameter value (traceur.def) :'// &
              "rho=",rho_q(iq)
          else
            write(*,*) ' Parameter value (default)     :'// &
              ' rho=', rho_q(iq)
          endif
      end subroutine get_tracdat

      end