source: trunk/LMDZ.GENERIC/libf/phystd/radioactive_tracers.F90 @ 3974

subroutine radioactive_tracers(ngrid,nlayer,nq,ptimestep,pq,zdqradio)

    USE tracer_h, only: half_life, top_prod, bot_prod, noms
    implicit none
    
    !==================================================================
    !     
    !     Purpose
    !     -------
    !     Calculates the decay of radioactive tracers
    !     
    !     Authors
    !     -------
    !     Maxime Maurice (14/04/2025)
    !     
    !==================================================================
    
    !     Arguments
    integer,intent(in) :: ngrid                      ! number of atmospheric columns
    integer,intent(in) :: nlayer                     ! number of atmospheric layers
    integer,intent(in) :: nq                         ! number of tracers
    real,intent(in)    :: ptimestep                  ! time interval
    real,intent(in)    :: pq(ngrid,nlayer,nq)        ! tracers (kg/kg)
    real,intent(out)   :: zdqradio(ngrid,nlayer,nq)  ! radioactive tracers tendencies

    !     Local variables
    integer iq

        zdqradio(1:ngrid,1:nlayer,1:nq) = 0
        do iq=1,nq
        if (half_life(iq) .ne. 0) then
            ! SOURCE
            ! bottom
            if (bot_prod(iq) .ne. 0) then
                zdqradio(1:ngrid,1,iq) = zdqradio(1:ngrid,1,iq) + bot_prod(iq) ! <== bottom production
            end if ! bot_prod != 0
            ! top
            if (top_prod(iq) .ne. 0) then
                zdqradio(1:ngrid,nlayer,iq) = zdqradio(1:ngrid,nlayer,iq) + top_prod(iq) ! <== top production (prod = log(2)/half_life)
            end if ! bot_prod != 0

            ! DECAY
            if (half_life(iq)/10 .le. ptimestep) then ! sanity test
                write(*,*)"timestep is more than half-life / 10 of tracer ", noms(iq)
                call abort_physic("radioactive_tracers","timestep is more than half-life / 10 of tracer",1)
            end if ! dt > 0.1 * half_life
            zdqradio(1:ngrid,1:nlayer,iq) = zdqradio(1:ngrid,1:nlayer,iq) - LOG(2.)/half_life(iq)*pq(1:ngrid,1:nlayer,iq) ! <== decay
        end if ! half_life(iq) != 0
        end do ! iq=1,nq

end subroutine radioactive_tracers