source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/lpj_fire.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 28.9 KB
Line 
1! =================================================================================================================================
2! MODULE        : lpj_fire
3!
4! CONTACT       : orchidee-help _at_ listes.ipsl.fr
5!
6! LICENCE       : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
7!
8!>\BRIEF         Calculates the probability of fire and its effects on the carbon cycle
9!!
10!!\n DESCRIPTION: None
11!!
12!! RECENT CHANGE(S): None
13!!
14!! SVN           :
15!! $HeadURL$
16!! $Date$
17!! $Revision$
18!! \n
19!_ ================================================================================================================================
20
21MODULE lpj_fire
22
23  ! modules used:
24  USE xios_orchidee
25  USE stomate_data
26  USE ioipsl_para 
27  USE pft_parameters
28  USE constantes
29
30  IMPLICIT NONE
31
32  ! private & public routines
33
34  PRIVATE
35  PUBLIC fire,fire_clear
36
37  LOGICAL, SAVE                   :: firstcall_fire = .TRUE.        !! first call
38!$OMP THREADPRIVATE(firstcall_fire)
39
40CONTAINS
41
42
43!! ================================================================================================================================
44!!  SUBROUTINE   : fire_clear
45!!
46!>\BRIEF        Set the firstcall_fire flag to .TRUE. and activate initialization
47!!
48!_ ================================================================================================================================
49
50  SUBROUTINE fire_clear
51    firstcall_fire = .TRUE.
52  END SUBROUTINE fire_clear
53
54
55!! ================================================================================================================================
56!! SUBROUTINE     : fire
57!!
58!>\BRIEF          Calculate fire index and fraction of area burned by comparing
59!! daily litter moisture with prescribed moisture of fire inhibition. If daily
60!! moisture is below the precribed threshold, allow fire disturbance and
61!! calculate CO_2 emissions from fire. The main algorithm follows Thonicke et al.
62!! (2001)
63!!
64!! DESCRIPTION  : Fire occurs when the three basic prerequisites are met:
65!! sufficient fuel on the ground, relatively dry fuel (fuel moisture lower than
66!! the moisture of extinction), and presence of ignition sources.  The flag
67!! ::disable_fire is used to use this fire module or not. While here in the module the
68!! ignition source is not explicitely represented. It's assumed
69!! that the ignition source is always available, i.e., if a sufficient amount
70!! of dead fuel exists with a moisture content below the moisture of fire
71!! inhibition, then both live and dead fuel will start to burn.\n
72!!
73!! The module completes the following tasks:
74!! 1. Calculates daily fire index, long term fire index and available ground
75!! litter for burning.\n
76!! \latexonly
77!! \input{lpj_fire_1.tex} 
78!! \endlatexonly
79!! Where, m is the daily litter moisture, m_e is the moisture of extinction, p(m)
80!! is probability of fire (i.e. the daily fire index).\n
81!!
82!! \latexonly
83!! \input{lpj_fire_2.tex} 
84!! \endlatexonly
85!! Where, s is the long term fire index and A(s) is the annual fire fraction.\n
86!! 2. Calculates annual fire fraction, then transform to daily fraction in the
87!! same time step of stomate.\n
88!! 3. Ground litter and grass live biomass (in growing season and except root and carbon
89!! reserve pool) are  burned at full fire fraction.\n
90!! 4. Fire fraction for tree PFT are compromised by prescribed fire resistence.
91!! Tree live biomass are consumed using this compromised fire fraction and consumption
92!! fraction for different biomass parts. In the case of activation of dynamic
93!! vegetation, tree individual density is updated after fire. \n
94!! 5. For all types of fuel (either ground litter or live biomass) that are
95!! burned, there is a certain fraction of “residue” that are not completely
96!! burned. The remaining part is transferred to (in case of biomass) or
97!! remains (in case of ground litter) as litter.\n
98!! 6. Update the biomass, litter or dead leaves pool after burning.
99!! 7. If the flag SPINUP_ANALYTIC is set to true, the matrix A is updated following
100!! Lardy (2011).
101!!
102!! RECENT CHANGE(S): April 2015: Black carbon calculations is removed. The black carbon was not conserved.
103!!
104!! MAIN OUTPUT VARIABLE(S): ::co2_fire or the carbon emitted into the atmosphere
105!! by fire, including both living and dead biomass (gC m^{-2} dtslow^{-1})$ @endtex
106!!
107!! REFERENCES    :
108!! - Thonicke K., Venevsky S., Sitch S., and Cramer W. (2001) The role of fire
109!! disturbance for global vegetation dynamics: coupling fire into a Dynamic
110!! Global Vegetation Model, Global Ecology & Biogeography, 10, 661-677.
111!! - Kuhlbusch et al. JGR 101, 23651-23665, 1996
112!! - Kuhlbusch & Crutzen, GBC 9, 491-501, 1995
113!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
114!! Environmental Modelling & Software (2011), doi:10.1016|j.envsoft.2011.05.016
115!!
116!! FLOWCHART     : None
117!! \n
118!_ ================================================================================================================================
119 
120  SUBROUTINE fire (npts, dt, litterpart, &
121       litterhum_daily, t2m_daily, lignin_struc,veget_cov_max, &
122       fireindex, firelitter, biomass, ind, &
123       litter, dead_leaves, bm_to_litter, &
124       co2_fire, MatrixA)
125
126  !! 0. Variable and parameter declarations
127
128    !! 0.1 Input variables
129   
130    INTEGER(i_std), INTENT(in)                                   :: npts                !! Domain size - number of pixels
131                                                                                        !! (unitless)   
132    REAL(r_std), INTENT(in)                                      :: dt                  !! Time step of the simulations for stomate
133                                                                                        !! (days)   
134    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in)           :: litterpart          !! [DISPENSABLE] Fraction of litter above
135                                                                                        !! the ground belonging to different PFTs
136                                                                                        !! ?? this variable is used but might be
137                                                                                        !! redundant (with value of always 1) ??
138                                                                                        !! To check but probably litterpart to be
139                                                                                        !! removed. Probably a residual from
140                                                                                        !! nat/agri before merge   
141    REAL(r_std), DIMENSION(npts), INTENT(in)                     :: litterhum_daily     !! Daily litter moisture (unitless, 0-1)   
142    REAL(r_std), DIMENSION(npts), INTENT(in)                     :: t2m_daily           !! Daily 2 meter temperature (K)   
143    REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(in)           :: lignin_struc        !! Ratio Lignin/Carbon in structural above
144                                                                                        !! and below ground litter
145                                                                                        !! @tex $(gC gC^{-1})$ @endtex
146    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                 :: veget_cov_max       !! Maximum fraction of vegetation type including
147                                                                                        !! non-biological fraction (0-1, unitless)
148
149    !! 0.2 Output variables
150   
151    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                :: co2_fire            !! Carbon emitted into the atmosphere by
152                                                                                        !! fire, including both living and dead
153                                                                                        !! biomass
154                                                                                        !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
155    !! 0.3 Modified variables
156   
157    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)              :: fireindex           !! Probability of fire; in the code means
158                                                                                        !! long term fire index (unitless, 0-1)   
159    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)              :: firelitter          !! Total natural litter
160                                                                                        !! available for burning above the ground
161                                                                                        !! @tex $(gC m^{-2})$ @endtex
162    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout)  :: biomass        !! Biomass @tex $(gC m^{-2})$ @endtex
163    REAl(r_std), DIMENSION(npts,nvm), INTENT(inout)                   :: ind            !! Density of individuals
164                                                                                        !! @tex $(m^{-2})$ @endtex
165    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter     !! Metabolic and structural litter, above
166                                                                                        !! and below ground
167                                                                                        !! @tex $(gC m^{-2})$ @endtex
168    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)        :: dead_leaves         !! Dead leaves on ground, per PFT,
169                                                                                        !! metabolic and structural
170                                                                                        !! @tex $(gC m^{-2})$ @endtex
171    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: bm_to_litter     !! Biomass entering the litter pool
172                                                                                        !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
173    REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA          !! Matrix containing the fluxes between the carbon pools
174                                                                                        !! here, it is called only if disable_fire = n
175                                                                                        !! once a day
176                                                                                        !! @tex $(gC.m^2.day^{-1})$ @endtex
177   
178
179    !! 0.4 Local variables
180
181    REAL(r_std), DIMENSION(npts)                                  :: fire_disturb       !! Actual fire disturbance fraction
182                                                                                        !! after consideration of inherent fire
183                                                                                        !! resistance of different PFTs
184                                                                                        !! (unitless, 0-1) 
185    REAL(r_std), DIMENSION(npts,nvm)                              :: firedeath          !! In case of activation of dynamic
186                                                                                        !! vegetation, the daily fraction of
187                                                                                        !! burned individuals (unitless, 0-1) 
188    REAL(r_std), DIMENSION(npts)                                  :: moistlimit         !! Moisture for fire inhibition per PFT                                                                                         !!
189                                                                                        !! (unitless, 0-1) ; temporary variable
190                                                                                        !! for each PFT in the loop over #PFT ?
191                                                                                        !! it's not fire inhibition, it's
192                                                                                        !! "moisture of extinction" by the original
193                                                                                        !! reference paper, and a more "fire
194                                                                                        !! professional name" it's better to change
195                                                                                        !! the name of moistlimit
196    REAL(r_std), DIMENSION(npts)                                  :: litter_above       !! Total aboveground litter per PFT                                                                                         !!
197                                                                                        !! @tex $(gC m^{-2})$ @endtex
198    REAL(r_std), DIMENSION(npts,nvm)                              :: fireindex_daily    !! Daily fire index (unitless, 0-1)
199    REAL(r_std), DIMENSION(npts, nvm)                             :: firefrac           !! Daily fire burning fraction on ground
200                                                                                        !! (unitless, 0-1)
201    REAL(r_std), DIMENSION(npts)                                  :: struc_residual     !! Fraction of structural litter not burned
202                                                                                        !! (thus residual), depending on strutural
203                                                                                        !! litter lignin content  (unitless, 0-1)
204    REAL(r_std), DIMENSION(npts)                                  :: residue            !! Fuel (either biomass or litter) not
205                                                                                        !! burned  @tex $(gC m^{-2})$ @endtex
206    REAL(r_std), DIMENSION(npts)                                  :: x                  !! Intermediate variable
207    REAL(r_std), DIMENSION(npts)                                  :: aff                !! Annual fire fraction (unitless, 0-1)
208    INTEGER(i_std)                                                :: j,k,m              !! Indeces
209!_ ================================================================================================================================
210
211    IF (printlev>=3) WRITE(numout,*) 'Entering fire'
212
213 !! 1. Initialization
214
215    IF ( firstcall_fire ) THEN
216
217       !! 1.1 Fraction to CO_2
218       !      What fraction of the plant biomass compartment, if burned, is transformed
219       !      into CO_2 released to the atmosphere?
220       !??! the fraction for heartabove seems too big, it's not clear if this value is correct.
221
222       !! 1.2 print control messages
223       IF (printlev >= 2) THEN
224          WRITE(numout,*) 'fire:'
225          WRITE(numout,*) '   > temporal inertia of fire probability (d): ',tau_fire
226         
227          WRITE(numout,*) '   > fraction of burned biomass that becomes CO2:'
228          WRITE(numout,*) '     leaves: ', co2frac(ileaf)
229          WRITE(numout,*) '     sap above ground: ', co2frac(isapabove)
230          WRITE(numout,*) '     sap below ground: ', co2frac(isapbelow)
231          WRITE(numout,*) '     heartwood above ground: ', co2frac(iheartabove)
232          WRITE(numout,*) '     heartwood below ground: ', co2frac(iheartbelow)
233          WRITE(numout,*) '     roots: ', co2frac(iroot)
234          WRITE(numout,*) '     fruits: ', co2frac(ifruit)
235          WRITE(numout,*) '     carbohydrate reserve: ', co2frac(icarbres)
236         
237          WRITE(numout,*) '   > critical litter quantity (gC/m**2): ',litter_crit
238          WRITE(numout,*) '   > We calculate a fire probability on agricultural ground, but'
239          WRITE(numout,*) '       the effective fire fraction is zero.'
240       END IF
241       firstcall_fire = .FALSE.
242
243    ENDIF
244
245    !! 1.4 Initialize output
246    co2_fire(:,:) = zero
247    firedeath(:,:) = zero
248    fireindex_daily(:,:) = zero
249    firefrac(:,:) = zero
250
251 !! 2. Determine fire probability
252 
253    ! Calculate probability that crops are not burned for the moment
254    ! Calculate long-term aboveground litter
255 
256    ! a long loop over PFT
257    DO j = 2,nvm !loop over #PFT
258
259       !! 2.1 Total above ground litter
260       !      Total litter above the ground for the PFT under consideration
261       litter_above(:) = litter(:,imetabolic,j,iabove,icarbon) + &
262            litter(:,istructural,j,iabove,icarbon)
263
264
265       !! 2.2 Calculate the PFT litter amount weighted by moisture
266       ! If litter moisture is higher than moisture of extinction, fire is not possible.
267       moistlimit(:) = zero
268
269       WHERE ( (t2m_daily(:) .GT. ZeroCelsius) .AND. (litter_above(:) .GT. min_stomate) )
270       !?? the calculation here is redundant? as the part before flam(j) is 1?
271       !?? see first comment.. litterpart to be removed ?
272          moistlimit(:) = &
273               ( litterpart(:,j,imetabolic)*litter(:,imetabolic,j,iabove,icarbon) + &
274               litterpart(:,j,istructural)*litter(:,istructural,j,iabove,icarbon)  ) / &
275               litter_above(:) * flam(j)
276
277       ENDWHERE
278
279       !! 2.3 Calculate daily fire index
280       !      Calculate daily fire index
281       !      \latexonly
282       !        \input{lpj_fire_1.tex} 
283       !      \endlatexonly
284       !      Where, m is the daily litter moisture, m_e is the moisture of extinction, p(m)
285       !      is probability of fire (i.e. the daily fire index).\n
286       WHERE ( moistlimit(:) .GT. min_stomate )
287          x(:) = litterhum_daily(:)/moistlimit(:)
288          fireindex_daily(:,j) = EXP( - pi * x(:) * x(:) )
289       ELSEWHERE
290          fireindex_daily(:,j) = zero
291       ENDWHERE
292
293       !! 2.4 Calculate long-term fire index
294       !      Calculate long-term fire index which is the mean probability of fire
295       fireindex(:,j) = ((tau_fire - dt) * fireindex(:,j) + (dt) * fireindex_daily(:,j)) / tau_fire
296
297       !! 2.5 Calculate long term aboveground litter that are available for burning
298       firelitter(:,j) = &
299            ( ( tau_fire-dt ) * firelitter(:,j) + dt * litter_above(:) ) / tau_fire
300
301  !! 3. Calculate fire fraction from litter and fireindex
302
303       !! 3.1 Fire fraction from long term fire index for natural PFTs
304       !  Transform the annual fire fraction to daily fraction. This is done by assuming that
305       !  each day the fire fraction is the same.
306       aff(:) = firefrac_func (npts, fireindex(:,j))
307
308       ! annual_fire_fraction = 1. - ( 1. - daily_fire_fraction )**365
309       ! Thus, daily_fire_fraction = 1. - ( 1. - annual_fire_fraction )**(1/365)
310       ! If annual firefrac<<1, then firefrac_daily = firefrac * dt/one_year
311       ! This approximation avoids numerical problems.
312       ! the variable 'un' is 1
313       IF(.NOT.disable_fire.AND.natural(j))THEN
314          WHERE ( aff(:) .GT. 0.1 )
315             firefrac(:,j) = un - ( un - aff(:) ) ** (dt/one_year)
316          ELSEWHERE
317             firefrac(:,j) = aff(:) * dt/one_year
318          ENDWHERE
319       ELSE
320          firefrac(:,j) = zero
321       ENDIF
322
323       ! No fire if litter is below critical value
324       WHERE ( firelitter(:,j) .LT. litter_crit )
325          firefrac(:,j) = zero
326       ENDWHERE
327
328       ! However, there is a minimum fire extent
329       firefrac(:,j) = MAX( 0.001_r_std * dt/one_year, firefrac(:,j) )
330
331       ! if FIRE_DISABLE flag is set no fire
332       IF (disable_fire) firefrac = zero
333       
334       !! 3.2 For agricultural ground, no fire is burned for the moment
335       IF ( .NOT. natural(j)) firefrac(:,j) = zero
336       
337  !! 4. Determine fire impact
338
339       ! Calculate fire disturbance fraction for each PFT, and fire emissions due
340       ! to grasses.
341   
342       !! 4.1 Tree and grass live biomass
343       ! Tree live biomass is not burned in fire.
344       ! However, in the dynamic vegetation module, tree individual density
345       ! will be updated after fire. The fraction of tree individuals that are
346       ! supposed to die in fire is the fire fraction multiplied by the tree PFT fire
347       ! resistance which reflect survivorship of different tree PFTs during fire. 
348       IF ( is_tree(j) ) THEN
349
350          !! 4.1.1 Disturban,ce factor for trees
351          !        Trees are disturbed over the whole year. A resistance factor is 
352          !        used to reflect survivorship of different tree PFTs during fire.
353          fire_disturb(:) = ( un - resist(j) ) * firefrac(:,j)
354
355       ELSE
356
357          !! 4.1.2 Disturbance factor for grasses
358          !        Grass is not disturbed outside the growing season thus grass biomass
359          !        is only burned during the growing season.
360          WHERE ( biomass(:,j,ileaf,icarbon) .GT. min_stomate )
361
362             fire_disturb(:) = ( un - resist(j) ) * firefrac(:,j)
363
364          ELSEWHERE
365
366             fire_disturb(:) = zero
367
368          ENDWHERE
369
370       ENDIF
371
372       !! 4.2 Burn live biomass
373       !      The burned part goes to either the atmposhere or litter
374       DO k = 1, nparts
375
376          ! for tree PFT, all biomass compartments are burned.
377          ! for grass biomass compartments, all are burned except root and carbon reserve
378          ! IF concerning PFT is tree; OR (the PFT is grass, but not root or carbon reserve biomass); then it's burned.
379          IF ( .NOT. ( ( .NOT. is_tree(j) ) .AND. ( ( k.EQ.iroot ) .OR. ( k.EQ.icarbres ) ) ) ) THEN
380
381             !! 4.2.1 Fraction to the atmosphere.
382             co2_fire(:,j) =  co2_fire(:,j)+ biomass(:,j,k,icarbon) * fire_disturb(:) * co2frac(k) / dt
383
384             !! 4.2.2 Determine the residue
385             !        Residual is expressed in gC/m^2
386             residue(:) = biomass(:,j,k,icarbon) * fire_disturb(:) * ( un - co2frac(k) )
387
388             !! 4.2.2.3 The rest (largest part) of the residue becomes litter.
389             bm_to_litter(:,j,k,icarbon) = bm_to_litter(:,j,k,icarbon) + residue(:)
390          ENDIF
391
392       ENDDO
393
394     
395       !! 4.3 Update biomass pool after burning
396       
397       !! 4.3.1 Decrease biomass amount except for grass roots and carbon reserve.
398       DO m = 1, nelements ! Loop over # elements
399         
400          DO k = 1, nparts
401
402             ! **2
403             IF ( .NOT. ( ( .NOT. is_tree(j) ) .AND. ( ( k.EQ.iroot ) .OR. ( k.EQ.icarbres) ) ) ) THEN
404               
405                biomass(:,j,k,m) = ( un - fire_disturb(:) ) * biomass(:,j,k,m)
406               
407             ENDIF
408
409          ENDDO
410         
411       END DO ! End Loop over # elements
412
413
414       !! 4.3.2 If vegetation is dynamic, then decrease the density of tree individuals.
415       IF ( (ok_dgvm .OR. .NOT.lpj_gap_const_mort) .AND. is_tree(j) ) THEN
416
417          firedeath(:,j) = fire_disturb(:) / dt
418
419          ind(:,j) = ( un - fire_disturb(:) ) * ind(:,j)
420
421       ENDIF
422
423    ENDDO      ! loop over #PFT
424
425  !! 5. Burn litter
426
427    !   All litter (either forest or grass) is burned throughout the whole year
428    !   Metabolic litter is totally burned. Burning fraction of structural litter is related
429    !   with its lignin content. The burned part either goes to the atmosphere
430    !   or remains in litter as unburned residue.
431    DO j = 2,nvm  !loop over #PFT
432
433       !! 5.1 Burn exposed metabolic litter
434       !      Exposed metabolic litter burns completely and goes directly into the atmosphere as CO2.
435
436       !! 5.1.1 CO2 flux from litter burning
437       !        Flux expressed in gC m^{-2}dtslow^{-1}
438       co2_fire(:,j) = co2_fire(:,j) + litter(:,imetabolic,j,iabove,icarbon) * &
439            firefrac(:,j) / dt
440
441       !! 5.1.2 Decrease metabolic litter
442       DO m = 1,nelements
443          litter(:,imetabolic,j,iabove,m) = litter(:,imetabolic,j,iabove,m) * &
444               ( un - firefrac(:,j) )
445       END DO
446   
447       !! 5.2 Burning exposed structural litter
448       !      The fraction that goes to the atmosphere is related to its lignin content. The remaining unburned
449       !      residues remain as litter.
450
451       !! 5.2.1 Incomplete burning of exposed structural litter 
452       !        Fraction of structural litter that does not burn completly. This fraction depends on lignin
453       !        content (lignin_struc).
454       struc_residual(:) = fire_resist_struct * lignin_struc(:,j,iabove)
455
456       !! 5.2.2 CO2 flux from litter burning
457       !  Flux expressed in gC m^{-2}dtslow^{-1}
458       co2_fire(:,j) = co2_fire(:,j) + &
459            litter(:,istructural,j,iabove,icarbon) * firefrac(:,j) * &
460            ( un - struc_residual(:) )/ dt
461
462       !! 5.2.3 Determine residue
463       !        The residue is litter that undergoes fire, but is not transformed into CO2
464!!       IF (ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
465!!          residue(:) = firefrac(:,j) * struc_residual(:)
466!!       ELSE
467          residue(:) = litter(:,istructural,j,iabove,icarbon) * firefrac(:,j) * &
468               struc_residual(:)
469!!       ENDIF
470
471       !! 5.2.6 TResidue remaining litter
472       !        The rest (largest part) of the residue remains litter. Remaining
473       !        litter is the sum of this and of the litter which has not undergone a fire.
474       litter(:,istructural,j,iabove,icarbon) = &
475            litter(:,istructural,j,iabove,icarbon) * ( un - firefrac(:,j) ) + &
476            residue(:)
477
478       !! 5.2.7 Update matrix A for analytical spin-up (only if SPINUP_ANALYTIC is activated)
479
480       IF (spinup_analytic) THEN
481         
482          ! litter structural above
483          MatrixA(:,j,istructural_above,istructural_above) =  MatrixA(:,j,istructural_above,istructural_above) - firefrac(:,j) &
484                                                                +  firefrac(:,j)* struc_residual(:)
485         
486          ! litter_metabolic above
487          MatrixA(:,j,imetabolic_above,imetabolic_above) = MatrixA(:,j,imetabolic_above,imetabolic_above) - firefrac(:,j)
488         
489          !DS! This is the exact formulation. The above is a simplified formulation.
490          !          MatrixA(:,j,istructural_above,istructural_above) =  MatrixA(:,j,istructural_above,istructural_above)*(1. - firefrac(:,j) + &
491          !                                                             firefrac(:,j)* struc_residual(:)
492          !          MatrixA(:,j,imetabolic_above,imetabolic_above) = MatrixA(:,j,imetabolic_above,imetabolic_above)*(1 - firefrac(:,j))
493         
494       ENDIF !(spinup_analytic)
495
496    ENDDO  !  loop over #PFT
497   
498   !! 5.3 Update exposed dead leaves (leaf litter) on the ground
499   !      Dead leaves are supposed to burn completely in fire even their structural part.
500    DO j = 2,nvm
501
502       DO k = 1, nlitt
503          dead_leaves(:,j,k) = dead_leaves(:,j,k) * ( un - firefrac(:,j) )
504       ENDDO
505
506    ENDDO
507
508  !! 6. write out history variables
509
510    ! output in 1/dtslow
511    firefrac(:,:) = firefrac(:,:) / dt
512
513    CALL xios_orchidee_send_field("firefrac",firefrac(:,:))
514    CALL xios_orchidee_send_field("firedeath",firedeath(:,:))
515
516    CALL histwrite_p (hist_id_stomate, 'FIREFRAC', itime, &
517         firefrac(:,:), npts*nvm, horipft_index)
518    CALL histwrite_p (hist_id_stomate, 'FIREDEATH', itime, &
519         firedeath(:,:), npts*nvm, horipft_index)
520
521    IF (printlev>=4) WRITE(numout,*) 'Leaving fire'
522
523  END SUBROUTINE fire
524
525
526!! ================================================================================================================================
527!! FUNCTION     : firefrac_func
528!!
529!>\BRIEF        Calculate the annual fire fraction using long term fire index
530!!
531!! DESCRIPTION  : None
532!!
533!! RECENT CHANGE(S): None
534!!
535!! RETURN VALUE : annual fire fraction (::firefrac_result)
536!!
537!! REFERENCES   :
538!! - Thonicke K., Venevsky S., Sitch S., and Cramer W. (2001) The role of fire
539!! disturbance for global vegetation dynamics: coupling fire into a Dynamic
540!! Global Vegetation Model, Global Ecology & Biogeography, 10, 661-677.\n
541!!
542!! FLOWCHART    : None
543!! \n
544!_ ================================================================================================================================
545   
546  FUNCTION firefrac_func (npts, x) RESULT (firefrac_result)
547
548!! 0. Variable and parameter declaration
549
550    !! 0.1 Input variables
551
552    INTEGER(i_std), INTENT(in)                 :: npts             !! Domain size (unitless)
553    REAL(r_std), DIMENSION(npts), INTENT(in)   :: x                !! fire index (unitless, 0-1)
554
555    !! 0.2 Output variables
556
557    REAL(r_std), DIMENSION(npts)               :: firefrac_result  !! fire fraction (unitless, 0-1)
558
559    !! 0.3 Modified variables
560
561    !! 0.4 Local variables
562    REAL(r_std), DIMENSION(npts)               :: xm1              !! intermediate variable
563!_ ================================================================================================================================
564
565!! 1. Function
566
567    xm1(:) = x(:) - 1.
568
569    ! this equation is from Thonicke et al. (2001) equation (9), it use the fire index as input to calculate annual fire fraction.
570    ! but with different parameters with the source literature.
571    !! \latexonly
572    !! \input{lpj_fire_2.tex} 
573    !! \endlatexonly
574    !! Where, s is the long term fire index and A(s) is the annual fire fraction.\n
575    !??! here we use a different parameter with K. Thonicke et al. (2001)
576    !??! it's not clear if these parameters are correct.
577    firefrac_result(:) = &
578         &     x(:) * EXP( xm1(:) / ( -firefrac_coeff(4)*xm1(:)*xm1(:)*xm1(:) + &
579         &     firefrac_coeff(3)*xm1(:)*xm1(:) + &
580         &     firefrac_coeff(2)*xm1(:) +  &     
581         &     firefrac_coeff(1) ) )
582
583  END FUNCTION firefrac_func
584
585END MODULE lpj_fire
Note: See TracBrowser for help on using the repository browser.