Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/stomate_lcchange.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: 59.5 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_lcchange
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Impact of land cover change on carbon stocks
10	!!
11	!!\n DESCRIPTION: None
12	!!
13	!! RECENT CHANGE(S): Including permafrost carbon
14	!!
15	!! REFERENCE(S) : None
16	!!
17	!! SVN :
18	!! $HeadURL$
19	!! $Date$
20	!! $Revision$
21	!! \n
22	!_ ================================================================================================================================
23
24
25	MODULE stomate_lcchange
26
27	! modules used:
28
29	USE ioipsl_para
30	USE stomate_data
31	USE pft_parameters
32	USE constantes
33	USE constantes_soil_var
34
35	IMPLICIT NONE
36
37	PRIVATE
38	PUBLIC lcchange_main
39	PUBLIC lcchange_deffire
40
41	CONTAINS
42
43
44	!! ================================================================================================================================
45	!! SUBROUTINE : lcchange_main
46	!!
47	!>\BRIEF Impact of land cover change on carbon stocks
48	!!
49	!! DESCRIPTION : This subroutine is always activate if VEGET_UPDATE>0Y in the configuration file, which means that the
50	!! vegetation map is updated regulary. lcchange_main is called from stomateLpj the first time step after the vegetation
51	!! map has been changed.
52	!! The impact of land cover change on carbon stocks is computed in this subroutine. The land cover change is written
53	!! by the difference of current and previous "maximal" coverage fraction of a PFT.
54	!! On the basis of this difference, the amount of 'new establishment'/'biomass export',
55	!! and increase/decrease of each component, are estimated.\n
56	!!
57	!! Main structure of lpj_establish.f90 is:
58	!! 1. Initialization
59	!! 2. Calculation of changes in carbon stocks and biomass by land cover change
60	!! 3. Update 10 year- and 100 year-turnover pool contents
61	!! 4. History
62	!!
63	!! RECENT CHANGE(S) : None
64	!!
65	!! MAIN OUTPUT VARIABLE(S) : ::prod10, ::prod100, ::flux10, ::flux100,
66	!! :: cflux_prod10 and :: cflux_prod100
67	!!
68	!! REFERENCES : None
69	!!
70	!! FLOWCHART :
71	!! \latexonly
72	!! \includegraphics[scale=0.5]{lcchange.png}
73	!! \endlatexonly
74	!! \n
75	!_ ================================================================================================================================
76
77
78	SUBROUTINE lcchange_main ( npts, dt_days, veget_cov_max_old, veget_cov_max_new, &
79	biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
80	co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,flux10,flux100, &
81	prod10, prod100, &
82	convflux, &
83	cflux_prod10, cflux_prod100, leaf_frac,&
84	npp_longterm, lm_lastyearmax, litter, litter_avail, litter_not_avail, &
85	carbon, &
86	deepC_a, deepC_s, deepC_p, &
87	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
88
89
90	IMPLICIT NONE
91
92	!! 0. Variable and parameter declaration
93
94	!! 0.1 Input variables
95
96	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
97	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
98	!! (days)
99	REAL(r_std), DIMENSION(nvm, nparts,nelements), INTENT(in) :: bm_sapl !! biomass of sapling
100	!! @tex ($gC individual^{-1}$) @endtex
101	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_cov_max_old!! Current "maximal" coverage fraction of a PFT (LAI
102	!! -> infinity) on ground
103	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_cov_max_new!! New "maximal" coverage fraction of a PFT (LAI ->
104	!! infinity) on ground (unitless)
105
106	!! 0.2 Output variables
107
108	REAL(r_std), DIMENSION(:,:), INTENT(out) :: convflux !! release during first year following land cover
109	!! change (npts, nwp)
110	REAL(r_std), DIMENSION(:,:), INTENT(out) :: cflux_prod10 !! total annual release from the 10 year-turnover
111	!! pool @tex ($gC m^{-2}$) @endtex
112	!! change (npts, nwp)
113	REAL(r_std), DIMENSION(:,:), INTENT(out) :: cflux_prod100 !! total annual release from the 100 year-
114	!! turnover pool @tex ($gC m^{-2}$) @endtex
115	!! change (npts, nwp)
116	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: turnover_daily !! Turnover rates
117
118
119	!! 0.3 Modified variables
120
121	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: biomass !! biomass @tex ($gC m^{-2}$) @endtex
122	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Number of individuals @tex ($m^{-2}$) @endtex
123	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! mean age (years)
124	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: senescence !! plant senescent (only for deciduous trees) Set
125	!! to .FALSE. if PFT is introduced or killed
126	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! Is pft there (unitless)
127	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or very
128	!! localized (unitless)
129	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit !! how many days ago was the beginning of the
130	!! growing season (days)
131	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm !! biomass uptaken
132	!! @tex ($gC m^{-2} day^{-1}$) @endtex
133	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! conversion of biomass to litter
134	!! @tex ($gC m^{-2} day^{-1}$) @endtex
135	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: cn_ind !! crown area of individuals
136	!! @tex ($m^{2}$) @endtex
137	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
138	!! pool after the annual release for each
139	!! compartment (10 + 1 : input from year of land
140	!! cover change) (npts,0:10,nwp)
141	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
142	!! pool after the annual release for each
143	!! compartment (100 + 1 : input from year of land
144	!! cover change) (npts,0:100,nwp)
145	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: flux10 !! annual release from the 10/100 year-turnover
146	!! pool compartments (npts,10,nwp)
147	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: flux100 !! annual release from the 100/100 year-turnover
148	!! pool compartments (npts,100,nwp)
149	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class
150	!! (unitless)
151	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
152	!! @tex ($gC m^{-2}$) @endtex
153	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm !! "long term" net primary productivity
154	!! @tex ($gC m^{-2} year^{-1}$) @endtex
155	REAL(r_std),DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout):: litter !! metabolic and structural litter, above and
156	!! below ground @tex ($gC m^{-2}$) @endtex
157	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(inout):: litter_avail
158	REAL(r_std), DIMENSION(npts,nlitt,nvm) , INTENT(inout):: litter_not_avail
159	REAL(r_std),DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
160	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
161	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
162	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
163	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_1hr
164	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_10hr
165	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_100hr
166	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_1000hr
167	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements) :: fuel_all_type
168	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements,4) :: fuel_type_frac
169
170	!! 0.4 Local variables
171
172	INTEGER(i_std) :: i, j, k, l, m !! indices (unitless)
173	REAL(r_std),DIMENSION(npts,nelements) :: bm_new !! biomass increase @tex ($gC m^{-2}$) @endtex
174	REAL(r_std),DIMENSION(npts,nparts,nelements) :: biomass_loss !! biomass loss @tex ($gC m^{-2}$) @endtex
175	REAL(r_std) :: above !! aboveground biomass @tex ($gC m^{-2}$) @endtex
176	REAL(r_std),DIMENSION(npts,nlitt,nlevs,nelements) :: dilu_lit !! Litter dilution @tex ($gC m^{-2}$) @endtex
177	REAL(r_std),DIMENSION(npts,ncarb) :: dilu_soil_carbon !! Soil Carbondilution @tex ($gC m^{-2}$) @endtex
178	REAL(r_std),DIMENSION(npts,ndeep,ncarb) :: dilu_soil_carbon_vertres !!vertically-resolved Soil Carbondilution (gC/mÂ²)
179
180	REAL(r_std),DIMENSION(nvm) :: delta_veg !! changes in "maximal" coverage fraction of PFT
181	REAL(r_std) :: delta_veg_sum !! sum of delta_veg
182	REAL(r_std),DIMENSION(npts,nvm) :: delta_ind !! change in number of individuals
183
184	!_ ================================================================================================================================
185
186	IF (printlev>=3) WRITE(numout,*) 'Entering lcchange_main'
187
188	!! 1. initialization
189
190	prod10(:,0,:) = zero
191	prod100(:,0,:) = zero
192	above = zero
193	convflux(:,:) = zero
194	cflux_prod10(:,:) = zero
195	cflux_prod100(:,:) = zero
196	delta_ind(:,:) = zero
197	delta_veg(:) = zero
198	dilu_soil_carbon_vertres(:,:,:) = zero
199
200	!! 3. calculation of changes in carbon stocks and biomass by land cover change\n
201
202	DO i = 1, npts ! Loop over # pixels - domain size
203
204	!! 3.1 initialization of carbon stocks\n
205	delta_veg(:) = veget_cov_max_new(i,:)-veget_cov_max_old(i,:)
206	delta_veg_sum = SUM(delta_veg,MASK=delta_veg.LT.0.)
207
208	dilu_lit(i,:,:,:) = zero
209	dilu_soil_carbon(i,:) = zero
210	biomass_loss(i,:,:) = zero
211
212	!! 3.2 if vegetation coverage decreases, compute dilution of litter, soil carbon, and biomass.\n
213	DO j=2, nvm
214	IF ( delta_veg(j) < -min_stomate ) THEN
215	dilu_lit(i,:,:,:) = dilu_lit(i,:,:,:) + delta_veg(j)*litter(i,:,j,:,:) / delta_veg_sum
216	biomass_loss(i,:,:) = biomass_loss(i,:,:) + biomass(i,j,:,:)*delta_veg(j) / delta_veg_sum
217	IF ( ok_pc ) THEN
218	dilu_soil_carbon_vertres(i,:,iactive)=dilu_soil_carbon_vertres(i,:,iactive) + &
219	delta_veg(j) * deepC_a(i,:,j) / delta_veg_sum
220	dilu_soil_carbon_vertres(i,:,islow)=dilu_soil_carbon_vertres(i,:,islow) + &
221	delta_veg(j) * deepC_s(i,:,j) / delta_veg_sum
222	dilu_soil_carbon_vertres(i,:,ipassive)=dilu_soil_carbon_vertres(i,:,ipassive) + &
223	delta_veg(j) * deepC_p(i,:,j) / delta_veg_sum
224	ELSE
225	dilu_soil_carbon(i,:) = dilu_soil_carbon(i,:) + delta_veg(j) * carbon(i,:,j) / delta_veg_sum
226	ENDIF
227	ENDIF
228	ENDDO
229
230	!! 3.3
231	DO j=2, nvm ! Loop over # PFTs
232
233	!! 3.3.1 The case that vegetation coverage of PFTj increases
234	IF ( delta_veg(j) > min_stomate) THEN
235
236	!! 3.3.1.1 Initial setting of new establishment
237	IF (veget_cov_max_old(i,j) .LT. min_stomate) THEN
238	IF (is_tree(j)) THEN
239
240	! cn_sapl(j)=0.5; stomate_data.f90
241	cn_ind(i,j) = cn_sapl(j)
242	ELSE
243	cn_ind(i,j) = un
244	ENDIF
245	ind(i,j)= delta_veg(j) / cn_ind(i,j)
246	PFTpresent(i,j) = .TRUE.
247	everywhere(i,j) = 1.
248	senescence(i,j) = .FALSE.
249	age(i,j) = zero
250
251	! large_value = 1.E33_r_std
252	when_growthinit(i,j) = large_value
253	leaf_frac(i,j,1) = 1.0
254	npp_longterm(i,j) = npp_longterm_init
255	lm_lastyearmax(i,j) = bm_sapl(j,ileaf,icarbon) * ind(i,j)
256	ENDIF
257	IF ( cn_ind(i,j) > min_stomate ) THEN
258	delta_ind(i,j) = delta_veg(j) / cn_ind(i,j)
259	ENDIF
260
261	!! 3.3.1.2 Update of biomass in each each carbon stock component
262	!! Update of biomass in each each carbon stock component (leaf, sapabove, sapbelow,
263	!> heartabove, heartbelow, root, fruit, and carbres)\n
264	DO k = 1, nparts ! loop over # carbon stock components, nparts = 8; stomate_constant.f90
265	DO l = 1,nelements ! loop over # elements
266
267	bm_new(i,l) = delta_ind(i,j) * bm_sapl(j,k,l)
268	IF (veget_cov_max_old(i,j) .GT. min_stomate) THEN
269
270	! in the case that bm_new is overestimated compared with biomass?
271	IF ((bm_new(i,l)/delta_veg(j)) > biomass(i,j,k,l)) THEN
272	bm_new(i,l) = biomass(i,j,k,l)*delta_veg(j)
273	ENDIF
274	ENDIF
275	biomass(i,j,k,l) = ( biomass(i,j,k,l) * veget_cov_max_old(i,j) + bm_new(i,l) ) / veget_cov_max_new(i,j)
276	co2_to_bm(i,j) = co2_to_bm(i,j) + (bm_new(i,icarbon)* dt_days) / (one_year * veget_cov_max_new(i,j))
277	END DO ! loop over # elements
278	ENDDO ! loop over # carbon stock components
279
280	!! 3.3.1.3 Calculation of dilution in litter, soil carbon, and input of litter
281	!! In this 'IF statement', dilu_* is zero. Formulas for litter and soil carbon
282	!! could be shortend?? Are the following formulas correct?
283
284	! Litter
285	litter(i,:,j,:,:)=(litter(i,:,j,:,:) * veget_cov_max_old(i,j) + &
286	dilu_lit(i,:,:,:) * delta_veg(j)) / veget_cov_max_new(i,j)
287	!gmjc available and not available litter for grazing
288	! only not available litter increase/decrease, available litter will not
289	! change, due to tree litter can not be eaten
290	IF (is_grassland_manag(j) .AND. is_grassland_grazed(j)) THEN
291	litter_avail(i,:,j) = litter_avail(i,:,j) * veget_cov_max_old(i,j) / veget_cov_max_new(i,j)
292	litter_not_avail(i,:,j) = litter(i,:,j,iabove,icarbon) - litter_avail(i,:,j)
293	ENDIF
294	!end gmjc
295	IF ( ok_pc ) THEN
296	deepC_a(i,:,j)=(deepC_a(i,:,j) * veget_cov_max_old(i,j) + &
297	dilu_soil_carbon_vertres(i,:,iactive) * delta_veg(j)) / veget_cov_max_new(i,j)
298	deepC_s(i,:,j)=(deepC_s(i,:,j) * veget_cov_max_old(i,j) + &
299	dilu_soil_carbon_vertres(i,:,islow) * delta_veg(j)) / veget_cov_max_new(i,j)
300	deepC_p(i,:,j)=(deepC_p(i,:,j) * veget_cov_max_old(i,j) + &
301	dilu_soil_carbon_vertres(i,:,ipassive) * delta_veg(j)) / veget_cov_max_new(i,j)
302	ELSE
303	! Soil carbon
304	carbon(i,:,j)=(carbon(i,:,j) * veget_cov_max_old(i,j) + dilu_soil_carbon(i,:) * delta_veg(j)) / veget_cov_max_new(i,j)
305	ENDIF
306
307	DO l = 1,nelements
308
309	! Litter input
310	bm_to_litter(i,j,isapbelow,l) = bm_to_litter(i,j,isapbelow,l) + &
311	& biomass_loss(i,isapbelow,l)*delta_veg(j) / veget_cov_max_new(i,j)
312	bm_to_litter(i,j,iheartbelow,l) = bm_to_litter(i,j,iheartbelow,l) + biomass_loss(i,iheartbelow,l) *delta_veg(j) &
313	& / veget_cov_max_new(i,j)
314	bm_to_litter(i,j,iroot,l) = bm_to_litter(i,j,iroot,l) + biomass_loss(i,iroot,l)*delta_veg(j) / veget_cov_max_new(i,j)
315	bm_to_litter(i,j,ifruit,l) = bm_to_litter(i,j,ifruit,l) + biomass_loss(i,ifruit,l)*delta_veg(j) / veget_cov_max_new(i,j)
316	bm_to_litter(i,j,icarbres,l) = bm_to_litter(i,j,icarbres,l) + &
317	& biomass_loss(i,icarbres,l) *delta_veg(j) / veget_cov_max_new(i,j)
318	bm_to_litter(i,j,ileaf,l) = bm_to_litter(i,j,ileaf,l) + biomass_loss(i,ileaf,l)*delta_veg(j) / veget_cov_max_new(i,j)
319
320	END DO
321
322	age(i,j)=age(i,j)*veget_cov_max_old(i,j)/veget_cov_max_new(i,j)
323
324	!! 3.3.2 The case that vegetation coverage of PFTj is no change or decreases
325	ELSE
326
327	!! 3.3.2.1 Biomass export
328	! coeff_lcchange_*: Coeff of biomass export for the year, decade, and century
329	above = biomass(i,j,isapabove,icarbon) + biomass(i,j,iheartabove,icarbon)
330	convflux(i,iwplcc) = convflux(i,iwplcc) - ( coeff_lcchange_1(j) * above * delta_veg(j) )
331	prod10(i,0,iwplcc) = prod10(i,0,iwplcc) - ( coeff_lcchange_10(j) * above * delta_veg(j) )
332	prod100(i,0,iwplcc) = prod100(i,0,iwplcc) - ( coeff_lcchange_100(j) * above * delta_veg(j) )
333
334	ENDIF ! End if PFT's coverage reduction
335
336	ENDDO ! Loop over # PFTs
337
338	!! 2.4 update 10 year-turnover pool content following flux emission
339	!! (linear decay (10%) of the initial carbon input)
340	DO l = 0, 8
341	m = 10 - l
342	cflux_prod10(i,iwplcc) = cflux_prod10(i,iwplcc) + flux10(i,m,iwplcc)
343	prod10(i,m,iwplcc) = prod10(i,m-1,iwplcc) - flux10(i,m-1,iwplcc)
344	flux10(i,m,iwplcc) = flux10(i,m-1,iwplcc)
345
346	IF (prod10(i,m,iwplcc) .LT. 1.0) prod10(i,m,iwplcc) = zero
347	ENDDO
348
349	cflux_prod10(i,iwplcc) = cflux_prod10(i,iwplcc) + flux10(i,1,iwplcc)
350	flux10(i,1,iwplcc) = 0.1 * prod10(i,0,iwplcc)
351	prod10(i,1,iwplcc) = prod10(i,0,iwplcc)
352
353	!! 2.5 update 100 year-turnover pool content following flux emission\n
354	DO l = 0, 98
355	m = 100 - l
356	cflux_prod100(i,iwplcc) = cflux_prod100(i,iwplcc) + flux100(i,m,iwplcc)
357	prod100(i,m,iwplcc) = prod100(i,m-1,iwplcc) - flux100(i,m-1,iwplcc)
358	flux100(i,m,iwplcc) = flux100(i,m-1,iwplcc)
359
360	IF (prod100(i,m,iwplcc).LT.1.0) prod100(i,m,iwplcc) = zero
361	ENDDO
362
363	cflux_prod100(i,iwplcc) = cflux_prod100(i,iwplcc) + flux100(i,1,iwplcc)
364	flux100(i,1,iwplcc) = 0.01 * prod100(i,0,iwplcc)
365	prod100(i,1,iwplcc) = prod100(i,0,iwplcc)
366	prod10(i,0,iwplcc) = zero
367	prod100(i,0,iwplcc) = zero
368
369
370
371	ENDDO ! Loop over # pixels - domain size
372
373	!! We redistribute the updated litter into four fuel classes, so that
374	!! the balance between aboveground litter and fuel is mainted. The subtraction
375	!! of fuel burned by land cover change fires from the fuel pool is made here.
376	fuel_all_type(:,:,:,:) = fuel_1hr(:,:,:,:) + fuel_10hr(:,:,:,:) + &
377	fuel_100hr(:,:,:,:) + fuel_1000hr(:,:,:,:)
378	fuel_type_frac(:,:,:,:,:) = 0.25
379	WHERE(fuel_all_type(:,:,:,:) > min_stomate)
380	fuel_type_frac(:,:,:,:,1) = fuel_1hr(:,:,:,:)/fuel_all_type(:,:,:,:)
381	fuel_type_frac(:,:,:,:,2) = fuel_10hr(:,:,:,:)/fuel_all_type(:,:,:,:)
382	fuel_type_frac(:,:,:,:,3) = fuel_100hr(:,:,:,:)/fuel_all_type(:,:,:,:)
383	fuel_type_frac(:,:,:,:,4) = fuel_1000hr(:,:,:,:)/fuel_all_type(:,:,:,:)
384	ENDWHERE
385	DO j=1,nvm
386	fuel_1hr(:,j,:,:) = litter(:,:,j,iabove,:) * fuel_type_frac(:,j,:,:,1)
387	fuel_10hr(:,j,:,:) = litter(:,:,j,iabove,:) * fuel_type_frac(:,j,:,:,2)
388	fuel_100hr(:,j,:,:) = litter(:,:,j,iabove,:) * fuel_type_frac(:,j,:,:,3)
389	fuel_1000hr(:,j,:,:) = litter(:,:,j,iabove,:) * fuel_type_frac(:,j,:,:,4)
390	END DO
391
392	!! 3. history
393	convflux = convflux/one_year*dt_days
394	cflux_prod10 = cflux_prod10/one_year*dt_days
395	cflux_prod100 = cflux_prod100/one_year*dt_days
396
397	IF (printlev>=4) WRITE(numout,*) 'Leaving lcchange_main'
398
399	END SUBROUTINE lcchange_main
400
401
402	!! The lcchange modelling including consideration of deforestation fires
403	SUBROUTINE lcchange_deffire ( npts, dt_days, veget_cov_max, veget_cov_max_new,&
404	biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
405	co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,flux10,flux100, &
406	prod10,prod100,&
407	convflux,&
408	cflux_prod10,cflux_prod100, leaf_frac,&
409	npp_longterm, lm_lastyearmax, litter, litter_avail, litter_not_avail, &
410	carbon,&
411	deepC_a, deepC_s, deepC_p,&
412	fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr,&
413	lcc,bafrac_deforest_accu,emideforest_litter_accu,emideforest_biomass_accu,&
414	deflitsup_total,defbiosup_total)
415
416	IMPLICIT NONE
417
418	!! 0. Variable and parameter declaration
419
420	!! 0.1 Input variables
421
422	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
423	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
424	!! (days)
425	REAL(r_std), DIMENSION(nvm, nparts,nelements), INTENT(in) :: bm_sapl !! biomass of sapling
426	!! @tex ($gC individual^{-1}$) @endtex
427	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: bafrac_deforest_accu !!cumulative deforestation fire burned fraction, unitless
428	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: emideforest_litter_accu !!cumulative deforestation fire carbon emissions from litter
429	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: emideforest_biomass_accu !!cumulative deforestation fire carbon emissions from tree biomass
430	REAL(r_std), DIMENSION(npts,nvm),INTENT(in) :: lcc !! land cover change happened at this day
431
432	!! 0.2 Output variables
433
434	REAL(r_std), DIMENSION(:,:), INTENT(out) :: convflux !! release during first year following land cover
435	!! change (npts,nwp)
436	REAL(r_std), DIMENSION(:,:), INTENT(out) :: cflux_prod10 !! total annual release from the 10 year-turnover
437	!! pool @tex ($gC m^{-2}$) @endtex
438	!! (npts,nwp)
439	REAL(r_std), DIMENSION(:,:), INTENT(out) :: cflux_prod100 !! total annual release from the 100 year-
440	!! turnover pool @tex ($gC m^{-2}$) @endtex
441	!! (npts,nwp)
442	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: turnover_daily !! Turnover rates
443	!! @tex ($gC m^{-2} day^{-1}$) @endtex
444
445	!! 0.3 Modified variables
446
447	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_cov_max !! "maximal" coverage fraction of a PFT (LAI ->
448	!! infinity) on ground (unitless)
449	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_cov_max_new !! new "maximal" coverage fraction of a PFT (LAI
450	!! -> infinity) on ground
451	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: biomass !! biomass @tex ($gC m^{-2}$) @endtex
452	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Number of individuals @tex ($m^{-2}$) @endtex
453	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! mean age (years)
454	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: senescence !! plant senescent (only for deciduous trees) Set
455	!! to .FALSE. if PFT is introduced or killed
456	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! Is pft there (unitless)
457	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or very
458	!! localized (unitless)
459	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit !! how many days ago was the beginning of the
460	!! growing season (days)
461	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm !! biomass uptaken
462	!! @tex ($gC m^{-2} day^{-1}$) @endtex
463	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! conversion of biomass to litter
464	!! @tex ($gC m^{-2} day^{-1}$) @endtex
465	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: cn_ind !! crown area of individuals
466	!! @tex ($m^{2}$) @endtex
467	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
468	!! pool after the annual release for each
469	!! compartment (10 + 1 : input from year of land
470	!! cover change) (npts,0:10,nwp)
471	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
472	!! pool after the annual release for each
473	!! compartment (100 + 1 : input from year of land
474	!! cover change) (npts,0:100,nwp)
475	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: flux10 !! annual release from the 10/100 year-turnover
476	!! pool compartments (npts,0:10,nwp)
477	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: flux100 !! annual release from the 10/100 year-turnover
478	!! pool compartments (npts,0:100,nwp)
479	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class
480	!! (unitless)
481	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
482	!! @tex ($gC m^{-2}$) @endtex
483	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm !! "long term" net primary productivity
484	!! @tex ($gC m^{-2} year^{-1}$) @endtex
485	REAL(r_std),DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout):: litter !! metabolic and structural litter, above and
486	!! below ground @tex ($gC m^{-2}$) @endtex
487	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(inout):: litter_avail
488	REAL(r_std), DIMENSION(npts,nlitt,nvm) , INTENT(inout):: litter_not_avail
489	REAL(r_std),DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
490	!! @tex ($gC m^{-2}$) @endtex
491	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
492	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
493	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
494
495	REAL(r_std),DIMENSION(npts,nvm), INTENT(inout) :: deflitsup_total
496	REAL(r_std),DIMENSION(npts,nvm), INTENT(inout) :: defbiosup_total
497	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_1hr
498	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_10hr
499	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_100hr
500	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements),INTENT(inout) :: fuel_1000hr
501
502	!! 0.4 Local variables
503
504	INTEGER(i_std) :: i, j, k, l, m, ilit, ipart !! indices (unitless)
505	REAL(r_std),DIMENSION(npts,nelements) :: bm_new !! biomass increase @tex ($gC m^{-2}$) @endtex
506	REAL(r_std),DIMENSION(npts,nparts,nelements) :: biomass_loss !! biomass loss @tex ($gC m^{-2}$) @endtex
507	REAL(r_std) :: above !! aboveground biomass @tex ($gC m^{-2}$) @endtex
508	REAL(r_std),DIMENSION(npts,nlitt,nlevs,nelements) :: dilu_lit !! Litter dilution @tex ($gC m^{-2}$) @endtex
509	REAL(r_std),DIMENSION(npts,ncarb) :: dilu_soil_carbon !! Soil Carbondilution @tex ($gC m^{-2}$) @endtex
510	REAL(r_std),DIMENSION(npts,ndeep,ncarb) :: dilu_soil_carbon_vertres !!vertically-resolved Soil Carbondilution (gC/mÂ²)
511
512	REAL(r_std),DIMENSION(nvm) :: delta_veg !! changes in "maximal" coverage fraction of PFT
513	REAL(r_std) :: delta_veg_sum !! sum of delta_veg
514	REAL(r_std),DIMENSION(npts,nvm) :: delta_ind !! change in number of individuals
515	REAL(r_std),DIMENSION(npts,nvm,nlitt) :: deforest_litter_surplus !! Surplus in ground litter for deforested land after
516	!! accounting for fire emissions
517	REAL(r_std),DIMENSION(npts,nvm,nparts) :: deforest_biomass_surplus !!Surplus in live biomass for deforested forest
518	!!after accounting for fire emissions
519	REAL(r_std),DIMENSION(npts,nvm,nlitt) :: deforest_litter_deficit
520	REAL(r_std),DIMENSION(npts,nvm,nparts) :: deforest_biomass_deficit
521	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements) :: fuel_all_type
522	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements,4) :: fuel_type_frac
523
524	REAL(r_std),DIMENSION(npts,nvm) :: pool_start_pft !! change in number of individuals
525	REAL(r_std),DIMENSION(npts) :: pool_start !! change in number of individuals
526	REAL(r_std),DIMENSION(npts,nvm) :: pool_end_pft !! change in number of individuals
527	REAL(r_std),DIMENSION(npts) :: pool_end !! change in number of individuals
528	REAL(r_std),DIMENSION(npts) :: outflux !! change in number of individuals
529	!_ ================================================================================================================================
530
531
532
533	pool_start_pft(:,:) = SUM(biomass(:,:,:,icarbon),DIM=3) &
534	+ SUM(SUM(litter(:,:,:,:,icarbon),DIM=2),DIM=3) &
535	+ SUM(carbon(:,:,:),DIM=2) &
536	+ SUM(bm_to_litter(:,:,:,icarbon),DIM=3) &
537	+ SUM(turnover_daily(:,:,:,icarbon),DIM=3)
538
539	pool_start(:) = SUM(pool_start_pft(:,:)*veget_cov_max(:,:),DIM=2) &
540	+ SUM(prod10(:,:,iwplcc),DIM=2) + SUM(prod100(:,:,iwplcc),DIM=2)
541
542
543	deforest_biomass_surplus(:,:,:) = zero
544	deforest_litter_surplus(:,:,:) = zero
545	deforest_biomass_deficit(:,:,:) = zero
546	deforest_litter_deficit(:,:,:) = zero
547
548	IF (printlev>=3) WRITE(numout,*) 'Entering lcchange_main'
549
550	!! 1. initialization
551
552	prod10(:,0,:) = zero
553	prod100(:,0,:) = zero
554	above = zero
555	convflux(:,:) = zero
556	cflux_prod10(:,:) = zero
557	cflux_prod100(:,:) = zero
558	delta_ind(:,:) = zero
559	delta_veg(:) = zero
560	dilu_soil_carbon_vertres(:,:,:) =zero
561	!! 2. calculation of changes in carbon stocks and biomass by land cover change\n
562
563	DO i = 1, npts ! Loop over # pixels - domain size
564
565	!! 2.1 initialization of carbon stocks\n
566	delta_veg(:) = veget_cov_max_new(i,:)-veget_cov_max(i,:)
567	delta_veg_sum = SUM(delta_veg,MASK=delta_veg.LT.0.) !note `delta_veg_sum` is a negative number
568
569	dilu_lit(i,:,:,:) = zero
570	dilu_soil_carbon(i,:) = zero
571	biomass_loss(i,:,:) = zero
572
573	!! 2.2 Compute dilution pool of litter, soil carbon, and biomass for
574	!! decreasing PFTs.
575	DO j=2, nvm
576	IF ( delta_veg(j) < -min_stomate ) THEN
577
578	! We make distinction between tree and grass because tree cover reduction might be due to fires.
579	! The litter that is burned in fire should be excluded from diluting litter pool.
580	IF (is_tree(j)) THEN
581	deforest_litter_surplus(i,j,:) = -1delta_veg(j)litter(i,:,j,iabove,icarbon) - emideforest_litter_accu(i,j,:,icarbon)
582
583	! Here we compensate the litter burned by deforestation fire if it's higher than the litter available for
584	! burning. It follows the same logic as biomass which is described below.
585	DO ilit = 1,nlitt
586	IF (deforest_litter_surplus(i,j,ilit) < zero) THEN
587	IF (veget_cov_max_new(i,j) < min_stomate) THEN
588	!WRITE (numout,*) 'Cumulative deforestation fire emission exceeds litter for point',i,',PFT ',j, &
589	! 'However the new veget_cov_max is zero, there is not remaining litter to be diluted'
590	!STOP
591	deforest_litter_deficit(i,j,ilit) = deforest_litter_surplus(i,j,ilit)
592
593	ELSE IF (litter(i,ilit,j,iabove,icarbon)*veget_cov_max_new(i,j) < -deforest_litter_surplus(i,j,ilit)) THEN
594	!WRITE (numout,*) 'Cumulative deforestation fire emission exceeds litter for point',i,',PFT ',j, &
595	! 'However the remaing litter is not engough for diluting'
596	!STOP
597	deforest_litter_deficit(i,j,ilit) = deforest_litter_surplus(i,j,ilit)
598	ELSE
599	litter(i,ilit,j,iabove,icarbon) = ( litter(i,ilit,j,iabove,icarbon)*veget_cov_max_new(i,j) &
600	+ deforest_litter_surplus(i,j,ilit) )/veget_cov_max_new(i,j)
601	END IF
602	ELSE
603	dilu_lit(i,ilit,iabove,icarbon) = dilu_lit(i,ilit,iabove,icarbon) -1 * deforest_litter_surplus(i,j,ilit)
604	END IF
605	END DO
606	dilu_lit(i,:,ibelow,:) = dilu_lit(i,:,ibelow,:) + delta_veg(j)*litter(i,:,j,ibelow,:)
607	ELSE
608	dilu_lit(i,:,:,:) = dilu_lit(i,:,:,:) + delta_veg(j)*litter(i,:,j,:,:)
609	END IF
610
611	IF (is_tree(j)) THEN
612	deforest_biomass_surplus(i,j,:) = -1delta_veg(j)biomass(i,j,:,icarbon) - emideforest_biomass_accu(i,j,:,icarbon)
613	! Here we check if the biomass burned by deforestation fires is higher than the amount
614	! that could be deforested, if yes, the extra burned biomass is compensated by the biomass
615	! that is not deforested. Here we assume that if this happens for one deforested tree PFT,
616	! it happens for all deforested tree PFTs, so that we don't assume this extra burned biomass
617	! could be compenstated by other tree PFTs.
618	DO ipart = 1,nparts
619	IF (deforest_biomass_surplus(i,j,ipart) < zero) THEN
620	IF (veget_cov_max_new(i,j) < min_stomate) THEN
621	!WRITE (numout,*) 'Cumulative deforestation fire emission exceeds biomass for point',i,',PFT ',j, &
622	! 'However the new veget_cov_max is zero, there is not remaining biomass to be diluted'
623	!STOP
624	deforest_biomass_deficit(i,j,ipart) = deforest_biomass_surplus(i,j,ipart)
625
626	ELSE IF (biomass(i,j,ipart,icarbon)*veget_cov_max_new(i,j) < -deforest_biomass_surplus(i,j,ipart)) THEN
627	!WRITE (numout,*) 'Cumulative deforestation fire emission exceeds biomass for point',i,',PFT ',j, &
628	! 'However the remaing biomass is not engough for diluting'
629	!STOP
630	deforest_biomass_deficit(i,j,ipart) = deforest_biomass_surplus(i,j,ipart)
631	ELSE
632	biomass(i,j,ipart,icarbon) = ( biomass(i,j,ipart,icarbon)*veget_cov_max_new(i,j) &
633	+ deforest_biomass_surplus(i,j,ipart) )/veget_cov_max_new(i,j)
634	END IF
635	ELSE
636	biomass_loss(i,ipart,icarbon) = biomass_loss(i,ipart,icarbon) -1 * deforest_biomass_surplus(i,j,ipart)
637	END IF
638	END DO
639	ELSE
640	biomass_loss(i,:,:) = biomass_loss(i,:,:) + biomass(i,j,:,:)*delta_veg(j)
641	END IF
642
643	!IF (ANY( deforest_biomass_surplus(i,j,:) .LT. 0.0 ) .OR. ANY( deforest_litter_surplus(i,j,:) .LT. 0.0 ) ) THEN
644	! STOP 'Negative biomass or litter surplus'
645	!ENDIF
646
647	IF ( ok_pc ) THEN
648	dilu_soil_carbon_vertres(i,:,iactive)=dilu_soil_carbon_vertres(i,:,iactive) + &
649	delta_veg(j) * deepC_a(i,:,j) / delta_veg_sum
650	dilu_soil_carbon_vertres(i,:,islow)=dilu_soil_carbon_vertres(i,:,islow) + &
651	delta_veg(j) * deepC_s(i,:,j) / delta_veg_sum
652	dilu_soil_carbon_vertres(i,:,ipassive)=dilu_soil_carbon_vertres(i,:,ipassive) + &
653	delta_veg(j) * deepC_p(i,:,j) / delta_veg_sum
654	ELSE
655	dilu_soil_carbon(i,:) = dilu_soil_carbon(i,:) + delta_veg(j) * carbon(i,:,j) / delta_veg_sum
656	ENDIF
657	ENDIF
658	ENDDO !nbpts
659
660
661	! Note here `biomass_loss` and `dilu_lit` will change their sign from negative to positive
662	IF ( delta_veg_sum < -min_stomate ) THEN
663	biomass_loss(i,:,:) = biomass_loss(i,:,:) / delta_veg_sum
664	dilu_lit(i,:,:,:) = dilu_lit(i,:,:,:) / delta_veg_sum
665	END IF
666
667
668	!! 2.3 Dilut the litter, soil carbon from decreasing PFTs to increasing ones.
669	!! Establish new biomass for increasing PFTs.
670	DO j=2, nvm ! Loop over # PFTs
671
672	!! 2.3.1 The case that vegetation coverage of PFTj increases
673	IF ( delta_veg(j) > min_stomate) THEN
674
675	!! 2.3.1.1 The PFTj increased from zero to non-zeor, we have to
676	!! initialize it by setting new establishment
677	IF (veget_cov_max(i,j) .LT. min_stomate) THEN
678	IF (is_tree(j)) THEN
679	cn_ind(i,j) = cn_sapl(j) ! cn_sapl(j)=0.5; stomate_data.f90
680	ELSE
681	cn_ind(i,j) = un
682	ENDIF
683
684	ind(i,j)= delta_veg(j) / cn_ind(i,j)
685	PFTpresent(i,j) = .TRUE.
686	everywhere(i,j) = 1.
687	senescence(i,j) = .FALSE.
688	age(i,j) = zero
689	when_growthinit(i,j) = large_value ! large_value = 1.E33_r_std
690	leaf_frac(i,j,1) = 1.0
691	npp_longterm(i,j) = npp_longterm_init
692	lm_lastyearmax(i,j) = bm_sapl(j,ileaf,icarbon) * ind(i,j)
693	ENDIF
694
695
696	! Calculate individual density increase because of coverage increase
697	IF ( cn_ind(i,j) > min_stomate ) THEN
698	delta_ind(i,j) = delta_veg(j) / cn_ind(i,j)
699	ENDIF
700	!! 2.3.1.2 The increase in `ind` should be companied by increase in
701	!! biomass, we do this by assuming increased `ind` are saplings.
702	DO k = 1, nparts ! loop over # carbon stock components, nparts = 8; stomate_constant.f90
703	DO l = 1,nelements ! loop over # elements
704	bm_new(i,l) = delta_ind(i,j) * bm_sapl(j,k,l)
705	IF (veget_cov_max(i,j) .GT. min_stomate) THEN
706	! Adjust bm_new equal to existing biomass if it's
707	! larger than the latter
708	IF ((bm_new(i,l)/delta_veg(j)) > biomass(i,j,k,l)) THEN
709	bm_new(i,l) = biomass(i,j,k,l)*delta_veg(j)
710	ENDIF
711	ENDIF
712	biomass(i,j,k,l) = ( biomass(i,j,k,l) * veget_cov_max(i,j) + bm_new(i,l) ) / veget_cov_max_new(i,j)
713	co2_to_bm(i,j) = co2_to_bm(i,j) + (bm_new(i,icarbon)* dt_days) / (one_year * veget_cov_max_new(i,j))
714	END DO ! loop over # elements
715	ENDDO ! loop over # carbon stock components
716
717	!! 2.3.1.3 Tow tasks are done here:
718	!! A. We transfer the litter and soil carbon from the
719	!! reduced PFTs to the increases PFTs.
720
721	! Litter
722	litter(i,:,j,:,:)=(litter(i,:,j,:,:) * veget_cov_max(i,j) + &
723	dilu_lit(i,:,:,:) * delta_veg(j)) / veget_cov_max_new(i,j)
724
725	!!######################This part needs to be discussed with JinFeng ############
726	!gmjc available and not available litter for grazing
727	! only not available litter increase/decrease, available litter will not
728	! change, due to tree litter can not be eaten
729	IF (is_grassland_manag(j) .AND. is_grassland_grazed(j)) THEN
730	litter_avail(i,:,j) = litter_avail(i,:,j) * veget_cov_max(i,j) / veget_cov_max_new(i,j)
731	litter_not_avail(i,:,j) = litter(i,:,j,iabove,icarbon) - litter_avail(i,:,j)
732	ENDIF
733	!end gmjc
734	!!###############################################################################
735
736	! Soil carbon
737	IF ( ok_pc ) THEN
738	deepC_a(i,:,j)=(deepC_a(i,:,j) * veget_cov_max(i,j) + &
739	dilu_soil_carbon_vertres(i,:,iactive) * delta_veg(j)) / veget_cov_max_new(i,j)
740	deepC_s(i,:,j)=(deepC_s(i,:,j) * veget_cov_max(i,j) + &
741	dilu_soil_carbon_vertres(i,:,islow) * delta_veg(j)) / veget_cov_max_new(i,j)
742	deepC_p(i,:,j)=(deepC_p(i,:,j) * veget_cov_max(i,j) + &
743	dilu_soil_carbon_vertres(i,:,ipassive) * delta_veg(j)) / veget_cov_max_new(i,j)
744	ELSE
745	carbon(i,:,j)=(carbon(i,:,j) * veget_cov_max(i,j) + dilu_soil_carbon(i,:) * delta_veg(j)) / veget_cov_max_new(i,j)
746	ENDIF
747
748	!! B. For the biomass pool of reducing PFTs, we cannot transfer them directly to the
749	!! increasing PFTs, because the latter ones are treated with new sapling estalishement
750	!! in section 2.3.1.2. So we assume the non-harvestable biomass of reducing PFTs will
751	!! go to litter pool via `bm_to_litter`, and these are further directly transferred to
752	!! the increasing PFTs.
753	!!
754	!! The non-harvestable parts are: isapbelow,iheartbelow,iroot,icarbres,ileaf,ifruit
755	!! Note that the icarbres,ileaf,ifruit could be burned in deforestation fires, the
756	!! emissions from these parts are already subtracted from `biomass_loss`, as done
757	!! in section 2.2. The harvestable biomass parts go to harvest pool and this will done
758	!! in the section for the reducing PFTs.
759	DO l = 1,nelements
760
761	bm_to_litter(i,j,isapbelow,l) = bm_to_litter(i,j,isapbelow,l) + &
762	& biomass_loss(i,isapbelow,l)*delta_veg(j) / veget_cov_max_new(i,j)
763	bm_to_litter(i,j,iheartbelow,l) = bm_to_litter(i,j,iheartbelow,l) + biomass_loss(i,iheartbelow,l) *delta_veg(j) &
764	& / veget_cov_max_new(i,j)
765	bm_to_litter(i,j,iroot,l) = bm_to_litter(i,j,iroot,l) + biomass_loss(i,iroot,l)*delta_veg(j) / veget_cov_max_new(i,j)
766	bm_to_litter(i,j,ifruit,l) = bm_to_litter(i,j,ifruit,l) + biomass_loss(i,ifruit,l)*delta_veg(j) / veget_cov_max_new(i,j)
767	bm_to_litter(i,j,icarbres,l) = bm_to_litter(i,j,icarbres,l) + &
768	& biomass_loss(i,icarbres,l) *delta_veg(j) / veget_cov_max_new(i,j)
769	bm_to_litter(i,j,ileaf,l) = bm_to_litter(i,j,ileaf,l) + biomass_loss(i,ileaf,l)*delta_veg(j) / veget_cov_max_new(i,j)
770	END DO
771
772	age(i,j)=age(i,j)*veget_cov_max(i,j)/veget_cov_max_new(i,j)
773
774	!! 2.3.2 The case that vegetation coverage of PFTj has no change or decreases.
775	ELSE
776
777	!! 2.3.2.1 Complete disappearing of PFTj, i.e., changes from non-zero
778	!! to zero.
779	IF ( veget_cov_max_new(i,j) .LT. min_stomate ) THEN
780	veget_cov_max_new(i,j)= zero
781	ind(i,j) = zero
782	biomass(i,j,:,:) = zero
783	PFTpresent(i,j) = .FALSE.
784	senescence(i,j) = .FALSE.
785	age(i,j) = zero
786	when_growthinit(i,j) = undef
787	everywhere(i,j) = zero
788	carbon(i,:,j) = zero
789	litter(i,:,j,:,:) = zero
790	litter_avail(i,:,j) = zero
791	litter_not_avail(i,:,j) = zero
792	bm_to_litter(i,j,:,:) = zero
793	turnover_daily(i,j,:,:) = zero
794	deepC_a(i,:,j) = zero
795	deepC_s(i,:,j) = zero
796	deepC_p(i,:,j) = zero
797	ENDIF
798	ENDIF ! The end the two cases: PFT-coverage reduction versus
799	! non-change-or-increase
800	ENDDO ! 2.3 Loop over # PFTs
801
802	!! 2.4 Biomass harvest and turnover of different harvest pools
803
804	!!?? Here we have some problem regarding grassland/cropland area dereasing,
805	!!?? Because their sapwood/heartwood aboveground are also treated as
806	!!?? wood products.
807
808	!! 2.4.1 We have already deforestation fire fluxes from sapwood/hearwood aboveground,
809	!! now we just assume the remaining unburned parts are harvested, as 10-year and
810	!! 100-year product pool.
811
812	print *,'delta_veg_sum',delta_veg_sum
813	print *,'prod10_in_lcc_before_assign',prod10(:,:,iwplcc)
814	print *,'biomass_loss',biomass_loss(:,:,:)
815	! Note before we divide biomass_loss by `delta_veg_sum` to convert it based on PFT area,
816	! Now we multiply it again by `delta_veg_sum` to convert it back based on grid cell area.
817	! Also note `delta_veg_sum` is negative, so we should multiply again by (-1)
818	above = (biomass_loss(i,isapabove,icarbon) + biomass_loss(i,iheartabove,icarbon))delta_veg_sum(-1)
819	convflux(i,iwplcc) = SUM(emideforest_biomass_accu(i,:,isapabove,icarbon)+emideforest_biomass_accu(i,:,iheartabove,icarbon))
820	prod10(i,0,iwplcc) = 0.4* above
821	prod100(i,0,iwplcc) = 0.6 * above
822	print *,'above_in_lcc_before_assign',above
823
824	!! 2.4.2 update 10 year-turnover pool content following flux emission
825	!! (linear decay (10%) of the initial carbon input)
826	DO l = 0, 8
827	m = 10 - l
828	cflux_prod10(i,iwplcc) = cflux_prod10(i,iwplcc) + flux10(i,m,iwplcc)
829	prod10(i,m,iwplcc) = prod10(i,m-1,iwplcc) - flux10(i,m-1,iwplcc)
830	flux10(i,m,iwplcc) = flux10(i,m-1,iwplcc)
831	IF (prod10(i,m,iwplcc) .LT. 1.0) prod10(i,m,iwplcc) = zero
832
833	ENDDO
834
835	cflux_prod10(i,iwplcc) = cflux_prod10(i,iwplcc) + flux10(i,1,iwplcc)
836	flux10(i,1,iwplcc) = 0.1 * prod10(i,0,iwplcc)
837	prod10(i,1,iwplcc) = prod10(i,0,iwplcc)
838
839	!! 3.5 update 100 year-turnover pool content following flux emission\n
840	DO l = 0, 98
841	m = 100 - l
842	cflux_prod100(i,iwplcc) = cflux_prod100(i,iwplcc) + flux100(i,m,iwplcc)
843	prod100(i,m,iwplcc) = prod100(i,m-1,iwplcc) - flux100(i,m-1,iwplcc)
844	flux100(i,m,iwplcc) = flux100(i,m-1,iwplcc)
845
846	IF (prod100(i,m,iwplcc).LT.1.0) prod100(i,m,iwplcc) = zero
847	ENDDO
848
849	cflux_prod100(i,iwplcc) = cflux_prod100(i,iwplcc) + flux100(i,1,iwplcc)
850	flux100(i,1,iwplcc) = 0.01 * prod100(i,0,iwplcc)
851	prod100(i,1,iwplcc) = prod100(i,0,iwplcc)
852	prod10(i,0,iwplcc) = zero
853	prod100(i,0,iwplcc) = zero
854
855	ENDDO ! Loop over # pixels - domain size
856	print *,'prod10_in_lcc_after_assign',prod10(:,:,iwplcc)
857
858	!!Jinfeng's grassland management module might should also be put here.
859
860	!! We redistribute the updated litter into four fuel classes, so that
861	!! the balance between aboveground litter and fuel is mainted. The subtraction
862	!! of fuel burned by land cover change fires from the fuel pool is made here.
863	fuel_all_type(:,:,:,:) = fuel_1hr(:,:,:,:) + fuel_10hr(:,:,:,:) + &
864	fuel_100hr(:,:,:,:) + fuel_1000hr(:,:,:,:)
865	fuel_type_frac(:,:,:,:,:) = 0.25
866	WHERE(fuel_all_type(:,:,:,:) > min_stomate)
867	fuel_type_frac(:,:,:,:,1) = fuel_1hr(:,:,:,:)/fuel_all_type(:,:,:,:)
868	fuel_type_frac(:,:,:,:,2) = fuel_10hr(:,:,:,:)/fuel_all_type(:,:,:,:)
869	fuel_type_frac(:,:,:,:,3) = fuel_100hr(:,:,:,:)/fuel_all_type(:,:,:,:)
870	fuel_type_frac(:,:,:,:,4) = fuel_1000hr(:,:,:,:)/fuel_all_type(:,:,:,:)
871	ENDWHERE
872	DO j=1,nvm
873	fuel_1hr(:,j,:,:) = litter(:,:,j,iabove,:) * fuel_type_frac(:,j,:,:,1)
874	fuel_10hr(:,j,:,:) = litter(:,:,j,iabove,:) * fuel_type_frac(:,j,:,:,2)
875	fuel_100hr(:,j,:,:) = litter(:,:,j,iabove,:) * fuel_type_frac(:,j,:,:,3)
876	fuel_1000hr(:,j,:,:) = litter(:,:,j,iabove,:) * fuel_type_frac(:,j,:,:,4)
877	END DO
878
879	!! 4. history
880
881	veget_cov_max(:,:) = veget_cov_max_new(:,:)
882	convflux(:,iwplcc) = convflux(:,iwplcc)/one_year*dt_days
883	cflux_prod10(:,iwplcc) = cflux_prod10(:,iwplcc)/one_year*dt_days
884	cflux_prod100(:,iwplcc) = cflux_prod100(:,iwplcc)/one_year*dt_days
885
886	pool_end_pft(:,:) = SUM(biomass(:,:,:,icarbon),DIM=3) &
887	+ SUM(SUM(litter(:,:,:,:,icarbon),DIM=2),DIM=3) &
888	+ SUM(carbon(:,:,:),DIM=2) &
889	+ SUM(bm_to_litter(:,:,:,icarbon),DIM=3) &
890	+ SUM(turnover_daily(:,:,:,icarbon),DIM=3)
891
892	pool_end(:) = SUM(pool_end_pft(:,:)*veget_cov_max(:,:),DIM=2) &
893	+ SUM(prod10(:,:,iwplcc),DIM=2) + SUM(prod100(:,:,iwplcc),DIM=2)
894
895
896	outflux(:) = SUM(SUM(emideforest_biomass_accu(:,:,:,icarbon),DIM=3),DIM=2) &
897	+ SUM(SUM(emideforest_litter_accu(:,:,:,icarbon),DIM=3),DIM=2) &
898	+ SUM(flux10(:,:,iwplcc),DIM=2) + SUM(flux100(:,:,iwplcc),DIM=2) &
899	- SUM(co2_to_bm(:,:)*veget_cov_max(:,:),DIM=2)
900
901	print *,"pool_start: ",pool_start(:)
902	print *,"pool_end: ",pool_end(:)
903	print *,"outflux: ",outflux(:)
904	print *,"pool_change: ",pool_start(:)-pool_end(:)
905	print *,'prod10_end_lcc',prod10(:,:,:)
906
907	deflitsup_total(:,:) = SUM(deforest_litter_surplus(:,:,:),dim=3)
908	defbiosup_total(:,:) = SUM(deforest_biomass_surplus(:,:,:),dim=3)
909
910	CALL histwrite (hist_id_stomate, 'dilu_lit_met', itime, &
911	dilu_lit(:,imetabolic,iabove,icarbon), npts, hori_index)
912	CALL histwrite (hist_id_stomate, 'dilu_lit_str', itime, &
913	dilu_lit(:,istructural,iabove,icarbon), npts, hori_index)
914
915
916	CALL histwrite (hist_id_stomate, 'SurpBioLEAF', itime, &
917	deforest_biomass_surplus(:,:,ileaf), npts*nvm, horipft_index)
918	CALL histwrite (hist_id_stomate, 'SurpBioRESERVE', itime, &
919	deforest_biomass_surplus(:,:,icarbres), npts*nvm, horipft_index)
920	CALL histwrite (hist_id_stomate, 'SurpBioFRUIT', itime, &
921	deforest_biomass_surplus(:,:,ifruit), npts*nvm, horipft_index)
922	CALL histwrite (hist_id_stomate, 'SurpBioSapABOVE', itime, &
923	deforest_biomass_surplus(:,:,isapabove), npts*nvm, horipft_index)
924	CALL histwrite (hist_id_stomate, 'SurpBioHeartABOVE', itime, &
925	deforest_biomass_surplus(:,:,iheartabove), npts*nvm, horipft_index)
926	CALL histwrite (hist_id_stomate, 'SurpBioSapBELOW', itime, &
927	deforest_biomass_surplus(:,:,isapbelow), npts*nvm, horipft_index)
928	CALL histwrite (hist_id_stomate, 'SurpBioHeartBELOW', itime, &
929	deforest_biomass_surplus(:,:,iheartbelow), npts*nvm, horipft_index)
930	CALL histwrite (hist_id_stomate, 'SurpBioROOT', itime, &
931	deforest_biomass_surplus(:,:,iroot), npts*nvm, horipft_index)
932	CALL histwrite (hist_id_stomate, 'SurpLitMET', itime, &
933	deforest_litter_surplus(:,:,imetabolic), npts*nvm, horipft_index)
934	CALL histwrite (hist_id_stomate, 'SurpLitSTR', itime, &
935	deforest_litter_surplus(:,:,istructural), npts*nvm, horipft_index)
936
937	CALL histwrite (hist_id_stomate, 'DefiBioLEAF', itime, &
938	deforest_biomass_deficit(:,:,ileaf), npts*nvm, horipft_index)
939	CALL histwrite (hist_id_stomate, 'DefiBioRESERVE', itime, &
940	deforest_biomass_deficit(:,:,icarbres), npts*nvm, horipft_index)
941	CALL histwrite (hist_id_stomate, 'DefiBioFRUIT', itime, &
942	deforest_biomass_deficit(:,:,ifruit), npts*nvm, horipft_index)
943	CALL histwrite (hist_id_stomate, 'DefiBioSapABOVE', itime, &
944	deforest_biomass_deficit(:,:,isapabove), npts*nvm, horipft_index)
945	CALL histwrite (hist_id_stomate, 'DefiBioHeartABOVE', itime, &
946	deforest_biomass_deficit(:,:,iheartabove), npts*nvm, horipft_index)
947	CALL histwrite (hist_id_stomate, 'DefiBioSapBELOW', itime, &
948	deforest_biomass_deficit(:,:,isapbelow), npts*nvm, horipft_index)
949	CALL histwrite (hist_id_stomate, 'DefiBioHeartBELOW', itime, &
950	deforest_biomass_deficit(:,:,iheartbelow), npts*nvm, horipft_index)
951	CALL histwrite (hist_id_stomate, 'DefiBioROOT', itime, &
952	deforest_biomass_deficit(:,:,iroot), npts*nvm, horipft_index)
953	CALL histwrite (hist_id_stomate, 'DefiLitMET', itime, &
954	deforest_litter_deficit(:,:,imetabolic), npts*nvm, horipft_index)
955	CALL histwrite (hist_id_stomate, 'DefiLitSTR', itime, &
956	deforest_litter_deficit(:,:,istructural), npts*nvm, horipft_index)
957
958
959	IF (printlev>=4) WRITE(numout,*) 'Leaving lcchange_main'
960
961	END SUBROUTINE lcchange_deffire
962
963
964	!SUBROUTINE lcc_neighbour_shift(ipts,neighbours,veget_cov_max,lcc,veget_cov_max_new)
965	! INTEGER(i_std), DIMENSION(npts,8), INTENT(in) :: neighbours !! indices of the 8 neighbours of each grid point
966	! !! (unitless);
967	! !! [1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW]
968
969
970	!END SUBROUTINE lcc_neighbour_shift
971
972	!print ,'end_biomass',SUM(SUM(biomass(:,:,:,icarbon),DIM=3)veget_cov_max(:,:),DIM=2)
973	!print ,'end_litter',SUM(SUM(SUM(litter(:,:,:,:,icarbon),DIM=2),DIM=3)veget_cov_max(:,:),DIM=2)
974	!print , 'end_soil',SUM(SUM(carbon(:,:,:),DIM=2)veget_cov_max(:,:),DIM=2)
975	!print ,'end_bm2lit',sum(SUM(bm_to_litter(:,:,:,icarbon),DIM=3)veget_cov_max(:,:),dim=2)
976	!print ,'end_turnover',sum(SUM(turnover_daily(:,:,:,icarbon),DIM=3)veget_cov_max(:,:),dim=2)
977	!print *,'end_prod10', SUM(prod10(:,:),DIM=2)
978	!print *,'end_prod100',SUM(prod100(:,:),DIM=2)
979
980	! !!block to check
981	! pool_end_pft(:,:) = SUM(biomass(:,:,:,icarbon),DIM=3) &
982	! + SUM(SUM(litter(:,:,:,:,icarbon),DIM=2),DIM=3) &
983	! + SUM(carbon(:,:,:),DIM=2) &
984	! + SUM(bm_to_litter(:,:,:,icarbon),DIM=3) &
985	! + SUM(turnover_daily(:,:,:,icarbon),DIM=3)
986	!
987	! pool_end(:) = SUM(pool_end_pft(:,:)*veget_cov_max(:,:),DIM=2) &
988	! + SUM(prod10(:,:),DIM=2) + SUM(prod100(:,:),DIM=2)
989	!
990	! outflux(:) = SUM(SUM(emideforest_biomass_accu(:,:,:,icarbon),DIM=3),DIM=2) &
991	! + SUM(SUM(emideforest_litter_accu(:,:,:,icarbon),DIM=3),DIM=2) &
992	! + SUM(flux10(:,:),DIM=2) + SUM(flux100,DIM=2) &
993	! - SUM(co2_to_bm(:,:)*veget_cov_max(:,:),DIM=2)
994	!
995	! print *,"pool_start: ",pool_start(:)
996	! print *,"pool_end: ",pool_end(:)
997	! print *,"outflux: ",outflux(:)
998	! print *,"pool_change: ",pool_start(:)-pool_end(:)
999	! !!end block to check
1000
1001
1002	END MODULE stomate_lcchange

Note: See TracBrowser for help on using the repository browser.

Download in other formats: