Context Navigation

stomate_glcchange_SinAgeC_fh.f90 @ 7346

Last change on this file since 7346 was 3306, checked in by albert.jornet, 8 years ago

Merge: from perso branch revision [3305/perso/albert.jornet/ORCHIDEE-MIC-GLUC-GM31]
Fix: multiple fixes in LUC

New: introduction of GLUC done by Chao
New: GM 3.1 done by Jinfeng
Fix: we pass "veget_lastyear" in the position of "veget_max_new" when calling stomate_main, which is wrong because veget_max_new in stomate.f90 should mean the new veget_max.
Fix: the line "veget_max_tmp(:,:) = veget_max(:,:)" is always executed in the file stomate_lpj.f90, which makes their difficiation when calling either lcchange_deffire or lcchange_main useless.
Fix: multiple fixes in LUC

File size: 142.9 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_lcchange_fh
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 This module is a copy of stomate_lcchange. It includes the forestry
10	! harvest.
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): Including permafrost carbon
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/perso/albert.jornet/ORCHIDEE-MICT/src_stomate/stomate_lcchange.f90 $
20	!! $Date: 2015-07-30 15:38:45 +0200 (Thu, 30 Jul 2015) $
21	!! $Revision: 2847 $
22	!! \n
23	!_ ================================================================================================================================
24
25
26	MODULE stomate_glcchange_SinAgeC_fh
27
28	! modules used:
29
30	USE ioipsl_para
31	USE stomate_data
32	USE pft_parameters
33	USE constantes
34	USE constantes_soil_var
35
36	IMPLICIT NONE
37
38	PRIVATE
39	PUBLIC gross_glcc_firstday_SinAgeC_fh, gross_glcchange_SinAgeC_fh
40
41	CONTAINS
42
43	! ================================================================================================================================
44	!! SUBROUTINE : harvest_forest
45	!!
46	!>\BRIEF : Handle forest harvest before its legacy is transferred to
47	! newly initialized youngest-age-class PFT.
48	!!
49	!>\DESCRIPTION
50	!_ ================================================================================================================================
51	!!++TEMP++ biomass,veget_frac are not used because the remaining biomass to be
52	!! harvested is calculated within the deforestation fire module.
53	SUBROUTINE harvest_forest (npts,ipts,ivm,biomass,frac, &
54	litter, deforest_biomass_remain,&
55	fuel_1hr,fuel_10hr,&
56	fuel_100hr,fuel_1000hr,&
57	lignin_struc,&
58	bm_to_litter_pro,convflux,prod10,prod100,&
59	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
60	fuel_1000hr_pro, lignin_content_pro)
61
62
63	IMPLICIT NONE
64
65	!! 0.1 Input variables
66	INTEGER, INTENT(in) :: npts
67	INTEGER, INTENT(in) :: ipts
68	INTEGER, INTENT(in) :: ivm
69	REAL(r_std), INTENT(in) :: frac !! the fraction of land covered by forest to be deforested
70	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
71	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
72	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
73	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
74	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
75	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! Vegetmax-weighted remaining litter on the ground for
76	!! deforestation region.
77	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
78	!! deforestation region.
79	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
80	!! above and below ground
81
82	!! 0.2 Modified variables
83	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro !! conversion of biomass to litter
84	!! @tex ($gC m^{-2} day^{-1}$) @endtex
85	REAL(r_std), DIMENSION(:), INTENT(inout) :: convflux !! release during first year following land cover
86	!! change
87
88	REAL(r_std), DIMENSION(npts,0:10), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
89	!! pool after the annual release for each
90	!! compartment (10 + 1 : input from year of land
91	!! cover change)
92	REAL(r_std), DIMENSION(npts,0:100), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
93	!! pool after the annual release for each
94	!! compartment (100 + 1 : input from year of land
95	!! cover change)
96
97	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: litter_pro
98	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1hr_pro
99	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_10hr_pro
100	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_100hr_pro
101	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1000hr_pro
102	REAL(r_std), DIMENSION(:),INTENT(inout) :: lignin_content_pro
103
104
105
106	!! 0.4 Local variables
107	REAL(r_std) :: above
108
109	! harvest of aboveground sap- and heartwood biomass after taking into
110	! account of deforestation fire
111	IF (allow_deforest_fire) THEN
112	above = deforest_biomass_remain(ipts,ivm,isapabove,icarbon)+ &
113	deforest_biomass_remain(ipts,ivm,iheartabove,icarbon)
114	convflux(ipts) = convflux(ipts) + 0
115	prod10(ipts,0) = prod10(ipts,0) + 0.4*above
116	prod100(ipts,0) = prod100(ipts,0) + 0.6*above
117	ELSE
118	above = (biomass(ipts,ivm,isapabove,icarbon)+ &
119	biomass(ipts,ivm,iheartabove,icarbon))*frac
120	convflux(ipts) = convflux(ipts) + coeff_lcchange_1(ivm) * above
121	prod10(ipts,0) = prod10(ipts,0) + coeff_lcchange_10(ivm) * above
122	prod100(ipts,0) = prod100(ipts,0) + coeff_lcchange_100(ivm) * above
123	ENDIF
124
125	! the transfer of dead biomass to litter
126	bm_to_litter_pro(isapbelow,:) = bm_to_litter_pro(isapbelow,:) + &
127	biomass(ipts,ivm,isapbelow,:)*frac
128	bm_to_litter_pro(iheartbelow,:) = bm_to_litter_pro(iheartbelow,:) + &
129	biomass(ipts,ivm,iheartbelow,:)*frac
130	bm_to_litter_pro(iroot,:) = bm_to_litter_pro(iroot,:) + &
131	biomass(ipts,ivm,iroot,:)*frac
132	bm_to_litter_pro(ifruit,:) = bm_to_litter_pro(ifruit,:) + &
133	biomass(ipts,ivm,ifruit,:)*frac
134	bm_to_litter_pro(icarbres,:) = bm_to_litter_pro(icarbres,:) + &
135	biomass(ipts,ivm,icarbres,:)*frac
136	bm_to_litter_pro(ileaf,:) = bm_to_litter_pro(ileaf,:) + &
137	biomass(ipts,ivm,ileaf,:)*frac
138
139	!update litter_pro
140	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
141	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
142	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
143	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
144	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
145	!don't forget to hanle litter lignin content
146	lignin_content_pro(:)= lignin_content_pro(:) + &
147	litter(ipts,istructural,ivm,:,icarbon)fraclignin_struc(ipts,ivm,:)
148
149	END SUBROUTINE harvest_forest
150
151	! ================================================================================================================================
152	!! SUBROUTINE : harvest_herb
153	!!
154	!>\BRIEF : Handle herbaceous PFT clearing before its legacy is transferred to
155	! newly initialized youngest-age-class PFT.
156	!!
157	!>\DESCRIPTION
158	!_ ================================================================================================================================
159	SUBROUTINE harvest_herb (ipts,ivm,biomass,veget_frac,bm_to_litter_pro)
160
161	IMPLICIT NONE
162
163	!! 0.1 Input variables
164	INTEGER, INTENT(in) :: ipts
165	INTEGER, INTENT(in) :: ivm
166	REAL(r_std), INTENT(in) :: veget_frac !! the fraction of land covered by herbaceous PFT to be cleared
167	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
168
169	!! 0.2 Modified variables
170	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro
171
172
173
174	! the transfer of dead biomass to litter
175	bm_to_litter_pro(:,:) = bm_to_litter_pro(:,:) + biomass(ipts,ivm,:,:)*veget_frac
176
177	END SUBROUTINE harvest_herb
178
179
180	! ================================================================================================================================
181	!! SUBROUTINE : initialize_proxy_pft
182	!!
183	!>\BRIEF Initialize a proxy new youngest age class PFT.
184	!!
185	!>\DESCRIPTION Initialize a proxy new youngest age class PFT that will be
186	!! merged with existing yongest age class, or fill the empty
187	!! niche of the youngest age class PFT.
188	!_ ================================================================================================================================
189	SUBROUTINE initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro, &
190	biomass_pro, co2_to_bm_pro, ind_pro, age_pro, &
191	senescence_pro, PFTpresent_pro, &
192	lm_lastyearmax_pro, everywhere_pro, npp_longterm_pro, &
193	leaf_frac_pro,leaf_age_pro)
194
195	IMPLICIT NONE
196
197	!! 0.1 Input variables
198	INTEGER, INTENT(in) :: ipts !!
199	INTEGER, INTENT(in) :: ipft_young_agec !! index of the concerned youngest-age-class PFT
200	REAL(r_std), INTENT(in) :: veget_max_pro !! fraction of grid cell land area that's to be occupied
201
202	!! 0.2 Modified variables
203	REAL(r_std), INTENT(inout) :: co2_to_bm_pro
204
205	!! 0.3 Output variables
206	REAL(r_std), DIMENSION(:,:), INTENT(out) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
207	REAL(r_std), DIMENSION(:), INTENT(out) :: leaf_frac_pro !! fraction of leaves in leaf age class
208	REAL(r_std), DIMENSION(:), INTENT(out) :: leaf_age_pro !! fraction of leaves in leaf age class
209	REAL(r_std), INTENT(out) :: age_pro, ind_pro, lm_lastyearmax_pro
210	REAL(r_std), INTENT(out) :: npp_longterm_pro
211	REAL(r_std), INTENT(out) :: everywhere_pro !! is the PFT everywhere in the grid box or very
212	LOGICAL, INTENT(out) :: senescence_pro !! plant senescent (only for deciduous trees) Set
213	!! to .FALSE. if PFT is introduced or killed
214	LOGICAL, INTENT(out) :: PFTpresent_pro !! Is pft there (unitless)
215
216	!! 0.4 Local variables
217	!REAL(r_std), DIMENSION(npts,nvm) :: when_growthinit !! how many days ago was the beginning of the
218	! !! growing season (days)
219
220	REAL(r_std), DIMENSION(nparts,nelements) :: bm_new !! biomass increase @tex ($gC m^{-2}$) @endtex
221	REAL(r_std) :: cn_ind,ind
222	INTEGER :: i,j,k,l
223
224	! -Note-
225	! This part of codes are copied from the original lcchange_main subroutine
226	! that initialize a new PFT.
227
228	i=ipts
229	j=ipft_young_agec
230
231	!! Initialization of some variables
232	leaf_frac_pro(:) = zero
233	leaf_age_pro(:) = zero
234
235	!! Initial setting of new establishment
236	IF (is_tree(j)) THEN
237	! cn_sapl(j)=0.5; stomate_data.f90
238	cn_ind = cn_sapl(j)
239	ELSE
240	cn_ind = un
241	ENDIF
242	ind = veget_max_pro / cn_ind
243	ind_pro = ind*veget_max_pro
244	PFTpresent_pro = .TRUE.
245	senescence_pro = .FALSE.
246	everywhere_pro = 1.*veget_max_pro
247	age_pro = zero
248
249	! large_value = 1.E33_r_std
250	! when_growthinit(i,j) = large_value
251	leaf_frac_pro(1) = 1.0 * veget_max_pro
252	leaf_age_pro(1) = 1.0 * veget_max_pro !This was not included in original lcchange_main subroutine
253	npp_longterm_pro = npp_longterm_init * veget_max_pro
254	lm_lastyearmax_pro = bm_sapl(j,ileaf,icarbon) * ind * veget_max_pro
255
256	!! Update of biomass in each each carbon stock component (leaf, sapabove, sapbelow,
257	!> heartabove, heartbelow, root, fruit, and carbres)\n
258	DO k = 1, nparts ! loop over # carbon stock components, nparts = 8; stomate_constant.f90
259	DO l = 1,nelements ! loop over # elements
260	biomass_pro(k,l) = ind * bm_sapl(j,k,l)
261	END DO ! loop over # elements
262	co2_to_bm_pro = co2_to_bm_pro + ind * bm_sapl(j,k,icarbon)
263	ENDDO ! loop over # carbon stock components
264
265	END SUBROUTINE initialize_proxy_pft
266
267	! ================================================================================================================================
268	!! SUBROUTINE sap_take
269	!!
270	!>\BRIEF : Take the sapling biomass of the new PFTs from the existing biomass, otherwise
271	! take from co2_to_bm
272	!!
273	!>\DESCRIPTION
274	!_ ================================================================================================================================
275	SUBROUTINE sap_take (ipts,ivma,veget_max,biomass_pro,biomass,co2_to_bm_pro)
276
277	INTEGER, INTENT(in) :: ipts !!
278	INTEGER, INTENT(in) :: ivma
279	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT (LAI ->
280	REAL(r_std), DIMENSION(:,:), INTENT(in) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
281
282	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
283	REAL(r_std), INTENT(inout) :: co2_to_bm_pro
284
285
286	REAL(r_std), DIMENSION(nparts,nelements) :: biomass_total !! biomass @tex ($gC m^{-2}$) @endtex
287	REAL(r_std) :: bm_org,bmpro_share
288	INTEGER :: i,ivm,ipart
289
290	biomass_total(:,:) = zero
291	bm_org = zero
292	bmpro_share = zero
293
294	DO i = 1,nagec_pft(ivma)
295	ivm = start_index(ivma)+i-1
296	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
297	biomass_total = biomass_total + biomass(ipts,ivm,:,:)*veget_max(ipts,ivm)
298	ENDIF
299	ENDDO
300
301	DO ipart = 1, nparts
302	IF (biomass_total(ipart,icarbon) .GT. biomass_pro(ipart,icarbon)) THEN
303	co2_to_bm_pro = co2_to_bm_pro - biomass_pro(ipart,icarbon)
304	!treat each PFT of the MTC
305	DO i = 1,nagec_pft(ivma)
306	ivm = start_index(ivma)+i-1
307	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
308	bm_org = biomass(ipts,ivm,ipart,icarbon) * veget_max(ipts,ivm)
309	bmpro_share = bm_org/biomass_total(ipart,icarbon) * biomass_pro(ipart,icarbon)
310	biomass(ipts,ivm,ipart,icarbon) = (bm_org - bmpro_share)/veget_max(ipts,ivm)
311	ENDIF
312	ENDDO
313	ENDIF
314	ENDDO
315
316	END SUBROUTINE
317
318	! ================================================================================================================================
319	!! SUBROUTINE collect_legacy_pft
320	!!
321	!>\BRIEF : Collect the legacy variables that are going to be included
322	! in the newly initialized PFT.
323	!!
324	!>\DESCRIPTION
325	!_ ================================================================================================================================
326	SUBROUTINE collect_legacy_pft(npts, ipts, ivma, glcc_pftmtc, &
327	biomass, bm_to_litter, carbon, litter, &
328	deepC_a, deepC_s, deepC_p, &
329	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
330	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
331	resp_maint, resp_growth, resp_hetero, co2_fire, &
332	def_fuel_1hr_remain, def_fuel_10hr_remain, &
333	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
334	deforest_litter_remain, deforest_biomass_remain, &
335	veget_max_pro, carbon_pro, lignin_struc_pro, litter_pro, &
336	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
337	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
338	bm_to_litter_pro, co2_to_bm_pro, gpp_daily_pro, &
339	npp_daily_pro, resp_maint_pro, resp_growth_pro, &
340	resp_hetero_pro, co2_fire_pro, &
341	convflux,prod10,prod100)
342
343	IMPLICIT NONE
344
345	!! 0.1 Input variables
346	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
347	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
348	INTEGER, INTENT(in) :: ivma !! Index for metaclass
349	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
350	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
351	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: bm_to_litter !! Transfer of biomass to litter
352	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
353	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: carbon !! carbon pool: active, slow, or passive
354	!! @tex ($gC m^{-2}$) @endtex
355	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_a !! Permafrost soil carbon (g/m**3) active
356	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
357	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
358	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! metabolic and structural litter, above and
359	!! below ground @tex ($gC m^{-2}$) @endtex
360	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
361	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
362	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
363	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
364	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
365	!! above and below ground
366	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_to_bm !! biomass uptaken
367	!! @tex ($gC m^{-2} day^{-1}$) @endtex
368	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gpp_daily !! Daily gross primary productivity
369	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
370	REAL(r_std), DIMENSION(:,:), INTENT(in) :: npp_daily !! Net primary productivity
371	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
372	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_maint !! Maintenance respiration
373	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
374	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_growth !! Growth respiration
375	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
376	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_hetero !! Heterotrophic respiration
377	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
378	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_fire !! Heterotrophic respiration
379	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
380	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_1hr_remain
381	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_10hr_remain
382	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_100hr_remain
383	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_1000hr_remain
384	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
385	!! deforestation region.
386	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
387	!! deforestation region.
388
389	!! 0.2 Output variables
390	REAL(r_std), DIMENSION(:), INTENT(out) :: carbon_pro
391	REAL(r_std), DIMENSION(:), INTENT(out) :: deepC_a_pro
392	REAL(r_std), DIMENSION(:), INTENT(out) :: deepC_s_pro
393	REAL(r_std), DIMENSION(:), INTENT(out) :: deepC_p_pro
394	REAL(r_std), DIMENSION(:), INTENT(out) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
395	!! above and below ground
396	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: litter_pro
397	REAL(r_std), DIMENSION(:,:), INTENT(out) :: fuel_1hr_pro
398	REAL(r_std), DIMENSION(:,:), INTENT(out) :: fuel_10hr_pro
399	REAL(r_std), DIMENSION(:,:), INTENT(out) :: fuel_100hr_pro
400	REAL(r_std), DIMENSION(:,:), INTENT(out) :: fuel_1000hr_pro
401	REAL(r_std), DIMENSION(:,:), INTENT(out) :: bm_to_litter_pro
402	REAL(r_std), INTENT(out) :: veget_max_pro, co2_to_bm_pro
403	REAL(r_std), INTENT(out) :: gpp_daily_pro, npp_daily_pro
404	REAL(r_std), INTENT(out) :: resp_maint_pro, resp_growth_pro
405	REAL(r_std), INTENT(out) :: resp_hetero_pro, co2_fire_pro
406
407	!! 0.3 Modified variables
408	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: convflux !! release during first year following land cover
409	!! change
410
411	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
412	!! pool after the annual release for each
413	!! compartment (10 + 1 : input from year of land
414	!! cover change)
415	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
416	!! pool after the annual release for each
417	!! compartment (100 + 1 : input from year of land
418	!! cover change)
419
420	!! 0.4 Local variables
421	REAL(r_std), DIMENSION(nlevs) :: lignin_content_pro
422	REAL(r_std) :: frac
423	INTEGER :: ivm
424
425
426	! All *_pro variables collect the legacy pools/fluxes of the ancestor
427	! PFTs for the receiving youngest age class. All *_pro variables
428	! represent the quantity weighted by the fraction of ancestor contributing
429	! PFTs.
430	! Exceptions:
431	! lignin_struc_pro:: the ratio of lignin content in structural litter.
432
433	veget_max_pro=zero
434	carbon_pro(:)=zero
435	deepC_a_pro(:)=zero
436	deepC_s_pro(:)=zero
437	deepC_p_pro(:)=zero
438	lignin_struc_pro(:)=zero
439	lignin_content_pro(:)=zero
440	litter_pro(:,:,:)=zero
441	fuel_1hr_pro(:,:)=zero
442	fuel_10hr_pro(:,:)=zero
443	fuel_100hr_pro(:,:)=zero
444	fuel_1000hr_pro(:,:)=zero
445	bm_to_litter_pro(:,:)=zero
446	co2_to_bm_pro=zero
447	gpp_daily_pro=zero
448	npp_daily_pro=zero
449	resp_maint_pro=zero
450	resp_growth_pro=zero
451	resp_hetero_pro=zero
452	co2_fire_pro=zero
453
454	DO ivm = 1,nvm
455	frac = glcc_pftmtc(ipts,ivm,ivma)
456	IF (frac>zero) THEN
457	veget_max_pro = veget_max_pro+frac
458
459	IF (is_tree(ivm)) THEN
460	IF (is_tree(start_index(ivma))) THEN
461	CALL harvest_forest (npts,ipts,ivm,biomass,frac, &
462	litter, deforest_biomass_remain,&
463	fuel_1hr,fuel_10hr,&
464	fuel_100hr,fuel_1000hr,&
465	lignin_struc,&
466	bm_to_litter_pro,convflux(:,iwphar),prod10(:,:,iwphar),prod100(:,:,iwphar),&
467	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
468	fuel_1000hr_pro, lignin_content_pro)
469	ELSE
470	CALL harvest_forest (npts,ipts,ivm,biomass,frac, &
471	litter, deforest_biomass_remain,&
472	fuel_1hr,fuel_10hr,&
473	fuel_100hr,fuel_1000hr,&
474	lignin_struc,&
475	bm_to_litter_pro,convflux(:,iwplcc),prod10(:,:,iwplcc),prod100(:,:,iwplcc),&
476	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
477	fuel_1000hr_pro, lignin_content_pro)
478	ENDIF
479	ELSE
480	CALL harvest_herb(ipts,ivm,biomass,frac, &
481	bm_to_litter_pro)
482	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
483	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
484	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
485	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
486	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
487	!don't forget to hanle litter lignin content
488	lignin_content_pro(:)= lignin_content_pro(:) + &
489	litter(ipts,istructural,ivm,:,icarbon)lignin_struc(ipts,ivm,:)frac
490	ENDIF
491
492	!! scalar variables to be accumulated and inherited
493	!! by the destination PFT
494	bm_to_litter_pro(:,:) = bm_to_litter_pro(:,:) + &
495	bm_to_litter(ipts,ivm,:,:)*frac
496	carbon_pro(:) = carbon_pro(:)+carbon(ipts,:,ivm)*frac
497	deepC_a_pro(:) = deepC_a_pro(:)+deepC_a(ipts,:,ivm)*frac
498	deepC_s_pro(:) = deepC_s_pro(:)+deepC_s(ipts,:,ivm)*frac
499	deepC_p_pro(:) = deepC_p_pro(:)+deepC_p(ipts,:,ivm)*frac
500	co2_to_bm_pro = co2_to_bm_pro + co2_to_bm(ipts,ivm)*frac
501
502	gpp_daily_pro = gpp_daily_pro + gpp_daily(ipts,ivm)*frac
503	npp_daily_pro = npp_daily_pro + npp_daily(ipts,ivm)*frac
504	resp_maint_pro = resp_maint_pro + resp_maint(ipts,ivm)*frac
505	resp_growth_pro = resp_growth_pro + resp_growth(ipts,ivm)*frac
506	resp_hetero_pro = resp_hetero_pro + resp_hetero(ipts,ivm)*frac
507	co2_fire_pro = co2_fire_pro + co2_fire(ipts,ivm)*frac
508	ENDIF
509	ENDDO
510
511	WHERE (litter_pro(istructural,:,icarbon) .GT. min_stomate)
512	lignin_struc_pro(:) = lignin_content_pro(:)/litter_pro(istructural,:,icarbon)
513	ENDWHERE
514
515	END SUBROUTINE collect_legacy_pft
516
517
518	! ================================================================================================================================
519	!! SUBROUTINE gross_lcchange
520	!!
521	!>\BRIEF : Apply gross land cover change.
522	!!
523	!>\DESCRIPTION
524	!_ ================================================================================================================================
525	SUBROUTINE gross_glcchange_SinAgeC_fh (npts, dt_days, harvest_matrix, &
526	glccSecondShift,glccPrimaryShift,glccNetLCC,&
527	def_fuel_1hr_remain, def_fuel_10hr_remain, &
528	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
529	deforest_litter_remain, deforest_biomass_remain, &
530	convflux, cflux_prod10, cflux_prod100, &
531	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
532	veget_max, prod10, prod100, flux10, flux100, &
533	PFTpresent, senescence, moiavail_month, moiavail_week, &
534	gpp_week, ngd_minus5, resp_maint, resp_growth, &
535	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
536	ind, lm_lastyearmax, everywhere, age, &
537	co2_to_bm, gpp_daily, co2_fire, &
538	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
539	gdd_m5_dormance, ncd_dormance, &
540	lignin_struc, carbon, leaf_frac, &
541	deepC_a, deepC_s, deepC_p, &
542	leaf_age, bm_to_litter, biomass, litter, &
543	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
544
545	IMPLICIT NONE
546
547	!! 0.1 Input variables
548
549	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
550	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
551	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccSecondShift !! the land-cover-change (LCC) matrix in case a gross LCC is
552	!! used.
553	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccPrimaryShift !! the land-cover-change (LCC) matrix in case a gross LCC is
554	!! used.
555	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
556	!! used.
557	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: harvest_matrix !!
558	!!
559
560	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1hr_remain
561	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_10hr_remain
562	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_100hr_remain
563	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1000hr_remain
564	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
565	!! deforestation region.
566	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
567	!! deforestation region.
568
569
570	!! 0.2 Output variables
571	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: convflux !! release during first year following land cover
572	!! change
573	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: cflux_prod10 !! total annual release from the 10 year-turnover
574	!! pool @tex ($gC m^{-2}$) @endtex
575	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: cflux_prod100 !! total annual release from the 100 year-
576	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
577	!! after considering the consistency between presribed
578	!! glcc matrix and existing vegetation fractions.
579	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
580	!! there are not enough fractions in the source PFTs
581	!! /vegetations to target PFTs/vegetations. I.e., these
582	!! fraction transfers are presribed in LCC matrix but
583	!! not realized.
584	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
585	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout):: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
586	!! i.e., the contribution of each PFT to the youngest age-class of MTC
587
588	!! 0.3 Modified variables
589	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT (LAI ->
590	!! infinity) on ground (unitless)
591	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
592	!! pool after the annual release for each
593	!! compartment (10 + 1 : input from year of land
594	!! cover change)
595	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
596	!! pool after the annual release for each
597	!! compartment (100 + 1 : input from year of land
598	!! cover change)
599	REAL(r_std), DIMENSION(npts,10,nwp), INTENT(inout) :: flux10 !! annual release from the 10/100 year-turnover
600	!! pool compartments
601	REAL(r_std), DIMENSION(npts,100,nwp), INTENT(inout) :: flux100 !! annual release from the 10/100 year-turnover
602	!! pool compartments
603	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
604	!! each pixel
605	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
606	!! for deciduous trees)
607	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
608	!! unitless)
609	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
610	!! (0 to 1, unitless)
611	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
612	!! @tex $(gC m^{-2} day^{-1})$ @endtex
613	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
614	!! -5 deg C (for phenology)
615	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
616	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
617	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
618	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
619	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
620	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
621	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
622	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
623	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
624	!! the growing season (days)
625	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
626	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
627	!! @tex $(m^{-2})$ @endtex
628	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
629	!! @tex ($gC m^{-2}$) @endtex
630	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
631	!! very localized (after its introduction) (?)
632	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
633	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
634	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
635	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
636	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
637	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_fire !! Fire carbon emissions
638	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
639	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
640	!! availability (days)
641	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
642	!! (for phenology) - this is written to the
643	!! history files
644	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
645	!! for crops
646	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
647	!! C (for phenology)
648	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
649	!! leaves were lost (for phenology)
650	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
651	!! above and below ground
652	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
653	!! @tex ($gC m^{-2}$) @endtex
654	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
655	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
656	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
657	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
658	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
659	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
660	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
661	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
662	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
663	!! below ground @tex ($gC m^{-2}$) @endtex
664	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
665	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
666	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
667	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
668
669	!! 0.4 Local variables
670	REAL(r_std), DIMENSION(nparts,nelements) :: bm_to_litter_pro !! conversion of biomass to litter
671	!! @tex ($gC m^{-2} day^{-1}$) @endtex
672	REAL(r_std), DIMENSION(nparts,nelements) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
673	REAL(r_std) :: veget_max_pro !! "maximal" coverage fraction of a PFT (LAI ->
674	!! infinity) on ground (unitless)
675	REAL(r_std), DIMENSION(ncarb) :: carbon_pro !! carbon pool: active, slow, or passive
676	!! @tex ($gC m^{-2}$) @endtex
677	REAL(r_std), DIMENSION(ndeep) :: deepC_a_pro !! Permafrost carbon pool: active, slow, or passive
678	!! @tex ($gC m^{-3}$) @endtex
679	REAL(r_std), DIMENSION(ndeep) :: deepC_s_pro !! Permafrost carbon pool: active, slow, or passive
680	!! @tex ($gC m^{-3}$) @endtex
681	REAL(r_std), DIMENSION(ndeep) :: deepC_p_pro !! Permafrost carbon pool: active, slow, or passive
682	!! @tex ($gC m^{-3}$) @endtex
683	REAL(r_std), DIMENSION(nlitt,nlevs,nelements) :: litter_pro !! metabolic and structural litter, above and
684	!! below ground @tex ($gC m^{-2}$) @endtex
685	REAL(r_std), DIMENSION(nlitt,nelements) :: fuel_1hr_pro
686	REAL(r_std), DIMENSION(nlitt,nelements) :: fuel_10hr_pro
687	REAL(r_std), DIMENSION(nlitt,nelements) :: fuel_100hr_pro
688	REAL(r_std), DIMENSION(nlitt,nelements) :: fuel_1000hr_pro
689	REAL(r_std), DIMENSION(nlevs) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
690	!! above and below ground
691	REAL(r_std), DIMENSION(nleafages) :: leaf_frac_pro !! fraction of leaves in leaf age class
692	REAL(r_std), DIMENSION(nleafages) :: leaf_age_pro !! fraction of leaves in leaf age class
693	LOGICAL :: PFTpresent_pro, senescence_pro !! Is pft there (unitless)
694	REAL(r_std) :: ind_pro, age_pro, lm_lastyearmax_pro, npp_longterm_pro
695	REAL(r_std) :: everywhere_pro
696	REAL(r_std) :: gpp_daily_pro, npp_daily_pro, co2_to_bm_pro
697	REAL(r_std) :: resp_maint_pro, resp_growth_pro
698	REAL(r_std) :: resp_hetero_pro, co2_fire_pro
699
700	INTEGER :: ipts,ivm,ivma,l,m,ipft_young_agec
701	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
702
703	REAL(r_std), DIMENSION(npts,nvmap) :: glcc_mtc !! Increase in fraction of each MTC in its youngest age-class
704	REAL(r_std), DIMENSION(npts,nvm) :: glccReal_tmp !! A temporary variable to hold glccReal
705	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf_final !!
706	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf_final !!
707	REAL(r_std), DIMENSION(npts) :: pf2yf_compen_sf2yf !!
708	REAL(r_std), DIMENSION(npts) :: sf2yf_compen_pf2yf !!
709	REAL(r_std), DIMENSION(npts,nvm) :: glcc_harvest !! Loss of fraction due to forestry harvest
710
711	WRITE(numout,*) 'Entering gross_lcchange_SinAgeC_fh'
712	glcc_harvest(:,:) = zero
713	glccReal_tmp(:,:) = zero
714
715	!! Some initialization
716	convflux(:,:)=zero
717	prod10(:,0,:) = zero
718	prod100(:,0,:) = zero
719	cflux_prod10(:,:) = zero
720	cflux_prod100(:,:) = zero
721
722	CALL gross_glcc_firstday_SinAgeC_fh(npts,veget_max,harvest_matrix, &
723	glccSecondShift,glccPrimaryShift,glccNetLCC,&
724	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
725	Deficit_pf2yf_final, Deficit_sf2yf_final, &
726	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
727
728	glcc_mtc(:,:) = SUM(glcc_pftmtc,DIM=2)
729	DO ipts=1,npts
730	! Note that we assume people don't intentionally change baresoil to
731	! vegetated land.
732	DO ivma = 2,nvmap
733	! we assume only the youngest age class receives the incoming PFT
734	! [chaoyuejoy@gmail.com 2015-08-04] This line is commented to allow
735	! the case of only single age class being handled.
736	IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
737	ipft_young_agec = start_index(ivma)
738
739	! 1. we accumulate the scalar variables that will be inherited
740	! note we don't handle the case of harvesting forest because
741	! we assume glcc_pftmtc(forest->forest) would be zero and this
742	! case won't occur as it's filtered by the condition of
743	! (frac>min_stomate)
744	CALL collect_legacy_pft(npts, ipts, ivma, glcc_pftmtc, &
745	biomass, bm_to_litter, carbon, litter, &
746	deepC_a, deepC_s, deepC_p, &
747	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
748	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
749	resp_maint, resp_growth, resp_hetero, co2_fire, &
750	def_fuel_1hr_remain, def_fuel_10hr_remain, &
751	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
752	deforest_litter_remain, deforest_biomass_remain, &
753	veget_max_pro, carbon_pro, lignin_struc_pro, litter_pro, &
754	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
755	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
756	bm_to_litter_pro, co2_to_bm_pro, gpp_daily_pro, &
757	npp_daily_pro, resp_maint_pro, resp_growth_pro, &
758	resp_hetero_pro, co2_fire_pro, &
759	convflux,prod10,prod100)
760
761	!++TEMP++
762	! Here we substract the outgoing fraction from the source PFT.
763	! If a too small fraction remains in this source PFT, then it is
764	! exhausted, we empty it. The subroutine 'empty_pft' might be
765	! combined with 'collect_legacy_pft', but now we just put it here.
766	DO ivm = 1,nvm
767	IF( glcc_pftmtc(ipts,ivm,ivma)>min_stomate ) THEN
768	veget_max(ipts,ivm) = veget_max(ipts,ivm)-glcc_pftmtc(ipts,ivm,ivma)
769	IF ( veget_max(ipts,ivm)<min_stomate ) THEN
770	CALL empty_pft(ipts, ivm, veget_max, biomass, ind, &
771	carbon, litter, lignin_struc, bm_to_litter, &
772	deepC_a, deepC_s, deepC_p, &
773	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
774	gpp_daily, npp_daily, gpp_week, npp_longterm, &
775	co2_to_bm, resp_maint, resp_growth, resp_hetero, &
776	lm_lastyearmax, leaf_frac, leaf_age, age, &
777	everywhere, PFTpresent, when_growthinit, &
778	senescence, gdd_from_growthinit, gdd_midwinter, &
779	time_hum_min, gdd_m5_dormance, ncd_dormance, &
780	moiavail_month, moiavail_week, ngd_minus5)
781	ENDIF
782	ENDIF
783	ENDDO
784
785	! 2. we establish a proxy PFT with the fraction of veget_max_pro,
786	! which is going to be either merged with existing target
787	! `ipft_young_agec` PFT, or fill the place if no existing target PFT
788	! exits.
789	CALL initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro, &
790	biomass_pro, co2_to_bm_pro, ind_pro, age_pro, &
791	senescence_pro, PFTpresent_pro, &
792	lm_lastyearmax_pro, everywhere_pro, npp_longterm_pro, &
793	leaf_frac_pro,leaf_age_pro)
794
795	CALL sap_take (ipts,ivma,veget_max,biomass_pro,biomass,co2_to_bm_pro)
796
797	! 3. we merge the newly initiazlized proxy PFT into existing one
798	! or use it to fill an empty PFT slot.
799	CALL add_incoming_proxy_pft(npts, ipts, ipft_young_agec, veget_max_pro,&
800	carbon_pro, litter_pro, lignin_struc_pro, bm_to_litter_pro, &
801	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
802	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
803	biomass_pro, co2_to_bm_pro, npp_longterm_pro, ind_pro, &
804	lm_lastyearmax_pro, age_pro, everywhere_pro, &
805	leaf_frac_pro, leaf_age_pro, PFTpresent_pro, senescence_pro, &
806	gpp_daily_pro, npp_daily_pro, resp_maint_pro, resp_growth_pro, &
807	resp_hetero_pro, co2_fire_pro, &
808	veget_max, carbon, litter, lignin_struc, bm_to_litter, &
809	deepC_a, deepC_s, deepC_p, &
810	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
811	biomass, co2_to_bm, npp_longterm, ind, &
812	lm_lastyearmax, age, everywhere, &
813	leaf_frac, leaf_age, PFTpresent, senescence, &
814	gpp_daily, npp_daily, resp_maint, resp_growth, &
815	resp_hetero, co2_fire)
816
817	ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
818
819	ENDDO
820	ENDDO
821
822	!! Update 10 year-turnover pool content following flux emission
823	!! (linear decay (10%) of the initial carbon input)
824	DO l = 0, 8
825	m = 10 - l
826	cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,m,:)
827	prod10(:,m,:) = prod10(:,m-1,:) - flux10(:,m-1,:)
828	flux10(:,m,:) = flux10(:,m-1,:)
829	WHERE (prod10(:,m,:) .LT. 1.0) prod10(:,m,:) = zero
830	ENDDO
831
832	cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,1,:)
833	flux10(:,1,:) = 0.1 * prod10(:,0,:)
834	prod10(:,1,:) = prod10(:,0,:)
835
836	!! 2.4.3 update 100 year-turnover pool content following flux emission\n
837	DO l = 0, 98
838	m = 100 - l
839	cflux_prod100(:,:) = cflux_prod100(:,:) + flux100(:,m,:)
840	prod100(:,m,:) = prod100(:,m-1,:) - flux100(:,m-1,:)
841	flux100(:,m,:) = flux100(:,m-1,:)
842
843	WHERE (prod100(:,m,:).LT.1.0) prod100(:,m,:) = zero
844	ENDDO
845
846	cflux_prod100(:,:) = cflux_prod100(:,:) + flux100(:,1,:)
847	flux100(:,1,:) = 0.01 * prod100(:,0,:)
848	prod100(:,1,:) = prod100(:,0,:)
849	prod10(:,0,:) = zero
850	prod100(:,0,:) = zero
851
852	convflux = convflux/one_year*dt_days
853	cflux_prod10 = cflux_prod10/one_year*dt_days
854	cflux_prod100 = cflux_prod100/one_year*dt_days
855
856	! Write out history files
857	CALL histwrite_p (hist_id_stomate, 'glcc_pft', itime, &
858	glcc_pft, npts*nvm, horipft_index)
859
860	glccReal_tmp(:,1:12) = glccReal
861	CALL histwrite_p (hist_id_stomate, 'glccReal', itime, &
862	glccReal_tmp, npts*nvm, horipft_index)
863
864	! Write out forestry harvest variables
865	DO ipts = 1,npts
866	DO ivm = 1,nvm
867	DO ivma = 1,nvmap
868	IF (is_tree(ivm) .AND. is_tree(start_index(ivma))) THEN
869	glcc_harvest(ipts,ivm) = glcc_harvest(ipts,ivm) + glcc_pftmtc(ipts,ivm,ivma)
870	ENDIF
871	ENDDO
872	ENDDO
873	ENDDO
874	CALL histwrite_p (hist_id_stomate, 'glcc_harvest', itime, &
875	glcc_harvest, npts*nvm, horipft_index)
876
877	glccReal_tmp(:,:) = zero
878	glccReal_tmp(:,1:4) = IncreDeficit
879	CALL histwrite_p (hist_id_stomate, 'IncreDeficit', itime, &
880	glccReal_tmp, npts*nvm, horipft_index)
881
882	glccReal_tmp(:,:) = zero
883	glccReal_tmp(:,1) = Deficit_pf2yf_final
884	glccReal_tmp(:,2) = Deficit_sf2yf_final
885	glccReal_tmp(:,3) = pf2yf_compen_sf2yf
886	glccReal_tmp(:,4) = sf2yf_compen_pf2yf
887
888	CALL histwrite_p (hist_id_stomate, 'DefiComForHarvest', itime, &
889	glccReal_tmp, npts*nvm, horipft_index)
890
891	DO ivma = 1, nvmap
892	WRITE(part_str,'(I2)') ivma
893	IF (ivma < 10) part_str(1:1) = '0'
894	CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_'//part_str(1:LEN_TRIM(part_str)), &
895	itime, glcc_pftmtc(:,:,ivma), npts*nvm, horipft_index)
896	ENDDO
897	END SUBROUTINE gross_glcchange_SinAgeC_fh
898
899
900	! ================================================================================================================================
901	!! SUBROUTINE : add_incoming_proxy_pft
902	!!
903	!>\BRIEF : Merge the newly incoming proxy PFT cohort with the exisiting
904	!! cohort.
905	!! \n
906	!
907	!_ ================================================================================================================================
908	SUBROUTINE add_incoming_proxy_pft(npts, ipts, ipft, veget_max_pro, &
909	carbon_pro, litter_pro, lignin_struc_pro, bm_to_litter_pro, &
910	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
911	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
912	biomass_pro, co2_to_bm_pro, npp_longterm_pro, ind_pro, &
913	lm_lastyearmax_pro, age_pro, everywhere_pro, &
914	leaf_frac_pro, leaf_age_pro, PFTpresent_pro, senescence_pro, &
915	gpp_daily_pro, npp_daily_pro, resp_maint_pro, resp_growth_pro, &
916	resp_hetero_pro, co2_fire_pro, &
917	veget_max, carbon, litter, lignin_struc, bm_to_litter, &
918	deepC_a, deepC_s, deepC_p, &
919	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
920	biomass, co2_to_bm, npp_longterm, ind, &
921	lm_lastyearmax, age, everywhere, &
922	leaf_frac, leaf_age, PFTpresent, senescence, &
923	gpp_daily, npp_daily, resp_maint, resp_growth, &
924	resp_hetero, co2_fire)
925
926	IMPLICIT NONE
927
928	!! 0.1 Input variables
929	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
930	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
931	INTEGER, INTENT(in) :: ipft
932	REAL(r_std), INTENT(in) :: veget_max_pro !! The land fraction of incoming new PFTs that are
933	!! the sum of all its ancestor PFTs
934	REAL(r_std), DIMENSION(:), INTENT(in) :: carbon_pro
935	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_a_pro
936	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_s_pro
937	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_p_pro
938	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: litter_pro
939	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_1hr_pro
940	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_10hr_pro
941	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_100hr_pro
942	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_1000hr_pro
943	REAL(r_std), DIMENSION(:,:), INTENT(in) :: bm_to_litter_pro
944	REAL(r_std), DIMENSION(:), INTENT(in) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
945	!! above and below ground
946	REAL(r_std), DIMENSION(:,:), INTENT(in) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
947	REAL(r_std), DIMENSION(:), INTENT(in) :: leaf_frac_pro !! fraction of leaves in leaf age class
948	REAL(r_std), DIMENSION(:), INTENT(in) :: leaf_age_pro !! fraction of leaves in leaf age class
949	REAL(r_std), INTENT(in) :: ind_pro, age_pro, lm_lastyearmax_pro
950	REAL(r_std), INTENT(in) :: npp_longterm_pro, co2_to_bm_pro
951	REAL(r_std), INTENT(in) :: everywhere_pro !! is the PFT everywhere in the grid box or very
952	LOGICAL, INTENT(in) :: PFTpresent_pro, senescence_pro !! Is pft there (unitless)
953
954	REAL(r_std), INTENT(in) :: gpp_daily_pro, npp_daily_pro
955	REAL(r_std), INTENT(in) :: resp_maint_pro, resp_growth_pro
956	REAL(r_std), INTENT(in) :: resp_hetero_pro, co2_fire_pro
957
958	!! 0.2 Output variables
959
960	!! 0.3 Modified variables
961
962	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
963	!! May sum to
964	!! less than unity if the pixel has
965	!! nobio area. (unitless, 0-1)
966
967	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
968	!! @tex ($gC m^{-2}$) @endtex
969	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
970	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
971	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
972	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
973	!! below ground @tex ($gC m^{-2}$) @endtex
974	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
975	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
976	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
977	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
978	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
979	!! above and below ground
980	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
981	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
982	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
983	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
984	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
985
986	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
987	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
988	!! @tex $(m^{-2})$ @endtex
989	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
990	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
991	!! each pixel
992	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
993	!! for deciduous trees)
994	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
995	!! @tex ($gC m^{-2}$) @endtex
996	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
997	!! very localized (after its introduction) (?)
998	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
999	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
1000
1001	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
1002	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1003	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
1004	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1005	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
1006	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1007	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
1008	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1009	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
1010	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1011	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_fire !! Heterotrophic respiration
1012	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1013
1014	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
1015	! !! unitless)
1016	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
1017	! !! (0 to 1, unitless)
1018	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
1019	! !! @tex $(gC m^{-2} day^{-1})$ @endtex
1020	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
1021	! !! -5 deg C (for phenology)
1022	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
1023	! !! the growing season (days)
1024	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
1025	! !! availability (days)
1026	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
1027	! !! (for phenology) - this is written to the
1028	! !! history files
1029	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
1030	! !! for crops
1031	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
1032	! !! C (for phenology)
1033	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
1034	! !! leaves were lost (for phenology)
1035
1036	!! 0.4 Local variables
1037
1038	INTEGER(i_std) :: iele !! Indeces(unitless)
1039	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
1040	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements) :: litter_old !! metabolic and structural litter, above and
1041	!! below ground @tex ($gC m^{-2}$) @endtex
1042	REAL(r_std) :: veget_old,veget_total
1043
1044
1045	! Back up some variables in case they're needed later
1046	litter_old(:,:,:,:,:) = litter(:,:,:,:,:)
1047
1048	!! General idea
1049	! The established proxy vegetation has a fraction of 'veget_max_pro'; the
1050	! existing iPFT has a fraction of veget_max(ipts,ipft).
1051	! Suppose we want to merge a scalar variable B, the value of B after merging
1052	! is (BiVi+BjVj)/(Vi+Vj), where Vi is the original veget_max, Vj is the
1053	! incoming veget_max. Note that in case Vi=0, this equation remains solid,
1054	! i.e. the veget_max after merging is Vj and B after merging is Bj. In other
1055	! words, the proxy vegetation "fills" up the empty niche of iPFT.
1056	! Also note that for many scalar variables our input value is Bj*Vj, which
1057	! is accumulated from multiple ancestor PFTs.
1058	veget_old = veget_max(ipts,ipft)
1059	veget_total = veget_old+veget_max_pro
1060
1061	!! Different ways of handling merging depending on nature of variables:
1062
1063	!! 1. Area-based scalar variables, use the equation above
1064	! biomass,carbon, litter, bm_to_litter, co2_to_bm, ind,
1065	! lm_lastyearmax, npp_longterm, lm_lastyearmax,
1066	! lignin_struc (ratio variable depending on area-based variable)
1067
1068	!! 2. Variables are tentatively handled like area-based variables:
1069	! leaf_frac, leaf_age,
1070
1071	!! 3. Variables that are overwritten by the newly initialized PFT:
1072	! PFTpresent, senescence
1073
1074	!! 4. Variables whose operation is uncertain and are not handled currently:
1075	! when_growthinit :: how many days ago was the beginning of the growing season (days)
1076	! gdd_from_growthinit :: growing degree days, since growthinit
1077	! gdd_midwinter, time_hum_min, gdd_m5_dormance, ncd_dormance,
1078	! moiavail_month, moiavail_week, ngd_minus5
1079
1080	!! 5. Variables that concern with short-term fluxes that do not apply in
1081	! this case:
1082	! gpp_daily, npp_daily etc.
1083
1084	! Add the coming veget_max_pro into existing veget_max
1085	veget_max(ipts,ipft) = veget_total
1086
1087	! Merge scalar variables which are defined on area basis
1088	carbon(ipts,:,ipft) = (veget_old * carbon(ipts,:,ipft) + &
1089	carbon_pro(:))/veget_total
1090	deepC_a(ipts,:,ipft) = (veget_old * deepC_a(ipts,:,ipft) + &
1091	deepC_a_pro(:))/veget_total
1092	deepC_s(ipts,:,ipft) = (veget_old * deepC_s(ipts,:,ipft) + &
1093	deepC_s_pro(:))/veget_total
1094	deepC_p(ipts,:,ipft) = (veget_old * deepC_p(ipts,:,ipft) + &
1095	deepC_p_pro(:))/veget_total
1096	litter(ipts,:,ipft,:,:) = (veget_old * litter(ipts,:,ipft,:,:) + &
1097	litter_pro(:,:,:))/veget_total
1098	fuel_1hr(ipts,ipft,:,:) = (veget_old * fuel_1hr(ipts,ipft,:,:) + &
1099	fuel_1hr_pro(:,:))/veget_total
1100	fuel_10hr(ipts,ipft,:,:) = (veget_old * fuel_10hr(ipts,ipft,:,:) + &
1101	fuel_10hr_pro(:,:))/veget_total
1102	fuel_100hr(ipts,ipft,:,:) = (veget_old * fuel_100hr(ipts,ipft,:,:) + &
1103	fuel_100hr_pro(:,:))/veget_total
1104	fuel_1000hr(ipts,ipft,:,:) = (veget_old * fuel_1000hr(ipts,ipft,:,:) + &
1105	fuel_1000hr_pro(:,:))/veget_total
1106
1107	WHERE (litter(ipts,istructural,ipft,:,icarbon) .GT. min_stomate)
1108	lignin_struc(ipts,ipft,:) = (veget_oldlitter_old(ipts,istructural,ipft,:,icarbon) &
1109	lignin_struc(ipts,ipft,:) + litter_pro(istructural,:,icarbon)* &
1110	lignin_struc_pro(:))/(veget_total*litter(ipts,istructural,ipft,:,icarbon))
1111	ENDWHERE
1112	bm_to_litter(ipts,ipft,:,:) = (veget_old * bm_to_litter(ipts,ipft,:,:) + &
1113	bm_to_litter_pro(:,:))/veget_total
1114
1115	biomass(ipts,ipft,:,:) = (biomass(ipts,ipft,:,:)*veget_old + &
1116	biomass_pro(:,:))/veget_total
1117	co2_to_bm(ipts,ipft) = (veget_old*co2_to_bm(ipts,ipft) + &
1118	co2_to_bm_pro)/veget_total
1119	ind(ipts,ipft) = (ind(ipts,ipft)*veget_old + ind_pro)/veget_total
1120	lm_lastyearmax(ipts,ipft) = (lm_lastyearmax(ipts,ipft)*veget_old + &
1121	lm_lastyearmax_pro)/veget_total
1122	npp_longterm(ipts,ipft) = (veget_old * npp_longterm(ipts,ipft) + &
1123	npp_longterm_pro)/veget_total
1124
1125	!CHECK: Here follows the original idea in DOFOCO, more strictly,
1126	! leas mass should be considered together. The same also applies on
1127	! leaf age.
1128	leaf_frac(ipts,ipft,:) = (leaf_frac(ipts,ipft,:)*veget_old + &
1129	leaf_frac_pro(:))/veget_total
1130	leaf_age(ipts,ipft,:) = (leaf_age(ipts,ipft,:)*veget_old + &
1131	leaf_age_pro(:))/veget_total
1132	age(ipts,ipft) = (veget_old * age(ipts,ipft) + &
1133	age_pro)/veget_total
1134
1135	! Everywhere deals with the migration of vegetation. Copy the
1136	! status of the most migrated vegetation for the whole PFT
1137	everywhere(ipts,ipft) = MAX(everywhere(ipts,ipft), everywhere_pro)
1138
1139	! Overwrite the original variables with that from newly initialized
1140	! proxy PFT
1141	PFTpresent(ipts,ipft) = PFTpresent_pro
1142	senescence(ipts,ipft) = senescence_pro
1143
1144	! This is to close carbon loop when writing history variables.
1145	gpp_daily(ipts,ipft) = (veget_old * gpp_daily(ipts,ipft) + &
1146	gpp_daily_pro)/veget_total
1147	npp_daily(ipts,ipft) = (veget_old * npp_daily(ipts,ipft) + &
1148	npp_daily_pro)/veget_total
1149	resp_maint(ipts,ipft) = (veget_old * resp_maint(ipts,ipft) + &
1150	resp_maint_pro)/veget_total
1151	resp_growth(ipts,ipft) = (veget_old * resp_growth(ipts,ipft) + &
1152	resp_growth_pro)/veget_total
1153	resp_hetero(ipts,ipft) = (veget_old * resp_hetero(ipts,ipft) + &
1154	resp_hetero_pro)/veget_total
1155	co2_fire(ipts,ipft) = (veget_old * co2_fire(ipts,ipft) + &
1156	co2_fire_pro)/veget_total
1157
1158	! Phenology- or time-related variables will be copied from original values if
1159	! there is already youngest-age-class PFT there, otherwise they're left
1160	! untouched, because 1. to initiliaze all new PFTs here is wrong and
1161	! phenology is not explicitly considered, so we cannot assign a value
1162	! to these variables. 2. We assume they will be correctly filled if
1163	! other variables are in place (e.g., non-zero leaf mass will lead to
1164	! onset of growing season). In this case, merging a newly initialized PFT
1165	! to an existing one is not the same as merging PFTs when they grow
1166	! old enough to exceed thresholds.
1167
1168	! gpp_week(ipts,ipft) = (veget_old * gpp_week(ipts,ipft) + &
1169	! gpp_week_pro)/veget_total
1170	! when_growthinit(ipts,ipft) = (veget_old * when_growthinit(ipts,ipft) + &
1171	! when_growthinit_pro)/veget_total
1172	! gdd_from_growthinit(ipts,ipft) = (veget_old * gdd_from_growthinit(ipts,ipft) + &
1173	! gdd_from_growthinit_pro)/veget_total
1174	! gdd_midwinter(ipts,ipft) = (veget_old * gdd_midwinter(ipts,ipft) + &
1175	! gdd_midwinter_pro)/veget_total
1176	! time_hum_min(ipts,ipft) = (veget_old * time_hum_min(ipts,ipft) + &
1177	! time_hum_min_pro)/veget_total
1178	! gdd_m5_dormance(ipts,ipft) = (veget_old * gdd_m5_dormance(ipts,ipft) + &
1179	! gdd_m5_dormance_pro)/veget_total
1180	! ncd_dormance(ipts,ipft) = (veget_old * ncd_dormance(ipts,ipft) + &
1181	! ncd_dormance_pro)/veget_total
1182	! moiavail_month(ipts,ipft) = (veget_old * moiavail_month(ipts,ipft) + &
1183	! moiavail_month_pro)/veget_total
1184	! moiavail_week(ipts,ipft) = (veget_old * moiavail_week(ipts,ipft) + &
1185	! moiavail_week_pro)/veget_total
1186	! ngd_minus5(ipts,ipft) = (veget_old * ngd_minus5(ipts,ipft) + &
1187	! ngd_minus5_pro)/veget_total
1188
1189
1190	END SUBROUTINE add_incoming_proxy_pft
1191
1192
1193	! ================================================================================================================================
1194	!! SUBROUTINE : empty_pft
1195	!!
1196	!>\BRIEF : Empty a PFT when,
1197	!! - it is exhausted because of land cover change.
1198	!! - it moves to the next age class
1199	!! \n
1200	!_ ================================================================================================================================
1201	SUBROUTINE empty_pft(ipts, ivm, veget_max, biomass, ind, &
1202	carbon, litter, lignin_struc, bm_to_litter, &
1203	deepC_a, deepC_s, deepC_p, &
1204	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1205	gpp_daily, npp_daily, gpp_week, npp_longterm, &
1206	co2_to_bm, resp_maint, resp_growth, resp_hetero, &
1207	lm_lastyearmax, leaf_frac, leaf_age, age, &
1208	everywhere, PFTpresent, when_growthinit, &
1209	senescence, gdd_from_growthinit, gdd_midwinter, &
1210	time_hum_min, gdd_m5_dormance, ncd_dormance, &
1211	moiavail_month, moiavail_week, ngd_minus5)
1212
1213	IMPLICIT NONE
1214
1215	!! 0.1 Input variables
1216	INTEGER, INTENT(in) :: ipts !! index for grid cell
1217	INTEGER, INTENT(in) :: ivm !! index for pft
1218
1219	!! 0.2 Output variables
1220
1221	!! 0.3 Modified variables
1222
1223	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1224	!! May sum to
1225	!! less than unity if the pixel has
1226	!! nobio area. (unitless, 0-1)
1227	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
1228	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
1229	!! @tex $(m^{-2})$ @endtex
1230	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
1231	!! @tex ($gC m^{-2}$) @endtex
1232	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
1233	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1234	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1235	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
1236	!! below ground @tex ($gC m^{-2}$) @endtex
1237	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
1238	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
1239	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
1240	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
1241	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1242	!! above and below ground
1243	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
1244	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1245	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
1246	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1247	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
1248	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1249	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
1250	!! @tex $(gC m^{-2} day^{-1})$ @endtex
1251	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
1252	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
1253	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
1254	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
1255	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1256	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
1257	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1258	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
1259	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1260	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
1261	!! @tex ($gC m^{-2}$) @endtex
1262	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
1263	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
1264	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
1265	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
1266	!! very localized (after its introduction) (?)
1267	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
1268	!! each pixel
1269	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
1270	!! the growing season (days)
1271	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
1272	!! for deciduous trees)
1273	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
1274	!! for crops
1275	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
1276	!! (for phenology) - this is written to the
1277	!! history files
1278	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
1279	!! availability (days)
1280	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
1281	!! C (for phenology)
1282	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
1283	!! leaves were lost (for phenology)
1284	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
1285	!! unitless)
1286	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
1287	!! (0 to 1, unitless)
1288	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
1289	!! -5 deg C (for phenology)
1290
1291	!! 0.4 Local variables
1292	INTEGER(i_std) :: iele !! Indeces(unitless)
1293	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
1294
1295	veget_max(ipts,ivm) = zero
1296	ind(ipts,ivm) = zero
1297	biomass(ipts,ivm,:,:) = zero
1298	litter(ipts,:,ivm,:,:) = zero
1299	fuel_1hr(ipts,ivm,:,:) = zero
1300	fuel_10hr(ipts,ivm,:,:) = zero
1301	fuel_100hr(ipts,ivm,:,:) = zero
1302	fuel_1000hr(ipts,ivm,:,:) = zero
1303	carbon(ipts,:,ivm) = zero
1304	deepC_a(ipts,:,ivm) = zero
1305	deepC_s(ipts,:,ivm) = zero
1306	deepC_p(ipts,:,ivm) = zero
1307	bm_to_litter(ipts,ivm,:,:) = zero
1308	DO ilev=1,nlevs
1309	lignin_struc(ipts,ivm,ilev) = zero
1310	ENDDO
1311	npp_longterm(ipts,ivm) = zero
1312	gpp_daily(ipts,ivm) = zero
1313	gpp_week(ipts,ivm) = zero
1314	resp_maint(ipts,ivm) = zero
1315	resp_growth(ipts,ivm) = zero
1316	resp_hetero(ipts,ivm) = zero
1317	npp_daily(ipts,ivm) = zero
1318	co2_to_bm(ipts,ivm) = zero
1319	lm_lastyearmax(ipts,ivm) = zero
1320	age(ipts,ivm) = zero
1321	leaf_frac(ipts,ivm,:) = zero
1322	leaf_age(ipts,ivm,:) = zero
1323	everywhere(ipts,ivm) = zero
1324	when_growthinit(ipts,ivm) = zero
1325	gdd_from_growthinit(ipts,ivm) = zero
1326	gdd_midwinter(ipts,ivm) = zero
1327	time_hum_min(ipts,ivm) = zero
1328	gdd_m5_dormance(ipts,ivm) = zero
1329	ncd_dormance(ipts,ivm) = zero
1330	moiavail_month(ipts,ivm) = zero
1331	moiavail_week(ipts,ivm) = zero
1332	ngd_minus5(ipts,ivm) = zero
1333	PFTpresent(ipts,ivm) = .FALSE.
1334	senescence(ipts,ivm) = .FALSE.
1335
1336	END SUBROUTINE empty_pft
1337
1338	! ================================================================================================================================
1339	!! SUBROUTINE : gross_lcc_firstday
1340	!!
1341	!>\BRIEF : When necessary, adjust input glcc matrix, and allocate it
1342	!! into different contributing age classes and receiving
1343	!! youngest age classes.
1344	!! \n
1345	!_ ================================================================================================================================
1346
1347	! Note: it has this name because this subroutine will also be called
1348	! the first day of each year to precalculate the forest loss for the
1349	! deforestation fire module.
1350	SUBROUTINE gross_glcc_firstday_SinAgeC_fh(npts,veget_max_org,harvest_matrix, &
1351	glccSecondShift,glccPrimaryShift,glccNetLCC,&
1352	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
1353	Deficit_pf2yf_final, Deficit_sf2yf_final, &
1354	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
1355
1356	IMPLICIT NONE
1357
1358	!! 0.1 Input variables
1359
1360	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
1361	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_org !! "maximal" coverage fraction of a PFT on the ground
1362	!! May sum to
1363	!! less than unity if the pixel has
1364	!! nobio area. (unitless, 0-1)
1365	REAL(r_std), DIMENSION(npts,12),INTENT(in) :: harvest_matrix !!
1366	!!
1367	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccSecondShift !! the land-cover-change (LCC) matrix in case a gross LCC is
1368	!! used.
1369	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccPrimaryShift !! the land-cover-change (LCC) matrix in case a gross LCC is
1370	!! used.
1371	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
1372	!! used.
1373
1374	!! 0.2 Output variables
1375	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
1376	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
1377	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
1378	!! after considering the consistency between presribed
1379	!! glcc matrix and existing vegetation fractions.
1380	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
1381	!! there are not enough fractions in the source PFTs
1382	!! /vegetations to target PFTs/vegetations. I.e., these
1383	!! fraction transfers are presribed in LCC matrix but
1384	!! not realized.
1385	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Deficit_pf2yf_final !!
1386	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Deficit_sf2yf_final !!
1387	REAL(r_std), DIMENSION(npts), INTENT(inout) :: pf2yf_compen_sf2yf !!
1388	REAL(r_std), DIMENSION(npts), INTENT(inout) :: sf2yf_compen_pf2yf !!
1389
1390
1391	!! 0.3 Modified variables
1392
1393	!! 0.4 Local variables
1394	REAL(r_std), DIMENSION (npts,12) :: glcc !! the land-cover-change (LCC) matrix in case a gross LCC is
1395	!! used.
1396	REAL(r_std), DIMENSION(npts,nvmap) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
1397	REAL(r_std), DIMENSION(npts,nagec_tree) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
1398	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_grass !! fraction of grass age-class groups, in sequence of old->young
1399	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_pasture !! fraction of pasture age-class groups, in sequence of old->young
1400	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_crop !! fraction of crop age-class groups, in sequence of old->young
1401
1402
1403	REAL(r_std), DIMENSION(npts,4) :: veget_4veg !! "maximal" coverage fraction of a PFT on the ground
1404	REAL(r_std), DIMENSION(npts) :: veget_tree !! "maximal" coverage fraction of a PFT on the ground
1405	REAL(r_std), DIMENSION(npts) :: veget_grass !! "maximal" coverage fraction of a PFT on the ground
1406	REAL(r_std), DIMENSION(npts) :: veget_pasture !! "maximal" coverage fraction of a PFT on the ground
1407	REAL(r_std), DIMENSION(npts) :: veget_crop !! "maximal" coverage fraction of a PFT on the ground
1408
1409	REAL(r_std), DIMENSION(npts,nvm) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1410	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_tmp !! "maximal" coverage fraction of a PFT on the ground
1411	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_old !! "maximal" coverage fraction of a PFT on the ground
1412	REAL(r_std), DIMENSION(npts,nvm) :: glcc_pft_tmp !! Loss of fraction in each PFT
1413
1414	! Different indexes for convenient local uses
1415	! We define the rules for gross land cover change matrix:
1416	! 1 forest->grass
1417	! 2 forest->pasture
1418	! 3 forest->crop
1419	! 4 grass->forest
1420	! 5 grass->pasture
1421	! 6 grass->crop
1422	! 7 pasture->forest
1423	! 8 pasture->grass
1424	! 9 pasture->crop
1425	! 10 crop->forest
1426	! 11 crop->grass
1427	! 12 crop->pasture
1428	INTEGER :: f2g=1, f2p=2, f2c=3
1429	INTEGER :: g2f=4, g2p=5, g2c=6, p2f=7, p2g=8, p2c=9, c2f=10, c2g=11, c2p=12
1430
1431	INTEGER, ALLOCATABLE :: indall_tree(:) !! Indices for all tree PFTs
1432	INTEGER, ALLOCATABLE :: indold_tree(:) !! Indices for old tree cohort only
1433	INTEGER, ALLOCATABLE :: indagec_tree(:,:) !! Indices for secondary tree cohorts,
1434	!! note the sequence is old->young.
1435	INTEGER, ALLOCATABLE :: indall_grass(:) !! Indices for all grass PFTs
1436	INTEGER, ALLOCATABLE :: indold_grass(:) !! Indices for old grasses only
1437	INTEGER, ALLOCATABLE :: indagec_grass(:,:) !! Indices for secondary grass cohorts
1438	!! note the sequence is old->young.
1439	INTEGER, ALLOCATABLE :: indall_pasture(:) !! Indices for all pasture PFTs
1440	INTEGER, ALLOCATABLE :: indold_pasture(:) !! Indices for old pasture only
1441	INTEGER, ALLOCATABLE :: indagec_pasture(:,:) !! Indices for secondary pasture cohorts
1442	!! note the sequence is old->young.
1443	INTEGER, ALLOCATABLE :: indall_crop(:) !! Indices for all crop PFTs
1444	INTEGER, ALLOCATABLE :: indold_crop(:) !! Indices for old crops only
1445	INTEGER, ALLOCATABLE :: indagec_crop(:,:) !! Indices for secondary crop cohorts
1446	!! note the sequence is old->young.
1447	INTEGER :: num_tree_sinagec,num_tree_mulagec,num_grass_sinagec,num_grass_mulagec, &
1448	num_pasture_sinagec,num_pasture_mulagec,num_crop_sinagec,num_crop_mulagec, &
1449	itree,itree2,igrass,igrass2,ipasture,ipasture2,icrop,icrop2,pf2yf,sf2yf
1450	INTEGER :: i,j,ivma,staind,endind,ivm
1451
1452
1453	REAL(r_std), DIMENSION(npts,12) :: glccDef !! Gross LCC deficit, negative values mean that there
1454	!! are not enough fractions in the source vegetations
1455	!! to the target ones as presribed by the LCC matrix.
1456	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf !!
1457	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf !!
1458	REAL(r_std), DIMENSION(npts) :: Surplus_pf2yf !!
1459	REAL(r_std), DIMENSION(npts) :: Surplus_sf2yf !!
1460	REAL(r_std), DIMENSION(npts,12) :: HmatrixReal !!
1461	INTEGER :: ipts
1462
1463
1464	!! 1. We first build all different indices that we are going to use
1465	!! in handling the PFT exchanges, three types of indices are built:
1466	!! - for all age classes
1467	!! - include only oldest age classes
1468	!! - include all age classes excpet the oldest ones
1469	! We have to build these indices because we would like to extract from
1470	! donating PFTs in the sequnce of old->young age classes, and add in the
1471	! receving PFTs only in the youngest-age-class PFTs. These indicies allow
1472	! us to know where the different age classes are.
1473
1474	num_tree_sinagec=0 ! number of tree PFTs with only one single age class
1475	! considered as the oldest age class
1476	num_tree_mulagec=0 ! number of tree PFTs having multiple age classes
1477	num_grass_sinagec=0
1478	num_grass_mulagec=0
1479	num_pasture_sinagec=0
1480	num_pasture_mulagec=0
1481	num_crop_sinagec=0
1482	num_crop_mulagec=0
1483
1484	!! 1.1 Calculate the number of PFTs for different MTCs and allocate
1485	!! the old and all indices arrays.
1486
1487	! [Note here the sequence to identify tree,pasture,grass,crop] is
1488	! critical. The similar sequence is used in the subroutine "calc_cover".
1489	! Do not forget to change the sequence there if you modify here.
1490	DO ivma =2,nvmap
1491	staind=start_index(ivma)
1492	IF (nagec_pft(ivma)==1) THEN
1493	IF (is_tree(staind)) THEN
1494	num_tree_sinagec = num_tree_sinagec+1
1495	ELSE IF (is_grassland_manag(staind)) THEN
1496	num_pasture_sinagec = num_pasture_sinagec+1
1497	ELSE IF (natural(staind)) THEN
1498	num_grass_sinagec = num_grass_sinagec+1
1499	ELSE
1500	num_crop_sinagec = num_crop_sinagec+1
1501	ENDIF
1502
1503	ELSE
1504	IF (is_tree(staind)) THEN
1505	num_tree_mulagec = num_tree_mulagec+1
1506	ELSE IF (is_grassland_manag(staind)) THEN
1507	num_pasture_mulagec = num_pasture_mulagec+1
1508	ELSE IF (natural(staind)) THEN
1509	num_grass_mulagec = num_grass_mulagec+1
1510	ELSE
1511	num_crop_mulagec = num_crop_mulagec+1
1512	ENDIF
1513	ENDIF
1514	ENDDO
1515
1516	!! Allocate index array
1517	! allocate all index
1518	ALLOCATE(indall_tree(num_tree_sinagec+num_tree_mulagec*nagec_tree))
1519	ALLOCATE(indall_grass(num_grass_sinagec+num_grass_mulagec*nagec_herb))
1520	ALLOCATE(indall_pasture(num_pasture_sinagec+num_pasture_mulagec*nagec_herb))
1521	ALLOCATE(indall_crop(num_crop_sinagec+num_crop_mulagec*nagec_herb))
1522
1523	! allocate old-ageclass index
1524	ALLOCATE(indold_tree(num_tree_sinagec+num_tree_mulagec))
1525	ALLOCATE(indold_grass(num_grass_sinagec+num_grass_mulagec))
1526	ALLOCATE(indold_pasture(num_pasture_sinagec+num_pasture_mulagec))
1527	ALLOCATE(indold_crop(num_crop_sinagec+num_crop_mulagec))
1528
1529	!! 1.2 Fill the oldest-age-class and all index arrays
1530	itree=0
1531	igrass=0
1532	ipasture=0
1533	icrop=0
1534	itree2=1
1535	igrass2=1
1536	ipasture2=1
1537	icrop2=1
1538	DO ivma =2,nvmap
1539	staind=start_index(ivma)
1540	IF (is_tree(staind)) THEN
1541	itree=itree+1
1542	indold_tree(itree) = staind+nagec_pft(ivma)-1
1543	DO j = 0,nagec_pft(ivma)-1
1544	indall_tree(itree2+j) = staind+j
1545	ENDDO
1546	itree2=itree2+nagec_pft(ivma)
1547	ELSE IF (natural(staind) .AND. .NOT. is_grassland_manag(staind)) THEN
1548	igrass=igrass+1
1549	indold_grass(igrass) = staind+nagec_pft(ivma)-1
1550	DO j = 0,nagec_pft(ivma)-1
1551	indall_grass(igrass2+j) = staind+j
1552	ENDDO
1553	igrass2=igrass2+nagec_pft(ivma)
1554	ELSE IF (is_grassland_manag(staind)) THEN
1555	ipasture = ipasture+1
1556	indold_pasture(ipasture) = staind+nagec_pft(ivma)-1
1557	DO j = 0,nagec_pft(ivma)-1
1558	indall_pasture(ipasture2+j) = staind+j
1559	ENDDO
1560	ipasture2=ipasture2+nagec_pft(ivma)
1561	ELSE
1562	icrop = icrop+1
1563	indold_crop(icrop) = staind+nagec_pft(ivma)-1
1564	DO j = 0,nagec_pft(ivma)-1
1565	indall_crop(icrop2+j) = staind+j
1566	ENDDO
1567	icrop2=icrop2+nagec_pft(ivma)
1568	ENDIF
1569	ENDDO
1570
1571	!! 1.3 Allocate and fill other age class index
1572
1573	! [chaoyuejoy@gmail.com 2015-08-05]
1574	! note that we treat the case of (num_tree_mulagec==0) differently. In this
1575	! case there is no distinction of age groups among tree PFTs. But we still
1576	! we want to use the "gross_lcchange" subroutine. In this case we consider
1577	! them as having a single age group. In the subroutines
1578	! of "type_conversion" and "cross_give_receive", only the youngest-age-group
1579	! PFTs of a given MTC or vegetation type could receive the incoming fractions.
1580	! To be able to handle this case with least amount of code change, we assign the index
1581	! of PFT between youngest and second-oldes (i.e., indagec_tree etc) the same as
1582	! those of oldest tree PFTs (or all tree PFTs because in this cases these two indices
1583	! are identical) . So that this case could be correctly handled in the subrountines
1584	! of "type_conversion" and "cross_give_receive". This treatment allows use
1585	! of gross land cover change subroutine with only one single age class. This single
1586	! age class is "simultanously the oldest and youngest age class". At the same
1587	! time, we also change the num_tree_mulagec as the same of num_crop_sinagec.
1588	! The similar case also applies in grass,pasture and crop.
1589
1590	IF (num_tree_mulagec .EQ. 0) THEN
1591	ALLOCATE(indagec_tree(num_tree_sinagec,1))
1592	indagec_tree(:,1) = indall_tree(:)
1593	num_tree_mulagec = num_tree_sinagec
1594	ELSE
1595	ALLOCATE(indagec_tree(num_tree_mulagec,nagec_tree-1))
1596	END IF
1597
1598	IF (num_grass_mulagec .EQ. 0) THEN
1599	ALLOCATE(indagec_grass(num_grass_sinagec,1))
1600	indagec_grass(:,1) = indall_grass(:)
1601	num_grass_mulagec = num_grass_sinagec
1602	ELSE
1603	ALLOCATE(indagec_grass(num_grass_mulagec,nagec_herb-1))
1604	END IF
1605
1606	IF (num_pasture_mulagec .EQ. 0) THEN
1607	ALLOCATE(indagec_pasture(num_pasture_sinagec,1))
1608	indagec_pasture(:,1) = indall_pasture(:)
1609	num_pasture_mulagec = num_pasture_sinagec
1610	ELSE
1611	ALLOCATE(indagec_pasture(num_pasture_mulagec,nagec_herb-1))
1612	END IF
1613
1614	IF (num_crop_mulagec .EQ. 0) THEN
1615	ALLOCATE(indagec_crop(num_crop_sinagec,1))
1616	indagec_crop(:,1) = indall_crop(:)
1617	num_crop_mulagec = num_crop_sinagec
1618	ELSE
1619	ALLOCATE(indagec_crop(num_crop_mulagec,nagec_herb-1))
1620	END IF
1621
1622	! fill the non-oldest age class index arrays when number of age classes
1623	! is more than 1.
1624	! [chaoyuejoy@gmail.com, 2015-08-05]
1625	! Note the corresponding part of code will be automatically skipped
1626	! when nagec_tree ==1 and/or nagec_herb ==1, i.e., the assginment
1627	! in above codes when original num_*_mulagec variables are zero will be retained.
1628	itree=0
1629	igrass=0
1630	ipasture=0
1631	icrop=0
1632	DO ivma = 2,nvmap
1633	staind=start_index(ivma)
1634	IF (nagec_pft(ivma) > 1) THEN
1635	IF (is_tree(staind)) THEN
1636	itree=itree+1
1637	DO j = 1,nagec_tree-1
1638	indagec_tree(itree,j) = staind+nagec_tree-j-1
1639	ENDDO
1640	ELSE IF (natural(staind) .AND. .NOT. is_grassland_manag(staind)) THEN
1641	igrass=igrass+1
1642	DO j = 1,nagec_herb-1
1643	indagec_grass(igrass,j) = staind+nagec_herb-j-1
1644	ENDDO
1645	ELSE IF (is_grassland_manag(staind)) THEN
1646	ipasture=ipasture+1
1647	DO j = 1,nagec_herb-1
1648	indagec_pasture(ipasture,j) = staind+nagec_herb-j-1
1649	ENDDO
1650	ELSE
1651	icrop=icrop+1
1652	DO j = 1,nagec_herb-1
1653	indagec_crop(icrop,j) = staind+nagec_herb-j-1
1654	ENDDO
1655	ENDIF
1656	ENDIF
1657	ENDDO
1658
1659
1660	! we make copies of original input veget_max
1661	! veget_max will be modified through different operations in order to
1662	! check various purposes, e.g., whether input glcc is compatible with
1663	! existing veget_max and how to allocate it etc.
1664	! veget_max_old will not be modified
1665	veget_max(:,:) = veget_max_org(:,:)
1666	veget_max_old(:,:) = veget_max_org(:,:)
1667
1668	!! 2. Calcuate the fractions covered by tree, grass, pasture and crops
1669	!! for each age class
1670
1671	!************************************************************************!
1672	!**block to calculate fractions for basic veg types and age classes *!
1673	! Note:
1674	! 1. "calc_cover" subroutine does not depend on how many age classes
1675	! there are in each MTC.
1676	! 2. Fraction of baresoil is excluded here. This means transformation
1677	! of baresoil to a vegetated PFT is excluded in gross land cover change.
1678	veget_mtc(:,:) = 0.
1679	vegagec_tree(:,:) = 0.
1680	vegagec_grass(:,:) = 0.
1681	vegagec_pasture(:,:) = 0.
1682	vegagec_crop(:,:) = 0.
1683
1684
1685	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
1686	vegagec_pasture,vegagec_crop)
1687
1688	veget_tree(:) = SUM(vegagec_tree(:,:),DIM=2)
1689	veget_grass(:) = SUM(vegagec_grass(:,:),DIM=2)
1690	veget_pasture(:) = SUM(vegagec_pasture(:,:),DIM=2)
1691	veget_crop(:) = SUM(vegagec_crop(:,:),DIM=2)
1692	itree=1
1693	igrass=2
1694	ipasture=3
1695	icrop=4
1696	veget_4veg(:,itree) = veget_tree(:)
1697	veget_4veg(:,igrass) = veget_grass(:)
1698	veget_4veg(:,ipasture) = veget_pasture(:)
1699	veget_4veg(:,icrop) = veget_crop(:)
1700	!**end block to calculate fractions for basic veg types and age classes *!
1701	!****************************************************************************!
1702
1703	!******************** block to handle forestry harvest **************
1704	!! 2B. Here we handle the forestry wood harvest
1705	! Rules:
1706	! 1. We take first from second oldest forest, then oldest forest
1707
1708	pf2yf=1 !primary to young forest conversion because of harvest
1709	sf2yf=2 !old secondary to young forest conversion because of harvest
1710
1711	!! Note that Deficit_pf2yf and Deficit_sf2yf are temporary, intermediate
1712	!! variables. The final deficits after mutual compensation are stored in
1713	!! Deficit_pf2yf_final and Deficit_sf2yf_final.
1714	Deficit_pf2yf(:) = zero
1715	Deficit_sf2yf(:) = zero
1716	Deficit_pf2yf_final(:) = zero
1717	Deficit_sf2yf_final(:) = zero
1718
1719	!! Note that both Surplus_pf2yf and Surplus_sf2yf and temporary intermediate
1720	!! variables, the final surplus after mutual compensation are not outputed.
1721	Surplus_pf2yf(:) = zero
1722	Surplus_sf2yf(:) = zero
1723
1724	!! Note in the naming of pf2yf_compen_sf2yf and sf2yf_compen_pf2yf, active
1725	!! tense is used.
1726	pf2yf_compen_sf2yf(:) = zero !primary->young conversion that compensates
1727	!the secondary->young conversion because of deficit
1728	!in the latter
1729	sf2yf_compen_pf2yf(:) = zero !seondary->young conversion that compensates
1730	!the primary->young conversion because of the deficit
1731	!in the latter
1732
1733
1734	!! Define the "real" harvest matrix after considering the mutual compenstation
1735	!! between primary->young and secondary->young transitions.
1736	HmatrixReal(:,:) = zero !Harvest matrix real, used to hold the
1737	!harvest matrix after considering the mutual
1738	!compensation between primary and old secondary
1739	!forest
1740
1741	! we sum together harvest from primary and secondary forest and consider
1742	! as all happening on parimary forest.
1743	HmatrixReal(:,1) = harvest_matrix(:,pf2yf) + harvest_matrix(:,sf2yf)
1744
1745	! Check the availability of forest fractions for harvest
1746	WHERE (veget_tree(:) .LE. HmatrixReal(:,1))
1747	Deficit_pf2yf_final(:) = veget_tree(:)-HmatrixReal(:,1)
1748	HmatrixReal(:,1) = veget_tree(:)
1749	ENDWHERE
1750
1751
1752	glcc_pft(:,:) = 0.
1753	glcc_pft_tmp(:,:) = 0.
1754	glcc_pftmtc(:,:,:) = 0.
1755
1756	!! Allocate harvest-caused out-going primary and secondary forest fraction
1757	!! into different primary and secondary forest PFTs.
1758	! [Note: here we need only glcc_pft, but not glcc_pft_tmp and glcc_pftmtc.
1759	! The latter two variables will be set to zero again when handling LCC in
1760	! later sections.]
1761	DO ipts=1,npts
1762	!pf2yf
1763	CALL type_conversion(ipts,pf2yf,HmatrixReal,veget_mtc, &
1764	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1765	1,nagec_herb, &
1766	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1767	ENDDO
1768
1769	! Because we use the container of type_conversion, now the glcc_pft_tmp
1770	! and glcc_pftmtc have wrong information (because harvest loss is assigned
1771	! on the newly created youngest-age-class pasture/crop MTCs). So they have
1772	! to be re-initialized to zero. Only the information in glcc_pft is what
1773	! we need.
1774	glcc_pft_tmp(:,:) = 0.
1775	glcc_pftmtc(:,:,:) = 0.
1776	!Here we need to put glcc_pft into glcc_pftmtc for forestry harvest.
1777	!The same MTC will be maintained when forest is harvested.
1778	DO ivm =1,nvm
1779	IF (is_tree(ivm)) THEN
1780	glcc_pftmtc(:,ivm,pft_to_mtc(ivm)) = glcc_pft(:,ivm)
1781	ENDIF
1782	ENDDO
1783	!**************** end block to handle forestry harvest **************
1784	veget_max_tmp(:,:) = veget_max(:,:)
1785
1786
1787	!************************************************************************!
1788	!**block to calculate fractions for basic veg types and age classes *!
1789	! Note:
1790	! 1. "calc_cover" subroutine does not depend on how many age classes
1791	! there are in each MTC.
1792	! 2. Fraction of baresoil is excluded here. This means transformation
1793	! of baresoil to a vegetated PFT is excluded in gross land cover change.
1794	veget_mtc(:,:) = 0.
1795	vegagec_tree(:,:) = 0.
1796	vegagec_grass(:,:) = 0.
1797	vegagec_pasture(:,:) = 0.
1798	vegagec_crop(:,:) = 0.
1799
1800
1801	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
1802	vegagec_pasture,vegagec_crop)
1803
1804	veget_tree(:) = SUM(vegagec_tree(:,:),DIM=2)
1805	veget_grass(:) = SUM(vegagec_grass(:,:),DIM=2)
1806	veget_pasture(:) = SUM(vegagec_pasture(:,:),DIM=2)
1807	veget_crop(:) = SUM(vegagec_crop(:,:),DIM=2)
1808	itree=1
1809	igrass=2
1810	ipasture=3
1811	icrop=4
1812	veget_4veg(:,itree) = veget_tree(:)
1813	veget_4veg(:,igrass) = veget_grass(:)
1814	veget_4veg(:,ipasture) = veget_pasture(:)
1815	veget_4veg(:,icrop) = veget_crop(:)
1816	!**end block to calculate fractions for basic veg types and age classes *!
1817	!****************************************************************************!
1818
1819	!! 3. Decompose the LCC matrix to different PFTs
1820	!! We do this through several steps:
1821	! 3.1 Check whether input LCC matrix is feasible with current PFT fractions
1822	! (i.e., the fractions of forest,grass,pasture and crops)
1823	! and if not, adjust the transfer matrix by compensating the deficits
1824	! using the surpluses.
1825	! 3.2 Allocate the decreasing fractions of tree/grass/pasture/crop to their
1826	! respective age classes, in the sequences of old->young.
1827	! 3.3 Allocate the incoming fractions of tree/grass/pasture/crop to their
1828	! respective youngest age classes. The incoming fractions are distributed
1829	! according to the existing fractions of youngest-age-class PFTs of the
1830	! same receiving vegetation type. If none of them exists, the incoming
1831	! fraction is distributed equally.
1832
1833	!! 3.1 Adjust LCC matrix if it's not feasible with current PFT fractions
1834
1835	glcc(:,:) = glccSecondShift+glccPrimaryShift+glccNetLCC
1836	IncreDeficit(:,:) = 0.
1837	glccReal(:,:) = 0.
1838	glccDef(:,:) = 0.
1839
1840	!to crop - sequence: p2c,g2c,f2c
1841	CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
1842	p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
1843	IncreDeficit)
1844
1845	!to pasture - sequence: g2p,c2p,f2p
1846	CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
1847	g2p,igrass,c2p,icrop,f2p,itree,ipasture, &
1848	IncreDeficit)
1849
1850	!to grass - sequence: p2g,c2g,f2g
1851	CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
1852	p2g,ipasture,c2g,icrop,f2g,itree,igrass, &
1853	IncreDeficit)
1854
1855	!to forest - sequence: c2f,p2f,g2f
1856	CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
1857	c2f,icrop,p2f,ipasture,g2f,igrass,itree, &
1858	IncreDeficit)
1859
1860	!! 3.2 & 3.3 Allocate LCC matrix to different PFTs/age-classes
1861
1862	! because we use veget_max as a proxy variable and it has been changed
1863	! when we derive the glccReal, so here we have to recover its original
1864	! values, which is veget_max_tmp after the forestry harvest.
1865	veget_max(:,:) = veget_max_tmp(:,:)
1866
1867	! Calculate again fractions for different age-classes.
1868	veget_mtc(:,:) = 0.
1869	vegagec_tree(:,:) = 0.
1870	vegagec_grass(:,:) = 0.
1871	vegagec_pasture(:,:) = 0.
1872	vegagec_crop(:,:) = 0.
1873
1874	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
1875	vegagec_pasture,vegagec_crop)
1876
1877	! We allocate in the sequences of old->young. Within the same age-class
1878	! group, we allocate in proportion with existing PFT fractions.
1879	DO ipts=1,npts
1880	!f2c
1881	CALL type_conversion(ipts,f2c,glccReal,veget_mtc, &
1882	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1883	nagec_tree,nagec_herb, &
1884	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1885	!f2p
1886	CALL type_conversion(ipts,f2p,glccReal,veget_mtc, &
1887	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
1888	nagec_tree,nagec_herb, &
1889	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1890	!f2g
1891	CALL type_conversion(ipts,f2g,glccReal,veget_mtc, &
1892	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
1893	nagec_tree,nagec_herb, &
1894	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1895	!g2c
1896	CALL type_conversion(ipts,g2c,glccReal,veget_mtc, &
1897	indold_grass,indagec_grass,indagec_crop,num_crop_mulagec, &
1898	nagec_herb,nagec_herb, &
1899	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1900	!g2p
1901	CALL type_conversion(ipts,g2p,glccReal,veget_mtc, &
1902	indold_grass,indagec_grass,indagec_pasture,num_pasture_mulagec, &
1903	nagec_herb,nagec_herb, &
1904	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1905	!g2f
1906	CALL type_conversion(ipts,g2f,glccReal,veget_mtc, &
1907	indold_grass,indagec_grass,indagec_tree,num_tree_mulagec, &
1908	nagec_herb,nagec_tree, &
1909	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1910	!p2c
1911	CALL type_conversion(ipts,p2c,glccReal,veget_mtc, &
1912	indold_pasture,indagec_pasture,indagec_crop,num_crop_mulagec, &
1913	nagec_herb,nagec_herb, &
1914	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1915	!p2g
1916	CALL type_conversion(ipts,p2g,glccReal,veget_mtc, &
1917	indold_pasture,indagec_pasture,indagec_grass,num_grass_mulagec, &
1918	nagec_herb,nagec_herb, &
1919	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1920	!p2f
1921	CALL type_conversion(ipts,p2f,glccReal,veget_mtc, &
1922	indold_pasture,indagec_pasture,indagec_tree,num_tree_mulagec, &
1923	nagec_herb,nagec_tree, &
1924	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1925	!c2p
1926	CALL type_conversion(ipts,c2p,glccReal,veget_mtc, &
1927	indold_crop,indagec_crop,indagec_pasture,num_pasture_mulagec, &
1928	nagec_herb,nagec_herb, &
1929	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1930	!c2g
1931	CALL type_conversion(ipts,c2g,glccReal,veget_mtc, &
1932	indold_crop,indagec_crop,indagec_grass,num_grass_mulagec, &
1933	nagec_herb,nagec_herb, &
1934	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1935	!c2f
1936	CALL type_conversion(ipts,c2f,glccReal,veget_mtc, &
1937	indold_crop,indagec_crop,indagec_tree,num_tree_mulagec, &
1938	nagec_herb,nagec_tree, &
1939	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1940	ENDDO
1941
1942	END SUBROUTINE gross_glcc_firstday_SinAgeC_fh
1943
1944
1945	SUBROUTINE cross_give_receive(ipts,frac_used,veget_mtc, &
1946	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
1947	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1948
1949
1950	IMPLICIT NONE
1951
1952	!! 0. Input variables
1953	INTEGER, INTENT(in) :: ipts
1954	REAL(r_std), INTENT(in) :: frac_used !! fraction that the giving PFTs are going to collectively give
1955	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
1956	INTEGER, DIMENSION(:), INTENT(in) :: indold_tree !! Indices for PFTs giving out fractions;
1957	!! here use old tree cohort as an example
1958	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_crop !! Indices for secondary basic-vegetation cohorts; The youngest age classes
1959	!! of these vegetations are going to receive fractions.
1960	!! here we use crop cohorts as an example
1961	INTEGER, INTENT(in) :: num_crop_mulagec !! number of crop MTCs with more than one age classes
1962	INTEGER, INTENT(in) :: nagec_receive !! number of age classes in the receiving basic types
1963	!! (i.e., tree, grass, pasture, crop), here we can use crop
1964	!! as an example, nagec_receive=nagec_herb
1965
1966	!! 1. Modified variables
1967	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1968	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fractions each PFT is going to lose
1969	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
1970	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
1971	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft_tmp !! a temporary variable to hold the fractions each PFT is going to lose
1972
1973	!! Local vriables
1974	INTEGER :: j,ipft, iyoung
1975	REAL(r_std) :: totalveg
1976
1977
1978	! Out final objective is to know glcc_pftmtc, i.e., the fraction from each PFT
1979	! to the youngest age group of each MTC. We separate this task into two steps:
1980	! 1. we allocate the total outgoing fraction into the same age-class PFTs of
1981	! the a basic-vegetation (for example, the same age-calss PFTs of forest);
1982	! 2. we further allocate the outgoing fraction of each age-class PFT to
1983	! the different receiving youngest age-class PFTs of the same basic-vegetation
1984	! type, for example, the youngest age-calss PFTs of cropland.
1985
1986	! glcc_pft_tmp used only as a temporary variable to store the value
1987	glcc_pft_tmp(ipts,indold_tree) = veget_max(ipts,indold_tree)/SUM(veget_max(ipts,indold_tree))*frac_used
1988	glcc_pft(ipts,indold_tree) = glcc_pft(ipts,indold_tree) + glcc_pft_tmp(ipts,indold_tree)
1989	!we have to remove the outgoing fraction from veget_max in order to use this information for next loop
1990	veget_max(ipts,indold_tree) = veget_max(ipts,indold_tree) - glcc_pft_tmp(ipts,indold_tree)
1991
1992	! when receiving basic-vegetation type has a single age group, it will be considered as
1993	! both old and young age group (thus recevie the fraction donation), otherwise the youngest
1994	! age group is always the final element of indagec_crop.
1995	IF (nagec_receive == 1) THEN
1996	iyoung = 1
1997	ELSE
1998	iyoung = nagec_receive - 1
1999	ENDIF
2000
2001	totalveg = 0.
2002	DO j=1,num_crop_mulagec
2003	totalveg = totalveg + veget_mtc(ipts,agec_group(indagec_crop(j,iyoung)))
2004	ENDDO
2005
2006	IF (totalveg>min_stomate) THEN
2007	DO j=1,num_crop_mulagec
2008	ipft = indagec_crop(j,iyoung)
2009	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) = glcc_pft_tmp(ipts,indold_tree) &
2010	*veget_mtc(ipts,agec_group(ipft))/totalveg
2011	ENDDO
2012	ELSE
2013	DO j=1,num_crop_mulagec
2014	ipft = indagec_crop(j,iyoung)
2015	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) = glcc_pft_tmp(ipts,indold_tree)/num_crop_mulagec
2016	ENDDO
2017	ENDIF
2018
2019	END SUBROUTINE cross_give_receive
2020
2021	! ================================================================================================================================
2022	!! SUBROUTINE : type_conversion
2023	!>\BRIEF : Allocate outgoing into different age classes and incoming into
2024	!! yongest age-class of receiving MTCs.
2025	!!
2026	!! REMARK : The current dummy variables give an example of converting forests
2027	!! to crops.
2028	!! \n
2029	!_ ================================================================================================================================
2030	SUBROUTINE type_conversion(ipts,f2c,glccReal,veget_mtc, &
2031	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
2032	nagec_giving,nagec_receive, &
2033	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
2034	iagec_start)
2035
2036	IMPLICIT NONE
2037
2038	!! Input variables
2039	INTEGER, INTENT(in) :: ipts,f2c
2040	REAL(r_std), DIMENSION(:,:), INTENT(in) :: glccReal !! The "real" glcc matrix that we apply in the model
2041	!! after considering the consistency between presribed
2042	!! glcc matrix and existing vegetation fractions.
2043	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
2044	INTEGER, DIMENSION(:), INTENT(in) :: indold_tree !! Indices for PFTs giving out fractions;
2045	!! here use old tree cohort as an example
2046	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_tree !! Indices for PFTs giving out fractions;
2047	!! here use old tree cohort as an example
2048	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_crop !! Indices for secondary basic-vegetation cohorts; The youngest age classes
2049	!! of these vegetations are going to receive fractions.
2050	!! here we use crop cohorts as an example
2051	INTEGER, INTENT(in) :: num_crop_mulagec !! number of crop MTCs with more than one age classes
2052	INTEGER, INTENT(in) :: nagec_giving !! number of age classes in the giving basic types
2053	!! (i.e., tree, grass, pasture, crop), here we can use tree
2054	!! as an example, nagec=nagec_tree
2055	INTEGER, INTENT(in) :: nagec_receive !! number of age classes in the receiving basic types
2056	!! (i.e., tree, grass, pasture, crop), here we can use crop
2057	!! as an example, nagec=nagec_herb
2058	INTEGER, OPTIONAL, INTENT(in) :: iagec_start !! starting index for iagec, this is added in order to handle
2059	!! the case of secondary forest harvest.
2060
2061	!! 1. Modified variables
2062	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
2063	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
2064	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fractions each PFT is going to lose
2065	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
2066	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft_tmp !! Loss of fraction in each PFT
2067
2068	!! Local vriables
2069	INTEGER :: j,iagec,iagec_start_proxy
2070	REAL(r_std) :: frac_begin,frac_used
2071	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
2072	IF (.NOT. PRESENT(iagec_start)) THEN
2073	iagec_start_proxy=1
2074	ELSE
2075	iagec_start_proxy=iagec_start
2076	ENDIF
2077
2078	! This subroutine handles the conversion from one basic-vegetation type
2079	! to another, by calling the subroutine cross_give_receive, which handles
2080	! allocation of giving-receiving fraction among the giving age classes
2081	! and receiving basic-vegetation young age classes.
2082	! We allocate in the sequences of old->young. Within the same age-class
2083	! group, we allocate in proportion with existing PFT fractions. The same
2084	! also applies in the receiving youngest-age-class PFTs, i.e., the receiving
2085	! total fraction is allocated according to existing fractions of
2086	! MTCs of the same basic vegetation type, otherwise it will be equally
2087	! distributed.
2088
2089	frac_begin = glccReal(ipts,f2c)
2090	DO WHILE (frac_begin>min_stomate)
2091	DO iagec=iagec_start_proxy,nagec_giving
2092	IF (vegagec_tree(ipts,iagec)>frac_begin) THEN
2093	frac_used = frac_begin
2094	ELSE IF (vegagec_tree(ipts,iagec)>min_stomate) THEN
2095	frac_used = vegagec_tree(ipts,iagec)
2096	ELSE
2097	frac_used = 0.
2098	ENDIF
2099
2100	IF (frac_used>min_stomate) THEN
2101	IF (iagec==1) THEN
2102	! Note that vegagec_tree is fractions of tree age-class groups in the
2103	! the sequence of old->young, so iagec==1 means that we're handling
2104	! first the oldest-age-group tree PFTs.
2105	CALL cross_give_receive(ipts,frac_used,veget_mtc, &
2106	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
2107	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
2108	ELSE
2109	! Note also the sequence of indagec_tree is from old->young, so by
2110	! increasing iagec, we're handling progressively the old to young
2111	! tree age-class PFTs.
2112	CALL cross_give_receive(ipts,frac_used,veget_mtc, &
2113	indagec_tree(:,iagec-1),indagec_crop,nagec_receive,num_crop_mulagec, &
2114	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
2115	ENDIF
2116	frac_begin = frac_begin-frac_used
2117	vegagec_tree(ipts,iagec)=vegagec_tree(ipts,iagec)-frac_used
2118	ENDIF
2119	ENDDO
2120	ENDDO
2121
2122	END SUBROUTINE type_conversion
2123
2124	! ================================================================================================================================
2125	!! SUBROUTINE : calc_cover
2126	!!
2127	!>\BRIEF Calculate coverage fraction for different age classes of forest,
2128	!! grass, pasture and crops and also for each metaclass. Note baresoil is excluded.
2129	!!
2130	!! DESCRIPTION :
2131	!!
2132	!!
2133	!! MAIN OUTPUT VARIABLE(S) :
2134	!!
2135	!! \n
2136	!_ ================================================================================================================================
2137	SUBROUTINE calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
2138	vegagec_pasture,vegagec_crop)
2139
2140
2141	IMPLICIT NONE
2142
2143	!! Input variables
2144	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
2145	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
2146
2147	!! Output variables
2148	REAL(r_std), DIMENSION(npts,nvmap), INTENT(inout) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
2149	REAL(r_std), DIMENSION(npts,nagec_tree), INTENT(inout) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
2150	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_grass !! fraction of grass age-class groups, in sequence of old->young
2151	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_pasture !! fraction of pasture age-class groups, in sequence of old->young
2152	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_crop !! fraction of crop age-class groups, in sequence of old->young
2153
2154	!! Local variables
2155	INTEGER(i_std) :: ivma,staind,endind,j !! indices (unitless)
2156
2157	! Calculate veget_max for MTCs
2158	DO ivma = 1,nvmap
2159	staind = start_index(ivma)
2160	IF (nagec_pft(ivma) == 1) THEN
2161	veget_mtc(:,ivma) = veget_max(:,staind)
2162	ELSE
2163	veget_mtc(:,ivma) = \
2164	SUM(veget_max(:,staind:staind+nagec_pft(ivma)-1),DIM=2)
2165	ENDIF
2166	ENDDO
2167
2168	! Calculate veget_max for each age class
2169	DO ivma = 2,nvmap !here we start with 2 to exclude baresoil (always PFT1)
2170	staind = start_index(ivma)
2171	endind = staind+nagec_pft(ivma)-1
2172
2173	! Single-age-class MTC goest to oldest age class.
2174	IF (nagec_pft(ivma) == 1) THEN
2175	IF (is_tree(staind)) THEN
2176	vegagec_tree(:,1) = vegagec_tree(:,1)+veget_max(:,staind)
2177	ELSE IF (is_grassland_manag(staind)) THEN
2178	vegagec_pasture(:,1) = vegagec_pasture(:,1)+veget_max(:,staind)
2179	ELSE IF (natural(staind)) THEN
2180	vegagec_grass(:,1) = vegagec_grass(:,1)+veget_max(:,staind)
2181	ELSE
2182	vegagec_crop(:,1) = vegagec_crop(:,1)+veget_max(:,staind)
2183	ENDIF
2184
2185	ELSE
2186	IF (is_tree(staind)) THEN
2187	DO j=1,nagec_tree
2188	vegagec_tree(:,j) = vegagec_tree(:,j)+veget_max(:,endind-j+1)
2189	ENDDO
2190	ELSE IF (is_grassland_manag(staind)) THEN
2191	DO j=1,nagec_herb
2192	vegagec_pasture(:,j) = vegagec_pasture(:,j)+veget_max(:,endind-j+1)
2193	ENDDO
2194	ELSE IF (natural(staind)) THEN
2195	DO j=1,nagec_herb
2196	vegagec_grass(:,j) = vegagec_grass(:,j)+veget_max(:,endind-j+1)
2197	ENDDO
2198	ELSE
2199	DO j=1,nagec_herb
2200	vegagec_crop(:,j) = vegagec_crop(:,j)+veget_max(:,endind-j+1)
2201	ENDDO
2202	ENDIF
2203	ENDIF
2204	ENDDO
2205
2206	END SUBROUTINE calc_cover
2207
2208	! Note this subroutine does not depend on how many age classes there are
2209	! in different MTCs.
2210	SUBROUTINE glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
2211	p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
2212	IncreDeficit)
2213
2214	IMPLICIT NONE
2215
2216	!! 0.1 Input variables
2217	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
2218	INTEGER, INTENT(in) :: p2c,ipasture,g2c,igrass,f2c,itree,icrop
2219	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glcc !! the land-cover-change (LCC) matrix in case a gross LCC is
2220	!! used.
2221
2222	!! 0.2 Output variables
2223
2224
2225	!! 0.3 Modified variables
2226	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: veget_4veg !! "maximal" coverage of tree/grass/pasture/crop
2227	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccDef !! Gross LCC deficit, negative values mean that there
2228	!! are not enough fractions in the source vegetations
2229	!! to the target ones as presribed by the LCC matrix.
2230	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
2231	!! after considering the consistency between presribed
2232	!! glcc matrix and existing vegetation fractions.
2233	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
2234	!! there are not enough fractions in the source PFTs
2235	!! /vegetations to target PFTs/vegetations. I.e., these
2236	!! fraction transfers are presribed in LCC matrix but
2237	!! not realized.
2238
2239	!! 0.4 Local variables
2240	REAL(r_std), DIMENSION(npts) :: tmpdef !! LCC deficits by summing up all the deficits to the
2241	!! the same target vegetation.
2242
2243
2244	!! 0. We first handle the cases where veget_4veg might be very small
2245	!tree
2246	WHERE(veget_4veg(:,itree) > min_stomate)
2247	glccDef(:,f2c) = veget_4veg(:,itree)-glcc(:,f2c)
2248	WHERE(veget_4veg(:,itree)>glcc(:,f2c))
2249	glccReal(:,f2c) = glcc(:,f2c)
2250	ELSEWHERE
2251	glccReal(:,f2c) = veget_4veg(:,itree)
2252	ENDWHERE
2253	ELSEWHERE
2254	glccReal(:,f2c) = 0.
2255	glccDef(:,f2c) = -1*glcc(:,f2c)
2256	ENDWHERE
2257
2258	!pasture
2259	WHERE(veget_4veg(:,ipasture) > min_stomate)
2260	glccDef(:,p2c) = veget_4veg(:,ipasture)-glcc(:,p2c)
2261	WHERE(veget_4veg(:,ipasture)>glcc(:,p2c))
2262	glccReal(:,p2c) = glcc(:,p2c)
2263	ELSEWHERE
2264	glccReal(:,p2c) = veget_4veg(:,ipasture)
2265	ENDWHERE
2266	ELSEWHERE
2267	glccReal(:,p2c) = 0.
2268	glccDef(:,p2c) = -1*glcc(:,p2c)
2269	ENDWHERE
2270
2271	!grass
2272	WHERE(veget_4veg(:,igrass) > min_stomate)
2273	glccDef(:,g2c) = veget_4veg(:,igrass)-glcc(:,g2c)
2274	WHERE(veget_4veg(:,igrass)>glcc(:,g2c))
2275	glccReal(:,g2c) = glcc(:,g2c)
2276	ELSEWHERE
2277	glccReal(:,g2c) = veget_4veg(:,igrass)
2278	ENDWHERE
2279	ELSEWHERE
2280	glccReal(:,g2c) = 0.
2281	glccDef(:,g2c) = -1*glcc(:,g2c)
2282	ENDWHERE
2283
2284	!! 1. Compensation sequence: pasture,grass,forest
2285	tmpdef(:) = glccDef(:,f2c)+glccDef(:,g2c)+glccDef(:,p2c)
2286	WHERE(glccDef(:,p2c)<0)
2287	WHERE(glccDef(:,g2c)<0)
2288	WHERE(glccDef(:,f2c)<0) ! 1 (-,-,-)
2289	IncreDeficit(:,icrop) = tmpdef(:)
2290	ELSEWHERE ! 2 (-,-,+)
2291	WHERE(tmpdef(:)>=min_stomate)
2292	glccReal(:,f2c) = glccReal(:,f2c)-glccDef(:,g2c)-glccDef(:,p2c)
2293	ELSEWHERE
2294	glccReal(:,f2c) = veget_4veg(:,itree)
2295	IncreDeficit(:,icrop) = tmpdef(:)
2296	ENDWHERE
2297	ENDWHERE
2298	ELSEWHERE
2299	WHERE(glccDef(:,f2c)<0) ! 3 (-,+,-)
2300	WHERE(tmpdef(:)>=min_stomate)
2301	glccReal(:,g2c) = glccReal(:,g2c)-glccDef(:,p2c)-glccDef(:,f2c)
2302	ELSEWHERE
2303	glccReal(:,g2c) = veget_4veg(:,igrass)
2304	IncreDeficit(:,icrop) = tmpdef(:)
2305	ENDWHERE
2306	ELSEWHERE ! 4 (-,+,+)
2307	WHERE(tmpdef(:)>=min_stomate)
2308	WHERE((glccDef(:,g2c)+glccDef(:,p2c))>=min_stomate)
2309	glccReal(:,g2c) = glccReal(:,g2c)-glccDef(:,p2c)
2310	ELSEWHERE
2311	glccReal(:,g2c) = veget_4veg(:,igrass)
2312	glccReal(:,f2c) = glccReal(:,f2c)-(glccDef(:,p2c)+glccDef(:,g2c))
2313	ENDWHERE
2314	ELSEWHERE
2315	glccReal(:,g2c) = veget_4veg(:,igrass)
2316	glccReal(:,f2c) = veget_4veg(:,itree)
2317	IncreDeficit(:,icrop) = tmpdef(:)
2318	ENDWHERE
2319	ENDWHERE
2320	ENDWHERE
2321	ELSEWHERE
2322	WHERE(glccDef(:,g2c)<0)
2323	WHERE(glccDef(:,f2c)<0) ! 5 (+,-,-)
2324	WHERE(tmpdef(:)>=min_stomate)
2325	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,g2c)-glccDef(:,f2c)
2326	ELSEWHERE
2327	IncreDeficit(:,icrop) = tmpdef(:)
2328	glccReal(:,p2c) = veget_4veg(:,ipasture)
2329	ENDWHERE
2330	ELSEWHERE ! 6 (+,-,+)
2331	WHERE(tmpdef(:)>=min_stomate)
2332	WHERE((glccDef(:,p2c)+glccDef(:,g2c))>=min_stomate)
2333	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,g2c)
2334	ELSEWHERE
2335	glccReal(:,p2c) = veget_4veg(:,ipasture)
2336	glccReal(:,f2c) = glccReal(:,f2c)-(glccDef(:,g2c)+glccDef(:,p2c))
2337	ENDWHERE
2338	ELSEWHERE
2339	IncreDeficit(:,icrop) = tmpdef(:)
2340	glccReal(:,p2c) = veget_4veg(:,ipasture)
2341	glccReal(:,f2c) = veget_4veg(:,itree)
2342	ENDWHERE
2343	ENDWHERE
2344	ELSEWHERE
2345	WHERE(glccDef(:,f2c)<0) ! 7 (+,+,-)
2346	WHERE(tmpdef(:)>=min_stomate)
2347	WHERE((glccDef(:,p2c)+glccDef(:,f2c))>=min_stomate)
2348	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,f2c)
2349	ELSEWHERE
2350	glccReal(:,p2c) = veget_4veg(:,ipasture)
2351	glccReal(:,g2c) = glccReal(:,g2c)-(glccDef(:,f2c)+glccDef(:,p2c))
2352	ENDWHERE
2353	ELSEWHERE
2354	IncreDeficit(:,icrop) = tmpdef(:)
2355	glccReal(:,g2c) = veget_4veg(:,igrass)
2356	glccReal(:,p2c) = veget_4veg(:,ipasture)
2357	ENDWHERE
2358	ELSEWHERE ! 8 (+,+,+)
2359	!do nothing
2360	ENDWHERE
2361	ENDWHERE
2362	ENDWHERE
2363	veget_4veg(:,itree) = veget_4veg(:,itree) - glccReal(:,f2c)
2364	veget_4veg(:,igrass) = veget_4veg(:,igrass) - glccReal(:,g2c)
2365	veget_4veg(:,ipasture) = veget_4veg(:,ipasture) - glccReal(:,p2c)
2366
2367	END SUBROUTINE glcc_compensation_full
2368
2369
2370
2371	!! This subroutine implements non-full compensation, is currently
2372	!! abandoned.
2373	SUBROUTINE glcc_compensation(npts,veget_4veg,glcc,glccDef, &
2374	p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
2375	IncreDeficit)
2376
2377	IMPLICIT NONE
2378
2379	!! 0.1 Input variables
2380	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
2381	REAL(r_std), DIMENSION(npts,4), INTENT(in) :: veget_4veg !! "maximal" coverage fraction of a PFT on the ground
2382	INTEGER, INTENT(in) :: p2c,ipasture,g2c,igrass,f2c,itree,icrop
2383
2384	!! 0.2 Output variables
2385
2386
2387	!! 0.3 Modified variables
2388	REAL(r_std), DIMENSION (npts,12),INTENT(inout) :: glcc !! the land-cover-change (LCC) matrix in case a gross LCC is
2389	!! used.
2390	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccDef !! Gross LCC deficit, negative values mean that there
2391	!! are not enough fractions in the source vegetations
2392	!! to the target ones as presribed by the LCC matrix.
2393	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
2394	!! there are not enough fractions in the source PFTs
2395	!! /vegetations to target PFTs/vegetations. I.e., these
2396	!! fraction transfers are presribed in LCC matrix but
2397	!! not realized.
2398
2399	!! 0.4 Local variables
2400	REAL(r_std), DIMENSION(npts) :: glccDef_all !! LCC deficits by summing up all the deficits to the
2401	!! the same target vegetation.
2402
2403
2404	WHERE(veget_4veg(:,itree) > min_stomate)
2405	glccDef(:,f2c) = veget_4veg(:,itree)-glcc(:,f2c)
2406	ELSEWHERE
2407	glccDef(:,f2c) = -1*glcc(:,f2c)
2408	glcc(:,f2c) = 0.
2409	ENDWHERE
2410
2411	WHERE(veget_4veg(:,ipasture) > min_stomate)
2412	glccDef(:,p2c) = veget_4veg(:,ipasture)-glcc(:,p2c)
2413	ELSEWHERE
2414	glccDef(:,p2c) = -1*glcc(:,p2c)
2415	glcc(:,p2c) = 0.
2416	ENDWHERE
2417
2418	WHERE(veget_4veg(:,igrass) > min_stomate)
2419	glccDef(:,g2c) = veget_4veg(:,igrass)-glcc(:,g2c)
2420	ELSEWHERE
2421	glccDef(:,g2c) = -1*glcc(:,g2c)
2422	glcc(:,g2c) = 0.
2423	ENDWHERE
2424
2425	glccDef_all(:) = glccDef(:,f2c)+glccDef(:,p2c)+glccDef(:,g2c)
2426
2427	! We allow the surpluses/deficits in p2c and g2c mutually compensating
2428	! for each other. If there are still deficits after this compensation,
2429	! they will be further compensated for by the surpluses from f2c (if there are any
2430	! surpluses). The ultimate deficits that cannot be compensated for
2431	! will be recorded and dropped.
2432
2433	! Because we assume the "pasture rule" is used, i.e., the crops
2434	! are supposed to come primarily from pastures and grasses, normally
2435	! we expect the deficits to occur in p2c or g2c rather than in f2c. But
2436	! if it happens that f2c has deficits while p2c or g2c has surpluse,
2437	! the surpluses will not be used to compensate for the f2c-deficits,
2438	! instead, we will just record and drop the f2c-deficits.
2439
2440	! In following codes for convenience we're not going to check
2441	! whether surpluses in f2c are enough to compensate for deficits
2442	! in p2c or g2c or both. Instead, we just add their deficits on top
2443	! of f2c. The issues of not-enough surpluses in f2c will be left for
2444	! the codes after this section to handle.
2445	WHERE (glccDef(:,p2c) < 0.)
2446	glcc(:,p2c) = veget_4veg(:,ipasture)
2447	WHERE (glccDef(:,g2c) < 0.)
2448	glcc(:,g2c) = veget_4veg(:,igrass)
2449	ELSEWHERE
2450	WHERE (glccDef(:,g2c)+glccDef(:,p2c) > min_stomate)
2451	glcc(:,g2c) = glcc(:,g2c)-glccDef(:,p2c)
2452	ELSEWHERE
2453	glcc(:,g2c) = veget_4veg(:,igrass)
2454	! whatever the case, we simply add the dificts to f2c
2455	glcc(:,f2c) = glcc(:,f2c)-glccDef(:,p2c)-glccDef(:,g2c)
2456	ENDWHERE
2457	ENDWHERE
2458
2459	ELSEWHERE
2460	WHERE(glccDef(:,g2c) < 0.)
2461	glcc(:,g2c) = veget_4veg(:,igrass)
2462	WHERE(glccDef(:,p2c)+glccDef(:,g2c) > min_stomate)
2463	glcc(:,p2c) = glcc(:,p2c)-glccDef(:,g2c)
2464	ELSEWHERE
2465	glcc(:,p2c) = veget_4veg(:,ipasture)
2466	! whatever the case, we simply add the dificts to f2c
2467	glcc(:,f2c) = glcc(:,f2c)-glccDef(:,p2c)-glccDef(:,g2c)
2468	ENDWHERE
2469	ELSEWHERE
2470	!Here p2c and g2c both show surplus, we're not going to check whether
2471	!glccDef(:,f2c) has negative values because we assume a "pasture rule"
2472	!is applied when constructing the gross LCC matrix, so deficits in
2473	!f2c will just be dropped but not be compensated for by the surpluses in
2474	!p2c or g2c.
2475	ENDWHERE
2476	ENDWHERE
2477
2478	! 1. We calculate again the f2c-deficit because f2c-glcc is adjusted in the
2479	! codes above as we allocated the deficits of p2c and g2c into f2c.
2480	! In cases where glccDef_all is less than zero, f2c-glcc will be larger
2481	! than available forest veget_max and we therefore limit the f2c-glcc to
2482	! available forest cover.
2483	! 2. There is (probably) a second case where glccDef_all is larger then zero,
2484	! but f2c-glcc is higher than veget_tree, i.e., Originally f2c is given a
2485	! high value that there is deficit in f2c but surpluses exist for p2c and g2c.
2486	! Normally we
2487	! assume this won't happen as explained above, given that a "pasture rule" was
2488	! used in constructing the gross LCC matrix. Nevertheless if this deos
2489	! happen, we will just drop the f2c deficit without being compensated
2490	! for by the surplus in p2c or g2c.
2491
2492	! we handle the 2nd case first
2493	WHERE(veget_4veg(:,itree) > min_stomate )
2494	WHERE(glccDef(:,f2c) < 0.)
2495	glcc(:,f2c) = veget_4veg(:,itree)
2496	WHERE (glccDef(:,p2c)+glccDef(:,g2c) > min_stomate)
2497	IncreDeficit(:,icrop) = glccDef(:,f2c)
2498	ELSEWHERE
2499	IncreDeficit(:,icrop) = glccDef_all(:)
2500	ENDWHERE
2501	ELSEWHERE
2502	WHERE(glccDef_all(:) < 0.) !handle the 1st case
2503	glcc(:,f2c) = veget_4veg(:,itree)
2504	IncreDeficit(:,icrop) = glccDef_all(:)
2505	ENDWHERE
2506	ENDWHERE
2507	ELSEWHERE
2508	WHERE(glccDef(:,p2c)+glccDef(:,g2c)>min_stomate)
2509	IncreDeficit(:,icrop) = glccDef(:,f2c)
2510	ELSEWHERE
2511	IncreDeficit(:,icrop) = glccDef_all(:)
2512	ENDWHERE
2513	ENDWHERE
2514
2515	END SUBROUTINE glcc_compensation
2516
2517
2518
2519	END MODULE stomate_glcchange_SinAgeC_fh

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_mict_peat_ch4/src_stomate/stomate_glcchange_SinAgeC_fh.f90 @ 7346

Download in other formats: