source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE/src_stomate/lpj_fire.f90 @ 2230

! calculate fire extent and impact on plants.
! This is treated on a pseudo-daily basis (fireindex has a long-term memory).
! We only take into account the natural litter.
! Grasses are totally burned.
! When the vegetation is dynamic, it also decreases the density of individuals.
! Fire burns litter on the ground.
!
! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/lpj_fire.f90,v 1.12 2010/04/06 15:44:01 ssipsl Exp $
! IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
MODULE lpj_fire

  ! modules used:

  USE ioipsl
  USE stomate_constants
  USE constantes_veg
  USE parallel

  IMPLICIT NONE

  ! private & public routines

  PRIVATE
  PUBLIC fire,fire_clear

  ! first call
  LOGICAL, SAVE                                                   :: firstcall = .TRUE.
  ! flag that disable fire
  LOGICAL, SAVE                                                   :: disable_fire

CONTAINS


  SUBROUTINE fire_clear
    firstcall = .TRUE.
  END SUBROUTINE fire_clear

  SUBROUTINE fire (npts, dt, litterpart, &
       litterhum_daily, t2m_daily, lignin_struc, &
       fireindex, firelitter, biomass, ind, &
       litter, dead_leaves, bm_to_litter, black_carbon, &
       co2_fire)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                             :: npts
    ! Time step in days
    REAL(r_std), INTENT(in)                                          :: dt
    ! fraction of litter above the ground belonging to different PFTs
    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in)              :: litterpart
    ! Daily litter moisture (between 0 and 1)
    REAL(r_std), DIMENSION(npts), INTENT(in)                         :: litterhum_daily
    ! Daily 2 meter temperature (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                         :: t2m_daily
    ! ratio Lignine/Carbon in structural litter, above and below ground,
    ! (gC/(m**2 of ground))
    REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(in)         :: lignin_struc

    ! 0.2 modified fields

    ! Probability of fire
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: fireindex
    ! Longer term total natural litter above the ground, gC/m**2 of natural ground
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: firelitter
    ! biomass (gC/m**2)
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)          :: biomass
    ! density of individuals (1/m**2)
    REAl(r_std), DIMENSION(npts,nvm), INTENT(inout)                 :: ind
    ! metabolic and structural litter, above and below ground (gC/m**2)
    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs), INTENT(inout):: litter
    ! dead leaves on ground, per PFT, metabolic and structural,
    !   in gC/(m**2 of ground)
    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)           :: dead_leaves
    ! conversion of biomass to litter (gC/(m**2 of average ground)) / day
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)          :: bm_to_litter
    ! black carbon on the ground (gC/(m**2 of total ground))
    REAL(r_std), DIMENSION(npts), INTENT(inout)                      :: black_carbon

    ! 0.3 output

    ! carbon emitted into the atmosphere by fire (living and dead biomass)
    ! (in gC/m**2/day)
    !NV devient 2D
    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                        :: co2_fire

    ! 0.4 local

    ! Time scale for memory of the fire index (days). Validated for one year in the DGVM.
    !MM Shilong ??
!!$    REAL(r_std), PARAMETER                                           :: tau_fire = 365.  ! GKtest
    REAL(r_std), PARAMETER                                           :: tau_fire = 30.  ! GKtest
    ! Critical litter quantity for fire
    REAL(r_std), PARAMETER                                           :: litter_crit = 200.
    ! fire perturbation
    REAL(r_std), DIMENSION(npts)                                     :: fire_disturb
    ! what fraction of the plants is burned each day?
    REAL(r_std), DIMENSION(npts,nvm)                                :: firedeath 
    ! What fraction of a burned plant compartment goes into the atmosphere
    !   (rest into litter)
    REAL(r_std), SAVE, DIMENSION(nparts)                             :: co2frac
    ! Moisture limit (critical moisture limit)
    REAL(r_std), DIMENSION(npts)                                     :: moistlimit
    ! total litter above the ground for a vegetation type (gC/m**2)
    REAL(r_std), DIMENSION(npts)                                     :: litter_above
    ! daily fire index
    REAL(r_std), DIMENSION(npts,nvm)                           :: fireindex_daily
    ! fire extent, on ground
    REAL(r_std), DIMENSION(npts, nvm)                          :: firefrac
    ! residual fraction of exposed structural litter, depending on lignin fraction
    REAL(r_std), DIMENSION(npts)                                     :: struc_residual
    ! residue, i.e. exposed carbon - volatilized carbon ( gC/(m**2 of ground))
    REAL(r_std), DIMENSION(npts)                                     :: residue
    ! fraction of residue transformed into black carbon
    REAL(r_std), DIMENSION(npts)                                     :: bcfrac
    ! intermediate variable
    REAL(r_std), DIMENSION(npts)                                     :: x
    ! annual fire fraction
    REAL(r_std), DIMENSION(npts)                                     :: aff
    ! index
    INTEGER(i_std)                                                  :: j,k,m

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering fire'

    !
    ! 1 Initializations
    !

    IF ( firstcall ) THEN

       !
       ! 1.1 What fraction of the plant compartment, if burned, is transformed into CO2
       !

       co2frac(ileaf) = .95
       co2frac(isapabove) = .95
       co2frac(isapbelow) = 0.
       co2frac(iheartabove) = 0.3
       co2frac(iheartbelow) = 0.
       co2frac(iroot) = 0.
       co2frac(ifruit) = .95
       co2frac(icarbres) = 0.95

       !
       ! 1.2 messages
       !

       WRITE(numout,*) 'fire:'

       WRITE(numout,*) '   > temporal inertia of fire probability (d): ',tau_fire

       WRITE(numout,*) '   > fraction of burned biomass that becomes CO2:'
       WRITE(numout,*) '     leaves: ', co2frac(ileaf)
       WRITE(numout,*) '     sap above ground: ', co2frac(isapabove)
       WRITE(numout,*) '     sap below ground: ', co2frac(isapbelow)
       WRITE(numout,*) '     heartwood above ground: ', co2frac(iheartabove)
       WRITE(numout,*) '     heartwood below ground: ', co2frac(iheartbelow)
       WRITE(numout,*) '     roots: ', co2frac(iroot)
       WRITE(numout,*) '     fruits: ', co2frac(ifruit)
       WRITE(numout,*) '     carbohydrate reserve: ', co2frac(icarbres)

       WRITE(numout,*) '   > critical litter quantity (gC/m**2): ',litter_crit
       WRITE(numout,*) '   > We calculate a fire probability on agricultural ground, but'
       WRITE(numout,*) '       the effective fire fraction is zero.'

       firstcall = .FALSE.
       !
       ! 1.3 read the flag that disable fire
       !
       !Config  Key  = FIRE_DISABLE
       !Config  Desc = no fire allowed
       !Config  Def  = n
       !Config  Help = With this variable, you can allow or not
       !Config         the estimation of CO2 lost by fire
       !
       disable_fire=.FALSE.
       CALL getin_p('FIRE_DISABLE', disable_fire)
    ENDIF

    !
    ! 1.4 Initialize output
    !

    co2_fire(:,:) = zero
    firedeath(:,:) = zero
    fireindex_daily(:,:) = zero
    firefrac(:,:) = zero

    !
    ! 2 Determine fire probability. We calculate this probability (and long-term litter)
    !   also for agricultural ground, but the fire fraction on agricultural ground is set
    !   to 0 for the moment.
    !


    DO j = 2,nvm

       ! total litter above the ground, for the vegetation type we are talking about

       litter_above(:) = litter(:,imetabolic,j,iabove) + &
            litter(:,istructural,j,iabove)

       !
       ! 2.1 calculate moisture limit. If it stays 0, this means that there is no litter
       !     on the ground, and this means that there can be no fore.
       !     Sum over different litter parts from the various PFTs, taking into account the 
       !       litter flamability which is a function of the PFT. 
       !     Difference to Stephen Sitch's approach: Daily litter, not annual mean.
       !     Reason: 1. seems more reasonable.
       !             2. easier to implement (otherwise, would need moisture limit
       !                from previous year)
       !

       moistlimit(:) = zero

       WHERE ( (t2m_daily(:) .GT. ZeroCelsius) .AND. (litter_above(:) .GT. min_stomate) )

          moistlimit(:) = &
               ( litterpart(:,j,imetabolic)*litter(:,imetabolic,j,iabove) + &
               litterpart(:,j,istructural)*litter(:,istructural,j,iabove)  ) / &
               litter_above(:) * flam(j)

       ENDWHERE

       !
       ! 2.2 daily fire index.
       !
       WHERE ( moistlimit(:) .GT. min_stomate )

          ! is a function of litter humidity. Very sensible to STOMATE's time step as
          ! with larger dt, one misses dry days with very high fireindex ( strongly
          ! nonlinear: exp(-x^2)! )

          x(:) = litterhum_daily(:)/moistlimit(:)
          fireindex_daily(:,j) = EXP( - pi * x(:) * x(:) )

       ELSEWHERE

          fireindex_daily(:,j) = zero

       ENDWHERE

       !
       ! 2.3 increase long-term fire index (mean probability)
       !

       fireindex(:,j) = ((tau_fire - dt) * fireindex(:,j) + (dt) * fireindex_daily(:,j)) / tau_fire

       !
       ! 2.4 litter influences fire intensity.
       !     We use longer-term litter to be consistent with the fire index.
       !

       firelitter(:,j) = &
            ( ( tau_fire-dt ) * firelitter(:,j) + dt * litter_above(:) ) / tau_fire


       !
       ! 3 Calculate fire fraction from litter and fireindex (i.e. basically drought)

       !
       ! 3.1 natural ground
       !

       !     This formulation has been developped for annual fire indices!
       !     original form: firefrac(i) = fireindex(i) * EXP( f(fireindex(i)) )
       !     Transform into daily fire fraction.

       ! annual fire fraction

       aff(:) = firefrac_func (npts, fireindex(:,j))

       ! transform from annual fraction to daily fraction. 
       ! annual fire fraction = 1. - (fraction of tree that survives each day) ** 365 =
       !                      = 1. - ( 1. - daily fire fraction )**365
       ! Thus, daily fire fraction = 1. - ( 1. - annual fire fraction )**(1/365)
       ! If annual firefrac<<1, then firefrac_daily = firefrac * dt/one_year
       ! This approximation avoids numerical problems.

       IF(.NOT.disable_fire.AND.natural(j))THEN
          WHERE ( aff(:) .GT. 0.1 )
             firefrac(:,j) = un - ( un - aff(:) ) ** (dt/one_year)
          ELSEWHERE
             firefrac(:,j) = aff(:) * dt/one_year
          ENDWHERE
       ELSE
          firefrac(:,j) = zero
       ENDIF

       ! No fire if litter is below critical value

       WHERE ( firelitter(:,j) .LT. litter_crit )
          firefrac(:,j) = zero
       ENDWHERE

       ! However, there is a minimum fire extent

       firefrac(:,j) = MAX( 0.001_r_std * dt/one_year, firefrac(:,j) )

       ! if FIRE_DISABLE flag is set no fire
       IF (disable_fire) firefrac=0
       !
       ! 3.2 agricultural ground: no fire for the moment
       !

       IF ( .NOT. natural(j)) firefrac(:,j) = zero

       !
       ! 4 Determine fire impact: calculate fire disturbance for each PFT
       !

       !
       ! 4.1 are we talking about a natural or an agricultural PFT?
       !

       !
       ! 4.2 fire disturbance
       !

       IF ( tree(j) ) THEN

          ! 4.2.1 Trees: always disturbed

          fire_disturb(:) = ( un - resist(j) ) * firefrac(:,j)

       ELSE

          ! 4.2.2 Grasses are not disturbed if they are not in their growing season

          WHERE ( biomass(:,j,ileaf) .GT. min_stomate )

             fire_disturb(:) = ( un - resist(j) ) * firefrac(:,j)

          ELSEWHERE

             fire_disturb(:) = zero

          ENDWHERE

       ENDIF

       !
       ! 4.3 litter and co2 created through fire on living biomass
       !

       ! biomass can go into litter or atmosphere, depending on what plant compartment
       ! we are talking about.

       DO k = 1, nparts

          ! grass roots and reserve survive.

          IF ( .NOT. ( ( .NOT. tree(j) ) .AND. ( ( k.EQ.iroot ) .OR. ( k.EQ.icarbres ) ) ) ) THEN

             ! 4.3.1 A fraction goes directly into the atmosphere.
             !       CO2 flux in gC/m**2 of total ground/day.
             !NV passe en 2D            
             co2_fire(:,j) =  co2_fire(:,j)+ biomass(:,j,k) * fire_disturb(:) * co2frac(k) / dt

             ! 4.3.2 Determine the residue, in gC/m**2 of ground.

             residue(:) = biomass(:,j,k) * fire_disturb(:) * ( un - co2frac(k) )
             !MM in SZ ???        residue(:) = fire_disturb(:) * ( un - co2frac(k) )

             ! 4.3.2.1 determine fraction of black carbon. Only for plant parts above the
             !         ground, i.e. when co2_frac > 0.
             !         A small part of the residue, which can be expressed as a function of
             !         the fraction of volatilized carbon (assimilated to co2frac here),
             !         becomes black carbon, and thus withdrawn from the soil carbon cycle (added
             !         to the "geologic carbon cycle we don't care about here).
             !         [Kuhlbusch et al. JGR 101, 23651-23665, 1996; Kuhlbusch & Crutzen, GBC 9,
             !          491-501, 1995].

             !             IF ( co2frac(k) .GT. min_stomate ) THEN
             IF ( co2frac(k) .GT. zero ) THEN
                bcfrac(:) = .3 / ( 1.3 ** ( 88.2 - 100.*co2frac(k) ) + 1. )
             ELSE
                bcfrac(:) = zero
             ENDIF

             ! 4.3.2.2 Add this fraction of the residue to the black carbon "reservoir", in
             !         gC/(m**2 of total ground).

             black_carbon(:) = &
                  black_carbon(:) + bcfrac(:) * residue(:)
             !MM in SZ                  black_carbon(:) + bcfrac(:) * residue(:) * biomass(:,j,k) 

             ! 4.3.2.3 The rest (largest part) of the residue becomes litter.
             bm_to_litter(:,j,k) = bm_to_litter(:,j,k) + residue(:) * ( 1 - bcfrac(:) )
             !MM in SZ   bm_to_litter(:,j,k) = bm_to_litter(:,j,k) + residue(:) * ( 1 - bcfrac(:) ) * biomass(:,j,k) 
          ENDIF  ! not for grass roots

       ENDDO

       !
       ! 4.4 new vegetation characteristics
       !

       ! 4.4.1 decrease biomass per m**2 of ground
       !       except for grass roots.

       DO k = 1, nparts

          IF ( .NOT. ( ( .NOT. tree(j) ) .AND. ( ( k.EQ.iroot ) .OR. ( k.EQ.icarbres) ) ) ) THEN

             biomass(:,j,k) = ( un - fire_disturb(:) ) * biomass(:,j,k)

          ENDIF

       ENDDO

       ! 4.4.2 If vegetation is dynamic, then decrease the density of tree
       !       individuals.

       IF ( control%ok_dgvm .AND. tree(j) ) THEN

          ! fraction of plants that dies each day.
          ! exact formulation: 1. - (1.-fire_disturb(:)) ** (1./dt)
          firedeath(:,j) = fire_disturb(:) / dt

          ind(:,j) = ( un - fire_disturb(:) ) * ind(:,j)

       ENDIF

    ENDDO      ! loop over PFTs

    !
    ! 5 A fraction of the litter is burned by the fire
    !

    DO j = 2,nvm

       !
       ! 5.1 exposed metabolic litter burns totally and goes directly into the atmosphere as
       !     CO2.
       !

       ! 5.1.1 CO2 flux, in gC/(m**2 of total ground)/day.

       co2_fire(:,j) = co2_fire(:,j) + litter(:,imetabolic,j,iabove) * &
            firefrac(:,j) / dt

       ! 5.1.2 decrease metabolic litter

       litter(:,imetabolic,j,iabove) = litter(:,imetabolic,j,iabove) * &
            ( un - firefrac(:,j) )

       !
       ! 5.2 exposed structural litter is not totally transformed into CO2.
       !

       ! 5.2.1 Fraction of exposed structural litter that does not
       !       burn totally should depend on lignin content (lignin_struc). VERY TENTATIVE!

       struc_residual(:) = 0.5 * lignin_struc(:,j,iabove)

       ! 5.2.2 CO2 flux, in gC/(m**2 of total ground)/day.

       co2_fire(:,j) = co2_fire(:,j) + &
            litter(:,istructural,j,iabove) * firefrac(:,j) * &
            ( un - struc_residual(:) )/ dt

       ! 5.2.3 determine residue (litter that undergoes fire, but is not transformed
       !       into CO2)

       residue(:) = litter(:,istructural,j,iabove) * firefrac(:,j) * &
            struc_residual(:)
       !MM in SZ        residue(:) = firefrac(:,j) * struc_residual(:)

       ! 5.2.4 determine fraction of black carbon in the residue.
       !       depends on volatilized fraction of carbon (see 4.3.2.1)

       bcfrac(:) = .3 / ( 1.3 ** ( 88.2 - 100.*(1.-struc_residual(:)) ) + 1. )

       ! 5.2.5 Add this fraction of the residue to the black carbon "reservoir", in
       !       gC/(m**2 of total ground).

       black_carbon(:) = &
            black_carbon(:) + bcfrac(:) * residue(:)
       !MM in SZ             black_carbon(:) + bcfrac(:) * residue(:) * litter(:,iwoody,j,iabove) 

       ! 5.2.6 The rest (largest part) of the residue remains litter. Remaining litter
       !       is the sum of this and of the litter which has not undergone a fire.

       litter(:,istructural,j,iabove) = &
            litter(:,istructural,j,iabove) * ( un - firefrac(:,j) ) + &
            residue(:) * ( un - bcfrac(:) )
       !MM in SZ            residue(:) * ( un - bcfrac(:) ) * litter(:,iwoody,j,iabove)

    ENDDO  !  ground

    !
    ! 5.3 diagnose fraction of leaves burned.
    !     exposed leaves are burned entirely, even their structural part
    !

    DO j = 2,nvm

       DO k = 1, nlitt
          dead_leaves(:,j,k) = dead_leaves(:,j,k) * ( un - firefrac(:,j) )
       ENDDO

    ENDDO

    !
    ! 6 history
    !

    ! output in 1/day
    firefrac(:,:) = firefrac(:,:) / dt

    CALL histwrite (hist_id_stomate, 'FIREFRAC', itime, &
         firefrac(:,:), npts*nvm, horipft_index)
    CALL histwrite (hist_id_stomate, 'FIREDEATH', itime, &
         firedeath(:,:), npts*nvm, horipft_index)

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving fire'

  END SUBROUTINE fire

  FUNCTION firefrac_func (npts, x) RESULT (firefrac_result)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                             :: npts
    ! fire index
    REAL(r_std), DIMENSION(npts), INTENT(in)                         :: x

    ! 0.2 result

    ! fire fraction
    REAL(r_std), DIMENSION(npts)                                     :: firefrac_result

    ! 0.3 local

    ! intermediate variable
    REAL(r_std), DIMENSION(npts)                                     :: xm1

    xm1(:) = x(:) - 1.

    firefrac_result(:) = &
         x(:) * EXP( xm1(:) / ( -.13*xm1(:)*xm1(:)*xm1(:) + .6*xm1(:)*xm1(:) + .8*xm1(:) + .45 ) )

  END FUNCTION firefrac_func

END MODULE lpj_fire