MODULE initracer_mod

implicit none

contains

SUBROUTINE initracer(ngrid,nq)

use tracer_mod
use comcstfi_h,     only: pi
use dust_param_mod, only: doubleq, submicron, dustbin
use callkeys_mod,   only: water, microphys, scavenging, rdstorm, &
                          topflows, photochem, callthermos, hdo, &
                          callnlte, nltemodel, dustinjection,    &
                          co2clouds, co2useh2o, meteo_flux

implicit none

! == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
!   subject:
!   --------
!   Initialization related to tracer
!   (transported dust, water, chemical species, ice...)
!
!   Name of the tracer
!
!   Test of dimension :
!   Initialize tracer related data in tracer_mod, using tracer names provided
!   by the dynamics in "infotrac"
!
!
!   author: F.Forget
!   ------
!    Modifs: Franck Montmessin, Sebastien Lebonnois (june 2003)
!            Ehouarn Millour (oct. 2008) identify tracers by their names
!            JB Clement (juillet 2025), upgrading to F90 with module
! == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =


integer, intent(in) :: ngrid ! number of atmospheric columns
integer, intent(in) :: nq    ! number of tracers

integer :: iq,ig,count
real    :: r0_lift, reff_lift, nueff_lift
real    :: r0_storm, reff_storm
! Ratio of small over large dust particles (used when both
!   doubleq and the submicron mode are active); In Montmessin
!   et al. (2002), a value of 25 has been deduced;
real, parameter                  :: popratio = 25.
character(len=30)                :: txt ! to store some text
character(len=30), dimension(nq) :: idtracers

!-----------------------------------------------------------------------
!  radius(nq)      ! aerosol particle radius (m)
!  rho_q(nq)       ! tracer densities (kg.m-3)
!  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
!-----------------------------------------------------------------------

! Sanity check: check that we are running with at least 1 tracer:
if (nq<1) then
    write(*,*) "initracer error, nq=",nq," must be >=1!"
    call abort_physic("initracer","nq<1",1)
endif
!------------------------------------------------------------
!         NAME and molar mass of the tracer
!------------------------------------------------------------

! Identify tracers by their names: (and set corresponding values of mmol)
! 0. initialize tracer indexes to zero:
igcm_dustbin(1:nq)=0
igcm_co2_ice=0
igcm_ccnco2_mass=0
igcm_ccnco2_number=0
igcm_ccnco2_meteor_mass=0
igcm_ccnco2_meteor_number=0
igcm_ccnco2_h2o_mass_ice=0
igcm_ccnco2_h2o_mass_ccn=0
igcm_ccnco2_h2o_number=0
igcm_dust_mass=0
igcm_dust_number=0
igcm_ccn_mass=0
igcm_ccn_number=0
igcm_dust_submicron=0
igcm_h2o_vap=0
igcm_h2o_ice=0
igcm_hdo_vap=0
igcm_hdo_ice=0
igcm_stormdust_mass=0
igcm_stormdust_number=0
igcm_topdust_mass=0
igcm_topdust_number=0
igcm_co2=0
igcm_co=0
igcm_o=0
igcm_o1d=0
igcm_o2=0
igcm_o3=0
igcm_h=0
igcm_d=0
igcm_hd=0
igcm_h2=0
igcm_od=0
igcm_do2=0
igcm_hdo2=0
igcm_oh=0
igcm_ho2=0
igcm_h2o2=0
igcm_ch4=0
igcm_n2=0
igcm_ar=0
igcm_ar_n2=0
igcm_n=0
igcm_no=0
igcm_no2=0
igcm_n2d=0
igcm_he=0
igcm_co2plus=0
igcm_oplus=0
igcm_o2plus=0
igcm_coplus=0
igcm_cplus=0
igcm_nplus=0
igcm_noplus=0
igcm_n2plus=0
igcm_hplus=0
igcm_hco2plus=0
igcm_hcoplus=0
igcm_h2oplus=0
igcm_h3oplus=0
igcm_ohplus=0
igcm_elec=0

! 1. find dust tracers
count=0
if (dustbin>0) then
    do iq=1,nq
        txt=" "
        write(txt,'(a4,i2.2)')'dust',count+1
        if (noms(iq) == txt) then
        count = count + 1
        idtracers(count) = txt
        igcm_dustbin(count)=iq
        mmol(iq)=100.
        endif
    enddo !do iq=1,nq
endif ! of if (dustbin>0)
if (doubleq) then
    do iq=1,nq
        if (noms(iq) == "dust_mass") then
            igcm_dust_mass=iq
            count = count + 1
            idtracers(count) = "dust_mass"
        endif
        if (noms(iq) == "dust_number") then
            igcm_dust_number=iq
            count = count + 1
            idtracers(count) = "dust_number"
        endif
    enddo
endif ! of if (doubleq)
if (microphys) then
    do iq=1,nq
        if (noms(iq) == "ccn_mass") then
            igcm_ccn_mass=iq
            count = count + 1
            idtracers(count) = "ccn_mass"
        endif
        if (noms(iq) == "ccn_number") then
            igcm_ccn_number=iq
            count = count + 1
            idtracers(count) = "ccn_number"
        endif
    enddo
endif ! of if (microphys)
if (submicron) then
    do iq=1,nq
        if (noms(iq) == "dust_submicron") then
            igcm_dust_submicron=iq
            mmol(iq)=100.
            count = count + 1
            idtracers(count) = "dust_submicron"
        endif
    enddo
endif ! of if (submicron)
if (rdstorm) then
    do iq=1,nq
        if (noms(iq) == "stormdust_mass") then
            igcm_stormdust_mass=iq
            count = count + 1
            idtracers(count) = "stormdust_mass"
        endif
        if (noms(iq) == "stormdust_number") then
            igcm_stormdust_number=iq
            count = count + 1
            idtracers(count) = "stormdust_number"
        endif
    enddo
endif ! of if (rdstorm)
if (topflows) then
    do iq=1,nq
        if (noms(iq) == "topdust_mass") then
            igcm_topdust_mass=iq
            count = count + 1
            idtracers(count) = "topdust_mass"
        endif
        if (noms(iq) == "topdust_number") then
            igcm_topdust_number=iq
            count = count + 1
            idtracers(count) = "topdust_number"
        endif
    enddo
endif ! of if (topflows)
! 2. find chemistry and water tracers
do iq=1,nq
    if (noms(iq) == "co2") then
        igcm_co2=iq
        mmol(igcm_co2)=44.
        count = count + 1
        idtracers(count) = "co2"
    endif
    if (noms(iq) == "co") then
        igcm_co=iq
        mmol(igcm_co)=28.
        count = count + 1
        idtracers(count) = "co"
    endif
    if (noms(iq) == "o") then
        igcm_o=iq
        mmol(igcm_o)=16.
        count = count + 1
        idtracers(count) = "o"
    endif
    if (noms(iq) == "o1d") then
        igcm_o1d=iq
        mmol(igcm_o1d)=16.
        count = count + 1
        idtracers(count) = "o1d"
    endif
    if (noms(iq) == "o2") then
        igcm_o2=iq
        mmol(igcm_o2)=32.
        count = count + 1
        idtracers(count) = "o2"
    endif
    if (noms(iq) == "o3") then
        igcm_o3=iq
        mmol(igcm_o3)=48.
        count = count + 1
        idtracers(count) = "o3"
    endif
    if (noms(iq) == "h") then
        igcm_h=iq
        mmol(igcm_h)=1.
        count = count + 1
        idtracers(count) = "h"
    endif
    if (noms(iq) == "h2") then
        igcm_h2=iq
        mmol(igcm_h2)=2.
        count = count + 1
        idtracers(count) = "h2"
    endif
    if (noms(iq) == "oh") then
        igcm_oh=iq
        mmol(igcm_oh)=17.
        count = count + 1
        idtracers(count) = "oh"
    endif
    if (noms(iq) == "ho2") then
        igcm_ho2=iq
        mmol(igcm_ho2)=33.
        count = count + 1
        idtracers(count) = "ho2"
    endif
    if (noms(iq) == "h2o2") then
        igcm_h2o2=iq
        mmol(igcm_h2o2)=34.
        count = count + 1
        idtracers(count) = "h2o2"
    endif
    if (noms(iq) == "n2") then
        igcm_n2=iq
        mmol(igcm_n2)=28.
        count = count + 1
        idtracers(count) = "n2"
    endif
    if (noms(iq) == "ch4") then
        igcm_ch4=iq
        mmol(igcm_ch4)=16.
        count = count + 1
        idtracers(count) = "ch4"
    endif
    if (noms(iq) == "ar") then
        igcm_ar=iq
        mmol(igcm_ar)=40.
        count = count + 1
        idtracers(count) = "ar"
    endif
    if (noms(iq) == "n") then
        igcm_n=iq
        mmol(igcm_n)=14.
        count = count + 1
        idtracers(count) = "n"
    endif
    if (noms(iq) == "no") then
        igcm_no=iq
        mmol(igcm_no)=30.
        count = count + 1
        idtracers(count) = "no"
    endif
    if (noms(iq) == "no2") then
        igcm_no2=iq
        mmol(igcm_no2)=46.
        count = count + 1
        idtracers(count) = "no2"
    endif
    if (noms(iq) == "n2d") then
        igcm_n2d=iq
        mmol(igcm_n2d)=28.
        count = count + 1
        idtracers(count) = "n2d"
    endif
    if (noms(iq) == "he") then
        igcm_he=iq
        mmol(igcm_he)=4.
        count = count + 1
        idtracers(count) = "he"
    endif
    if (noms(iq) == "co2plus") then
        igcm_co2plus=iq
        mmol(igcm_co2plus)=44.
        count = count + 1
        idtracers(count) = "co2plus"
    endif
    if (noms(iq) == "oplus") then
        igcm_oplus=iq
        mmol(igcm_oplus)=16.
        count = count + 1
        idtracers(count) = "oplus"
    endif
    if (noms(iq) == "o2plus") then
        igcm_o2plus=iq
        mmol(igcm_o2plus)=32.
        count = count + 1
        idtracers(count) = "o2plus"
    endif
    if (noms(iq) == "coplus") then
        igcm_coplus=iq
        mmol(igcm_coplus)=28.
        count = count + 1
        idtracers(count) = "coplus"
    endif
    if (noms(iq) == "cplus") then
        igcm_cplus=iq
        mmol(igcm_cplus)=12.
        count = count + 1
        idtracers(count) = "cplus"
    endif
    if (noms(iq) == "nplus") then
        igcm_nplus=iq
        mmol(igcm_nplus)=14.
        count = count + 1
        idtracers(count) = "nplus"
    endif
    if (noms(iq) == "noplus") then
        igcm_noplus=iq
        mmol(igcm_noplus)=30.
        count = count + 1
        idtracers(count) = "noplus"
    endif
    if (noms(iq) == "n2plus") then
        igcm_n2plus=iq
        mmol(igcm_n2plus)=28.
        count = count + 1
        idtracers(count) = "n2plus"
    endif
    if (noms(iq) == "hplus") then
        igcm_hplus=iq
        mmol(igcm_hplus)=1.
        count = count + 1
        idtracers(count) = "hplus"
    endif
    if (noms(iq) == "hco2plus") then
        igcm_hco2plus=iq
        mmol(igcm_hco2plus)=45.
        count = count + 1
        idtracers(count) = "hco2plus"
    endif
    if (noms(iq) == "hcoplus") then
        igcm_hcoplus=iq
        mmol(igcm_hcoplus)=29.
        count = count + 1
        idtracers(count) = "hcoplus"
    endif
    if (noms(iq) == "h2oplus") then
        igcm_h2oplus=iq
        mmol(igcm_h2oplus)=18.
        count = count + 1
        idtracers(count) = "h2oplus"
    endif
    if (noms(iq) == "h3oplus") then
        igcm_h3oplus=iq
        mmol(igcm_h3oplus)=19.
        count = count + 1
        idtracers(count) = "h3oplus"
    endif
    if (noms(iq) == "ohplus") then
        igcm_ohplus=iq
        mmol(igcm_ohplus)=17.
        count = count + 1
        idtracers(count) = "ohplus"
    endif
    if (noms(iq) == "elec") then
        igcm_elec=iq
        mmol(igcm_elec)=1./1822.89
        count = count + 1
        idtracers(count) = "elec"
    endif
    if (noms(iq) == "h2o_vap") then
        igcm_h2o_vap=iq
        mmol(igcm_h2o_vap)=18.
        count = count + 1
        idtracers(count) = "h2o_vap"
    endif
    if (noms(iq) == "hdo_vap") then
        igcm_hdo_vap=iq
        mmol(igcm_hdo_vap)=19.
        count = count + 1
        idtracers(count) = "hdo_vap"
    endif
    if (noms(iq) == "od") then
        igcm_od=iq
        mmol(igcm_od)=18.
        count = count + 1
        idtracers(count) = "od"
    endif
    if (noms(iq) == "d") then
        igcm_d=iq
        mmol(igcm_d)=2.
        count = count + 1
        idtracers(count) = "d"
    endif
    if (noms(iq) == "hd") then
        igcm_hd=iq
        mmol(igcm_hd)=3.
        count = count + 1
        idtracers(count) = "hd"
    endif
    if (noms(iq) == "do2") then
        igcm_do2=iq
        mmol(igcm_do2)=34.
        count = count + 1
        idtracers(count) = "do2"
    endif
    if (noms(iq) == "hdo2") then
        igcm_hdo2=iq
        mmol(igcm_hdo2)=35.
        count = count + 1
        idtracers(count) = "hdo2"
    endif
    if (noms(iq) == "co2_ice") then
        igcm_co2_ice=iq
        mmol(igcm_co2_ice)=44.
        count = count + 1
        idtracers(count) = "co2_ice"
    endif
    if (noms(iq) == "h2o_ice") then
        igcm_h2o_ice=iq
        mmol(igcm_h2o_ice)=18.
        count = count + 1
        idtracers(count) = "h2o_ice"
    endif
    if (noms(iq) == "hdo_ice") then
        igcm_hdo_ice=iq
        mmol(igcm_hdo_ice)=19.
        count = count + 1
        idtracers(count) = "hdo_ice"
    endif
    ! Other stuff: e.g. for simulations using co2 + neutral gaz
    if (noms(iq) == "Ar_N2") then
        igcm_ar_n2=iq
        mmol(igcm_ar_n2)=30.
        count = count + 1
        idtracers(count) = "Ar_N2"
    endif
    if (co2clouds) then
        if (noms(iq) == "ccnco2_mass") then
            igcm_ccnco2_mass=iq
            count = count + 1
            idtracers(count) = "ccnco2_mass"
        endif
        if (noms(iq) == "ccnco2_number") then
            igcm_ccnco2_number=iq
            count = count + 1
            idtracers(count) = "ccnco2_number"
        endif
        if (meteo_flux) then
            if (noms(iq) == "ccnco2_meteor_mass") then
                igcm_ccnco2_meteor_mass=iq
                count = count + 1
                idtracers(count) = "ccnco2_meteor_mass"
            endif
            if (noms(iq) == "ccnco2_meteor_number") then
                igcm_ccnco2_meteor_number=iq
                count = count + 1
                idtracers(count) = "ccnco2_meteor_number"
            endif
        end if
        if (co2useh2o) then
            if (noms(iq) == "ccnco2_h2o_number") then
                igcm_ccnco2_h2o_number=iq
                count = count + 1
                idtracers(count) = "ccnco2_h2o_number"
            endif
            if (noms(iq) == "ccnco2_h2o_mass_ice") then
                igcm_ccnco2_h2o_mass_ice=iq
                count = count + 1
                idtracers(count) = "ccnco2_h2o_mass_ice"
            endif
            if (noms(iq) == "ccnco2_h2o_mass_ccn") then
                igcm_ccnco2_h2o_mass_ccn=iq
                count = count + 1
                idtracers(count) = "ccnco2_h2o_mass_ccn"
            endif
        end if
    endif
enddo ! of do iq=1,nq

! check that we identified all tracers:
if (count /= nq) then
    write(*,*) "initracer: found only",count,"tracers but",nq,"are expected according to traceur.def!"
    write(*,*) "initracer: the only identified tracers are:"
    do iq = 1,count
        write(*,*) '      ',iq,trim(idtracers(iq))
    enddo
    call abort_physic("initracer","tracer mismatch",1)
else
    write(*,*) "initracer: found all expected tracers, namely:"
    do iq = 1,nq
        write(*,*) '      ',iq,' ',trim(noms(iq))
    enddo
endif

! Additional test required for HDO
! We need to compute some things for H2O before HDO
if (hdo) then
    if (igcm_h2o_vap>igcm_hdo_vap) then
        call abort_physic("initracer","Tracer H2O must be initialized before HDO",1)
    else if ((nqfils(igcm_h2o_ice)<1) .or. (nqfils(igcm_h2o_vap)<1)) then
        write(*,*) "HDO must be transported as a son of H2O: complete traceur.def"
        call abort_physic("initracer","adapt your tracer.def",1)
    else if ((igcm_hdo_ice < nq-2) .or. (igcm_hdo_vap < nq-2)) then
        write(*,*) "The isotopes (HDO) must be placed at the end of the list in traceur.def"
        call abort_physic("initracer","adapt your tracer.def",1)
    endif
endif

!------------------------------------------------------------
!     Initialize tracers .... (in tracer_mod)
!------------------------------------------------------------
! start by setting everything to (default) zero
rho_q(1:nq)=0     ! tracer density (kg.m-3)
radius(1:nq)=0.   ! tracer particle radius (m)
alpha_lift(1:nq) =0.  ! tracer saltation vertical flux/horiz flux ratio (m-1)
alpha_devil(1:nq)=0.  ! tracer lifting coefficient by dust devils


! some reference values
rho_dust=2500.  ! Mars dust density (kg.m-3)
rho_ice=920.    ! Water ice density (kg.m-3)
rho_ice_co2=1650.

if (doubleq) then
    ! verify that we indeed have dust_mass and dust_number tracers
    if (igcm_dust_mass == 0) then
           write(*,*) "initracer: error!"
            write(*,*) "  Cannot use doubleq option without a dust_mass tracer!"
           call abort_physic("initracer","doubleq issue",1)
    endif
    if (igcm_dust_number == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use doubleq option without a dust_number tracer!"
        call abort_physic("initracer","doubleq issue",1)
    endif

    ! "doubleq" technique
    ! -------------------
    ! transport of mass and number mixing ratio)
    ! iq=1: Q mass mixing ratio, iq=2: N number mixing ratio

    if( (nq<2).or.(water.and.(nq<4)) .or.(hdo.and.(nq<6) )) then
        write(*,*)'initracer: nq is too low : nq=', nq
        write(*,*)'water= ',water,' doubleq= ',doubleq
    end if

    nueff_lift = 0.5
    varian=sqrt(log(1.+nueff_lift))

    rho_q(igcm_dust_mass)=rho_dust
    rho_q(igcm_dust_number)=rho_dust

    ! Intermediate calcul for computing geometric mean radius r0
    ! as a function of mass and number mixing ratio Q and N
    ! (r0 = (r3n_q * Q/ N)^(1/3))
    r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust)

    ! Intermediate calcul for computing effective radius reff
    ! from geometric mean radius r0
    ! (reff = ref_r0 * r0)
    ref_r0 = exp(2.5*varian**2)

    ! lifted dust :
    ! '''''''''''
    reff_lift = 3.0e-6 !3.e-6 !Effective radius of lifted dust (m)
    alpha_devil(igcm_dust_mass)=9.e-9   !  dust devil lift mass coeff
    ! alpha_lift(igcm_dust_mass)=3.0e-15  !  Lifted mass coeff

    ! default lifting settings
    ! -- GCM: alpha_lift not zero because large-scale lifting by default
    ! -- MESOSCALE: alpha_lift zero because no lifting at all in mesoscale by default
#ifdef MESOSCALE
    alpha_lift(igcm_dust_mass)=0.0
#else
    alpha_lift(igcm_dust_mass)=1.e-6 !1.e-6 !Lifted mass coeff
    if (dustinjection >= 1) then
        reff_lift = 3.0e-6 ! Effective radius of lifted dust (m)
        alpha_lift(igcm_dust_mass)=(4/3.)*reff_lift*rho_dust/2.4
    endif
#endif

    r0_lift = reff_lift/ref_r0
    alpha_devil(igcm_dust_number)=r3n_q*alpha_devil(igcm_dust_mass)/r0_lift**3
    alpha_lift(igcm_dust_number)=r3n_q*alpha_lift(igcm_dust_mass)/r0_lift**3

    radius(igcm_dust_mass) = reff_lift
    radius(igcm_dust_number) = reff_lift

    write(*,*) "initracer: doubleq_param reff_lift:", reff_lift
    write(*,*) "initracer: doubleq_param nueff_lift:", nueff_lift
    write(*,*) "initracer: doubleq_param alpha_lift:",alpha_lift(igcm_dust_mass)
    ! ----------------------------------------------------------------------
    ! rocket dust storm scheme
    ! lifting tracer stormdust using same distribution than
    ! normal dust
    if (rdstorm) then
        reff_storm=3.e-6 ! reff_lift !3.e-6
        r0_storm=reff_storm/ref_r0
        rho_q(igcm_stormdust_mass)=rho_dust
        rho_q(igcm_stormdust_number)=rho_dust

        alpha_devil(igcm_stormdust_mass)=9.e-9   ! dust devil lift mass coeff
        alpha_lift(igcm_stormdust_mass)=4./3./2.4*reff_storm*rho_dust

        write(*,*) 'alpha_lift(rds):',alpha_lift(igcm_stormdust_mass)

        alpha_devil(igcm_stormdust_number)=r3n_q*alpha_devil(igcm_stormdust_mass)/r0_storm**3
        alpha_lift(igcm_stormdust_number)=r3n_q*alpha_lift(igcm_stormdust_mass)/r0_storm**3

        radius(igcm_stormdust_mass) = reff_storm
        radius(igcm_stormdust_number) = reff_storm
    end if !(rdstorm)
    ! ----------------------------------------------------------------------
    ! mountain top dust flows scheme
    ! you need a radius value for topdust to active its sedimentation
    ! we take the same value as for the normal dust
    if (topflows) then
        rho_q(igcm_topdust_mass)=rho_dust
        rho_q(igcm_topdust_number)=rho_dust
        radius(igcm_topdust_mass) = 3.e-6
        radius(igcm_topdust_number) = 3.e-6
    end if !(topflows)
    ! ----------------------------------------------------------------------

else

    ! initialize varian, which may be used (e.g. by surfacearea)
    ! even with conrath dust
    nueff_lift = 0.5
    varian=sqrt(log(1.+nueff_lift))

    if (dustbin>1) then
        print*,'initracer: STOP! properties of dust need to be set in initracer!!!'
        call abort_physic("initracer","dustbin properties issue",1)

    else if (dustbin == 1) then

!       This will be used for 1 dust particle size:
!       ------------------------------------------
        radius(igcm_dustbin(1))=3.e-6
        alpha_lift(igcm_dustbin(1))=0.0e-6
        alpha_devil(igcm_dustbin(1))=7.65e-9
        rho_q(igcm_dustbin(1))=rho_dust

    endif
end if    ! (doubleq)


! Scavenging of dust particles by H2O clouds:
! ------------------------------------------
! Initialize the two tracers used for the CCNs
if (water.and.doubleq.and.scavenging) then
    radius(igcm_ccn_mass) = radius(igcm_dust_mass)
    alpha_lift(igcm_ccn_mass) = 1e-30
    alpha_devil(igcm_ccn_mass) = 1e-30
    rho_q(igcm_ccn_mass) = rho_dust

    radius(igcm_ccn_number) = radius(igcm_ccn_mass)
    alpha_lift(igcm_ccn_number) = alpha_lift(igcm_ccn_mass)
    alpha_devil(igcm_ccn_number) = alpha_devil(igcm_ccn_mass)
    rho_q(igcm_ccn_number) = rho_q(igcm_ccn_mass)
endif ! of if (water.AND.doubleq.AND.scavenging)

! Submicron dust mode:
! --------------------
if (submicron) then
    radius(igcm_dust_submicron)=0.1e-6
    rho_q(igcm_dust_submicron)=rho_dust
    if (doubleq) then
!        If doubleq is also active, we use the population ratio:
        alpha_lift(igcm_dust_submicron) = alpha_lift(igcm_dust_number)*popratio*rho_q(igcm_dust_submicron)*4./3.*pi*radius(igcm_dust_submicron)**3.
        alpha_devil(igcm_dust_submicron)=1.e-30
    else
        alpha_lift(igcm_dust_submicron)=1e-6
        alpha_devil(igcm_dust_submicron)=1.e-30
    endif ! (doubleq)
end if  ! (submicron)

! Initialization for water vapor
! ------------------------------
if (water) then
    radius(igcm_h2o_vap)=0.
    alpha_lift(igcm_h2o_vap) =0.
    alpha_devil(igcm_h2o_vap)=0.
    if (water.and.(nq>=2)) then
        radius(igcm_h2o_ice)=3.e-6
        rho_q(igcm_h2o_ice)=rho_ice
        alpha_lift(igcm_h2o_ice) =0.
        alpha_devil(igcm_h2o_ice)=0.
    else if(water.and.(nq<2)) then
        write(*,*) 'nq is too low : nq=', nq
        write(*,*) 'water= ',water
    endif
end if  ! (water)

! Initialization for hdo vapor
! ------------------------------
if (hdo) then
    radius(igcm_hdo_vap)=0.
    alpha_lift(igcm_hdo_vap) =0.
    alpha_devil(igcm_hdo_vap)=0.
    if (water.and.(nq>=2)) then
        radius(igcm_hdo_ice)=3.e-6
        rho_q(igcm_hdo_ice)=rho_ice
        alpha_lift(igcm_hdo_ice) =0.
        alpha_devil(igcm_hdo_ice)=0.
    else if (hdo.and.(nq<6)) then
        write(*,*) 'nq is too low : nq=', nq
        write(*,*) 'hdo= ',hdo
    endif
end if  ! (hdo)

! Output for records:
! ~~~~~~~~~~~~~~~~~~
write(*,*)
Write(*,*) '******** initracer : dust transport parameters :'
write(*,*) 'alpha_lift = ', alpha_lift
write(*,*) 'alpha_devil = ', alpha_devil
write(*,*) 'radius  = ', radius
if(doubleq) then
    write(*,*) 'reff_lift (um) =  ', reff_lift
    write(*,*) 'size distribution variance  = ', varian
    write(*,*) 'r3n_q , ref_r0 : ', r3n_q , ref_r0
end if

!
! some extra (possibly redundant) sanity checks for tracers:
! ---------------------------------------------------------
if ((.not.doubleq).and.(dustbin>0)) then
    ! verify that we indeed have 'dustbin' dust tracers
    count=0
    do iq=1,dustbin
        if (igcm_dustbin(iq)/=0) then
            count = count + 1
        endif
    enddo
    if (count/=dustbin) then
        write(*,*) "initracer: error!"
        write(*,*) "  dustbin is set to ",dustbin," but we only have the following dust tracers:"
        do iq=1,count
            write(*,*)"   ",trim(noms(igcm_dustbin(iq)))
        enddo
        call abort_physic("initracer","dustbin issue",1)
    endif
endif

if (water) then
    ! verify that we indeed have h2o_vap and h2o_ice tracers
    if (igcm_h2o_vap == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use water option without an h2o_vap tracer!"
        call abort_physic("initracer","water cycle issue",1)
    endif
    if (igcm_h2o_ice == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use water option without an h2o_ice tracer!"
        call abort_physic("initracer","water cycle issue",1)
    endif
endif

if (hdo) then
    ! verify that we indeed have hdo_vap and hdo_ice tracers
    if (igcm_hdo_vap == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use hdo option without an hdo_vap tracer!"
        call abort_physic("initracer","hdo cycle issue",1)
    endif
    if (igcm_hdo_ice == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use hdo option without an hdo_ice tracer!"
        call abort_physic("initracer","hdo cycle issue",1)
    endif
endif

if (microphys) then
    ! verify that we indeed have ccn_number and ccn_mass tracers
    if (igcm_ccn_number == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use microphys option without a ccn_number tracer!"
        call abort_physic("initracer","microphys issue",1)
    endif
    if (igcm_ccn_mass == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use microphys option without a ccn_mass tracer!"
        call abort_physic("initracer","microphys issue",1)
    endif
endif

if (co2clouds) then
    !verify that we have co2_ice and co2 tracers
    if (igcm_co2 == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use co2 clouds option without a co2 tracer!"
        call abort_physic("initracer","co2 clouds issue",1)
    end if
    if (igcm_co2_ice == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use co2 clouds option without a co2_ice tracer!"
        call abort_physic("initracer","co2 clouds issue",1)
    end if
    if (igcm_ccnco2_number == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use co2 clouds option without a ccnco2_number tracer!"
        call abort_physic("initracer","co2 clouds issue",1)
    end if
    if (igcm_ccnco2_mass == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use co2 clouds option without a ccnco2_mass tracer!"
        call abort_physic("initracer","co2 clouds issue",1)
    end if
    if (co2useh2o) then
        if (igcm_ccnco2_h2o_number == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  Cannot use co2 clouds option without a ccnco2_h2o_number tracer!"
            call abort_physic("initracer","co2 clouds issue",1)
        end if
        if (igcm_ccnco2_h2o_mass_ice == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  Cannot use co2 clouds option without a ccnco2_h2o_mass_ice tracer!"
            call abort_physic("initracer","co2 clouds issue",1)
        end if
        if (igcm_ccnco2_h2o_mass_ccn == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  Cannot use co2 clouds option without a ccnco2_h2o_mass_ccn tracer!"
            call abort_physic("initracer","co2 clouds issue",1)
        end if
    end if
    if (meteo_flux) then
        if (igcm_ccnco2_meteor_number == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  Cannot use co2 clouds option without a ccnco2_meteor_number tracer!"
            call abort_physic("initracer","co2 clouds issue",1)
        end if
        if (igcm_ccnco2_meteor_mass == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  Cannot use co2 clouds option without a ccnco2_h2o_mass_ice tracer!"
            call abort_physic("initracer","co2 clouds issue",1)
        end if
        if (igcm_ccnco2_meteor_mass == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  Cannot use co2 clouds option without a ccnco2_meteor_mass tracer!"
            call abort_physic("initracer","co2 clouds issue",1)
        end if
    end if
endif

if (rdstorm) then
    ! verify that we indeed have stormdust_mass and stormdust_number tracers
    if (igcm_stormdust_mass == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use rdstorm option without a stormdust_mass tracer!"
        call abort_physic("initracer","rdstorm issue",1)
    endif
    if (igcm_stormdust_number == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use rdstorm option without a stormdust_number tracer!"
        call abort_physic("initracer","rdstorm issue",1)
    endif
endif

if (topflows) then
    ! verify that we indeed have topdust_mass and topdust_number tracers
    if (igcm_topdust_mass == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use topflows option without a topdust_mass tracer!"
        call abort_physic("initracer","topflows issue",1)
    endif
    if (igcm_topdust_number == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use topflows option without a topdust_number tracer!"
        call abort_physic("initracer","topflows issue",1)
    endif
endif

if (callnlte) then ! NLTE requirements
    if (nltemodel>=1) then
        ! check that co2, co, o and n2 tracers are available
        if (igcm_co2 == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  with nltemodel>0, we need the co2 tracer!"
            call abort_physic("initracer","missing co2 tracer",1)
        endif
        if (igcm_co == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  with nltemodel>0, we need the co tracer!"
            call abort_physic("initracer","missing co tracer",1)
        endif
        if (igcm_o == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  with nltemodel>0, we need the o tracer!"
            call abort_physic("initracer","missing o tracer",1)
        endif
        if (igcm_n2 == 0) then
            write(*,*) "initracer: error!"
            write(*,*) "  with nltemodel>0, we need the n2 tracer!"
            call abort_physic("initracer","missing n2 tracer",1)
        endif
    endif
endif

if (scavenging) then
    ! verify that we indeed have ccn_mass and ccn_number tracers
    if (igcm_ccn_mass == 0 .and. igcm_ccnco2_mass == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use scavenging option without a ccn_mass or ccnco2_mass tracer!"
        call abort_physic("initracer","scavenging issue",1)
    endif
    if (igcm_ccn_number == 0 .and. igcm_ccnco2_number == 0 ) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use scavenging option without a ccn_number or ccnco2_number tracer!"
        call abort_physic("initracer","scavenging issue",1)
    endif
endif ! of if (scavenging)

if (photochem .or. callthermos) then
    ! verify that we indeed have the chemistry tracers
    if (igcm_co2 == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a co2 tracer!"
        call abort_physic("initracer","missing co2 tracer",1)
    endif
    if (igcm_co == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a co tracer!"
        call abort_physic("initracer","missing co tracer",1)
    endif
    if (igcm_o == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a o tracer!"
        call abort_physic("initracer","missing o tracer",1)
    endif
    if (igcm_o1d == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a o1d tracer!"
        call abort_physic("initracer","missing o1d tracer",1)
    endif
    if (igcm_o2 == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without an o2 tracer!"
        call abort_physic("initracer","missing o2 tracer",1)
    endif
    if (igcm_o3 == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without an o3 tracer!"
        call abort_physic("initracer","missing o3 tracer",1)
    endif
    if (igcm_h == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a h tracer!"
        call abort_physic("initracer","missing h tracer",1)
    endif
    if (igcm_h2 == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a h2 tracer!"
        call abort_physic("initracer","missing h2 tracer",1)
    endif
    if (igcm_oh == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without an oh tracer!"
        call abort_physic("initracer","missing oh tracer",1)
    endif
    if (igcm_ho2 == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a ho2 tracer!"
        call abort_physic("initracer","missing ho2 tracer",1)
    endif
    if (igcm_h2o2 == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a h2o2 tracer!"
        call abort_physic("initracer","missing h2o2 tracer",1)
    endif
    if (igcm_n2 == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without a n2 tracer!"
        call abort_physic("initracer","missing n2 tracer",1)
    endif
    if (igcm_ar == 0) then
        write(*,*) "initracer: error!"
        write(*,*) "  Cannot use chemistry option without an ar tracer!"
        call abort_physic("initracer","missing ar tracer",1)
    endif
endif ! of if (photochem .or. callthermos)

! Initialisation for CO2 clouds
if (co2clouds) then
    radius(igcm_ccnco2_mass) = radius(igcm_dust_mass)
    alpha_lift(igcm_ccnco2_mass) = 1e-30
    alpha_devil(igcm_ccnco2_mass) = 1e-30
    rho_q(igcm_ccnco2_mass) = rho_dust
    radius(igcm_ccnco2_number) = radius(igcm_ccnco2_mass)
    alpha_lift(igcm_ccnco2_number) = alpha_lift(igcm_ccnco2_mass)
    alpha_devil(igcm_ccnco2_number) = alpha_devil(igcm_ccnco2_mass)
    rho_q(igcm_ccnco2_number) = rho_q(igcm_ccnco2_mass)

    radius(igcm_co2)=0.
    alpha_lift(igcm_co2) =0.
    alpha_devil(igcm_co2)=0.
    radius(igcm_co2_ice)=1.e-8
    rho_q(igcm_co2_ice)=rho_ice_co2
    alpha_lift(igcm_co2_ice) =0.
    alpha_devil(igcm_co2_ice)=0.
endif

END SUBROUTINE initracer

END MODULE initracer_mod