Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_stomate/stomate_lcchange.f90 @ 5491

Last change on this file since 5491 was 3771, checked in by albert.jornet, 8 years ago
Merge: new Grassland Management Module (GRM) from revision [3759/perso/jinfeng.chang/ORCHIDEE-MICT/ORCHIDEE/ORCHIDEE]. Done by Jinfeng.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 58.1 KB

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

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

Download in other formats: