! 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_data USE parallel USE pft_parameters USE constantes 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 ! 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 ! 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 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) = 1. - ( 1. - 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(:) = ( 1. - 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(:) = ( 1. - 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(:) * ( 1. - co2frac(k) ) !MM in SZ ??? residue(:) = fire_disturb(:) * ( 1. - 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(:) = bcfrac_coeff(1) / ( bcfrac_coeff(2) ** ( bcfrac_coeff(3) - 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) = ( 1. - 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) = ( 1. - 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) * & ( 1. - 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) * & ( 1. - 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(:) = bcfrac_coeff(1)/ ( bcfrac_coeff(2) ** ( bcfrac_coeff(3) - 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) * ( 1. - firefrac(:,j) ) + & residue(:) * ( 1. - bcfrac(:) ) !MM in SZ residue(:) * ( 1. - 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) * ( 1. - 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 ) ) x(:) * EXP( xm1(:) / ( -firefrac_coeff(4)*xm1(:)*xm1(:)*xm1(:) + firefrac_coeff(3)*xm1(:)*xm1(:) + firefrac_coeff(2)*xm1(:) + firefrac_coeff(1) ) ) END FUNCTION firefrac_func END MODULE lpj_fire