Context Navigation

stomate_glcchange_fh.f90 @ 5491

Last change on this file since 5491 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: 209.4 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_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_fh, gross_glcchange_fh, age_class_distr
40
41	CONTAINS
42
43	! ================================================================================================================================
44	!! SUBROUTINE : age_class_distr
45	!!
46	!>\BRIEF Redistribute biomass, litter, soilcarbon and water across
47	!! the age classes
48	!!
49	!! DESCRIPTION : Following growth, the trees from an age class may have become
50	!! too big to belong to this age class. The biomass, litter, soilcarbon and
51	!! soil water then need to be moved from one age class to the next age class.
52	!!
53	!! RECENT CHANGE(S) :
54	!!
55	!! MAIN OUTPUT VARIABLE(S) :
56	!!
57	!! REFERENCES : None
58	!!
59	!! FLOWCHART :
60	!! \n
61	!_ ================================================================================================================================
62
63	SUBROUTINE age_class_distr(npts, lalo, resolution, bound_spa, &
64	biomass, veget_max, ind, &
65	lm_lastyearmax, leaf_frac, co2_to_bm, &
66	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
67	everywhere, litter, carbon, lignin_struc, &
68	deepC_a, deepC_s, deepC_p, &
69	bm_to_litter, PFTpresent, when_growthinit,&
70	senescence, npp_longterm, gpp_daily, leaf_age, &
71	gdd_from_growthinit, gdd_midwinter, time_hum_min, gdd_m5_dormance, &
72	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
73	gpp_week, resp_maint, resp_growth, npp_daily)
74
75	IMPLICIT NONE
76
77	!! 0. Variable and parameter declaration
78
79	!! 0.1 Input variables
80
81	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
82	REAL(r_std),DIMENSION(npts,2),INTENT(in) :: lalo !! Geographical coordinates (latitude,longitude)
83	!! for pixels (degrees)
84	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution !! Resolution at each grid point (m)
85	!! [1=E-W, 2=N-S]
86
87	!! 0.2 Output variables
88
89
90	!! 0.3 Modified variables
91
92	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
93	!! each pixel
94	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
95	!! for deciduous trees)
96	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
97	!! unitless)
98	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
99	!! (0 to 1, unitless)
100	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
101	!! @tex $(gC m^{-2} day^{-1})$ @endtex
102	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
103	!! -5 deg C (for phenology)
104	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
105	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
106	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
107	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
108	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
109	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
110	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
111	!! the growing season (days)
112	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
113	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
114	!! @tex $(m^{-2})$ @endtex
115	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
116	!! May sum to
117	!! less than unity if the pixel has
118	!! nobio area. (unitless, 0-1)
119	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
120	!! @tex ($gC m^{-2}$) @endtex
121	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
122	!! very localized (after its introduction) (?)
123	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
124	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
125	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
126	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
127	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
128	!! availability (days)
129	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
130	!! (for phenology) - this is written to the
131	!! history files
132	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
133	!! for crops
134	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
135	!! C (for phenology)
136	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
137	!! leaves were lost (for phenology)
138	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
139	!! above and below ground
140	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
141	!! @tex ($gC m^{-2}$) @endtex
142	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
143	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
144	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
145	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
146	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
147	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
148	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
149	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
150	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
151	!! below ground @tex ($gC m^{-2}$) @endtex
152	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
153	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
154	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
155	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
156	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bound_spa !! Spatial age class boundaries.
157
158	!! 0.4 Local variables
159
160	INTEGER(i_std) :: ipts,ivm,igroup !! Indeces(unitless)
161	INTEGER(i_std) :: iele,ipar,ipft !! Indeces(unitless)
162	INTEGER(i_std) :: iagec,imbc,icirc !! Indeces(unitless)
163	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
164	INTEGER(i_std) :: ivma !! Indeces(unitless)
165	REAL(r_std) :: share_expanded !! Share of the veget_max of the existing vegetation
166	!! within a PFT over the total veget_max following
167	!! expansion of that PFT (unitless, 0-1)
168	!! @tex $(ind m^{-2})$ @endtex
169	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
170	!! mass balance chech for this routine
171	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
172	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
173	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
174	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start and end pool of this routine
175	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
176	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! Start and end pool of this routine
177	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
178	REAL(r_std), DIMENSION(nelements) :: temp_start !! Start and end pool of this routine
179	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
180	REAL(r_std), DIMENSION(nelements) :: temp_end !! Start and end pool of this routine
181	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
182	REAL(r_std), DIMENSION(nlitt,nlevs) :: litter_weight_expanded !! The fraction of litter on the expanded
183	!! PFT.
184	!! @tex $-$ @endtex
185	REAL(r_std), DIMENSION(npts,nvm) :: woodmass !! Woodmass of individuals (gC)
186	REAL(r_std), DIMENSION(npts,nvm) :: soilcarbon !!
187	REAL(r_std), DIMENSION(npts,nvm) :: agec_indicator !!
188	CHARACTER(LEN=80) :: data_filename
189
190	!_ ================================================================================================================================
191
192	IF (printlev.GE.3) WRITE(numout,*) 'Entering age class distribution'
193
194	!CALL getin_p('AgeC_Threshold_File',data_filename)
195	!CALL slowproc_read_data(npts, lalo, resolution, bound_spa, data_filename, 'matrix')
196
197	IF (.NOT. use_bound_spa) THEN
198	DO ipts = 1,npts
199	bound_spa(ipts,:) = age_class_bound(:)
200	ENDDO
201	ENDIF
202
203	!! 1. Initialize
204
205	woodmass(:,:) = biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
206	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)
207	soilcarbon(:,:) = -1 *(SUM(carbon(:,:,:),DIM=2) + &
208	SUM(SUM(litter(:,:,:,:,icarbon),DIM=2),DIM=3))
209
210	!! 1.2 Initialize check for mass balance closure
211	! The mass balance is calculated at the end of this routine
212	! in section 3. Initial biomass and harvest pool all other
213	! relevant pools were just set to zero.
214	pool_start(:,:,:) = zero
215	DO iele = 1,nelements
216
217	! co2_to_bm
218	pool_start(:,:,iele) = pool_start(:,:,iele) + co2_to_bm(:,:)
219
220	! Biomass pool + bm_to_litter
221	DO ipar = 1,nparts
222	pool_start(:,:,iele) = pool_start(:,:,iele) + &
223	(biomass(:,:,ipar,iele) + bm_to_litter(:,:,ipar,iele)) * &
224	veget_max(:,:)
225	ENDDO
226
227	! Litter pool (gC m-2) * (m2 m-2)
228	DO ilit = 1,nlitt
229	DO ilev = 1,nlevs
230	pool_start(:,:,iele) = pool_start(:,:,iele) + &
231	litter(:,ilit,:,ilev,iele) * veget_max(:,:)
232	ENDDO
233	ENDDO
234
235	! Soil carbon (gC m-2) * (m2 m-2)
236	DO icarb = 1,ncarb
237	pool_start(:,:,iele) = pool_start(:,:,iele) + &
238	carbon(:,icarb,:) * veget_max(:,:)
239	ENDDO
240
241	ENDDO
242
243
244	!! 2. Handle the merge of PFTs when one age class moves to the next one.
245
246	! Following growth, the value of age-class indicator variable
247	! from an age class may have become too big to stay
248	! in this age class. The biomass, litter, soilcarbon and soil
249	! water then need to be moved from one age class to the next age class.
250	DO ipts = 1,npts
251	! This loops over all the MTCs that we have ignoring age classes
252	DO ivma=1,nvmap
253	ivm=start_index(ivma)
254
255	! If we only have a single age class for this
256	! PFT, we can skip it.
257	IF(nagec_pft(ivma) .EQ. 1)CYCLE
258
259	IF(is_tree(ivm)) THEN
260	agec_indicator(:,:) = woodmass(:,:)
261	ELSE
262	agec_indicator(:,:) = soilcarbon(:,:)
263	ENDIF ! is_tree(ivm)
264
265	CALL check_merge_same_MTC(ipts, ivma, agec_indicator, bound_spa, &
266	biomass, veget_max, ind, &
267	lm_lastyearmax, leaf_frac, co2_to_bm, &
268	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
269	everywhere, litter, carbon, lignin_struc, &
270	deepC_a, deepC_s, deepC_p, &
271	bm_to_litter, PFTpresent, when_growthinit,&
272	senescence, npp_longterm, gpp_daily, leaf_age, &
273	gdd_from_growthinit, gdd_midwinter, time_hum_min, gdd_m5_dormance, &
274	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
275	gpp_week, resp_maint, resp_growth, npp_daily)
276
277	ENDDO ! Looping over MTCs
278	ENDDO ! loop over #pixels - domain size
279
280
281	!! 3. Mass balance closure
282
283	!! 3.1 Calculate components of the mass balance
284	pool_end(:,:,:) = zero
285	DO iele = 1,nelements
286
287	! co2_to_bm
288	pool_end(:,:,iele) = pool_end(:,:,iele) + co2_to_bm(:,:)
289
290	! Biomass pool + bm_to_litter
291	DO ipar = 1,nparts
292	pool_end(:,:,iele) = pool_end(:,:,iele) + &
293	(biomass(:,:,ipar,iele) + bm_to_litter(:,:,ipar,iele)) * &
294	veget_max(:,:)
295	ENDDO
296
297	! Litter pool (gC m-2) * (m2 m-2)
298	DO ilit = 1,nlitt
299	DO ilev = 1,nlevs
300	pool_end(:,:,iele) = pool_end(:,:,iele) + &
301	litter(:,ilit,:,ilev,iele) * veget_max(:,:)
302	ENDDO
303	ENDDO
304
305	! Soil carbon (gC m-2) * (m2 m-2)
306	DO icarb = 1,ncarb
307	pool_end(:,:,iele) = pool_end(:,:,iele) + &
308	carbon(:,icarb,:) * veget_max(:,:)
309	ENDDO
310	ENDDO
311
312	!! 3.2 Calculate mass balance
313	check_intern(:,:,iatm2land,icarbon) = zero
314	check_intern(:,:,iland2atm,icarbon) = -un * zero
315	check_intern(:,:,ilat2out,icarbon) = zero
316	check_intern(:,:,ilat2in,icarbon) = -un * zero
317	check_intern(:,:,ipoolchange,icarbon) = -un * (pool_end(:,:,icarbon) - pool_start(:,:,icarbon))
318	closure_intern = zero
319	DO imbc = 1,nmbcomp
320	closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + check_intern(:,:,imbc,icarbon)
321	ENDDO
322
323	!! 3.3 Write outcome of the check
324	! Sum over ivm because of age class redistribution
325	DO ipts = 1,npts
326	IF (SUM(closure_intern(ipts,:,icarbon)) .LT. min_stomate .AND. &
327	SUM(closure_intern(ipts,:,icarbon)) .GT. -min_stomate) THEN
328	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure: age_class_distr', ipts
329	ELSE
330	WRITE(numout,*) 'Error: mass balance is not closed in age_class_distr'
331	WRITE(numout,*) ' Difference, ipts, ', ipts, SUM(closure_intern(ipts,:,icarbon))
332	ENDIF
333	ENDDO
334
335	IF (printlev.GE.4) WRITE(numout,*) 'Leaving age class distribution'
336
337	END SUBROUTINE age_class_distr
338
339
340	SUBROUTINE check_merge_same_MTC(ipts, ivma, woodmass, bound_spa, &
341	biomass, veget_max, ind, &
342	lm_lastyearmax, leaf_frac, co2_to_bm, &
343	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
344	everywhere, litter, carbon, lignin_struc, &
345	deepC_a, deepC_s, deepC_p, &
346	bm_to_litter, PFTpresent, when_growthinit,&
347	senescence, npp_longterm, gpp_daily, leaf_age, &
348	gdd_from_growthinit, gdd_midwinter, time_hum_min, gdd_m5_dormance, &
349	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
350	gpp_week, resp_maint, resp_growth, npp_daily)
351
352	IMPLICIT NONE
353
354	!! 0. Variable and parameter declaration
355
356	!! 0.1 Input variables
357
358	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
359	INTEGER, INTENT(in) :: ivma !!
360	REAL(r_std), DIMENSION(:,:), INTENT(in) :: woodmass !! Woodmass of individuals (gC)
361	REAL(r_std), DIMENSION(:,:), INTENT(in) :: bound_spa !!
362
363	!! 0.2 Output variables
364
365
366	!! 0.3 Modified variables
367
368	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
369	!! each pixel
370	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
371	!! for deciduous trees)
372	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
373	!! unitless)
374	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
375	!! (0 to 1, unitless)
376	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
377	!! @tex $(gC m^{-2} day^{-1})$ @endtex
378	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
379	!! -5 deg C (for phenology)
380	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
381	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
382	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
383	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
384	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
385	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
386	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
387	!! the growing season (days)
388	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
389	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
390	!! @tex $(m^{-2})$ @endtex
391	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
392	!! May sum to
393	!! less than unity if the pixel has
394	!! nobio area. (unitless, 0-1)
395	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
396	!! @tex ($gC m^{-2}$) @endtex
397	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
398	!! very localized (after its introduction) (?)
399	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
400	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
401	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
402	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
403	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
404	!! availability (days)
405	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
406	!! (for phenology) - this is written to the
407	!! history files
408	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
409	!! for crops
410	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
411	!! C (for phenology)
412	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
413	!! leaves were lost (for phenology)
414	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
415	!! above and below ground
416	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
417	!! @tex ($gC m^{-2}$) @endtex
418	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
419	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
420	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
421	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
422	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
423	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
424	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
425	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
426	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
427	!! below ground @tex ($gC m^{-2}$) @endtex
428	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
429	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
430	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
431	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
432
433	!! 0.4 Local variables
434
435	INTEGER(i_std) :: iele,ipar,ipft !! Indeces(unitless)
436	INTEGER(i_std) :: iagec,imbc,icirc !! Indeces(unitless)
437	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
438	REAL(r_std) :: share_expanded !! Share of the veget_max of the existing vegetation
439	!! within a PFT over the total veget_max following
440	!! expansion of that PFT (unitless, 0-1)
441	!! @tex $(ind m^{-2})$ @endtex
442	REAL(r_std), DIMENSION(nlitt,nlevs) :: litter_weight_expanded !! The fraction of litter on the expanded
443	!! PFT.
444
445
446	!_ ================================================================================================================================
447
448	!! 1 Check if the trees still belong to this age class
449	! Note that the term age class is used but that the classes used in the
450	! code are not defined on an age criterion. Instead the biomass or
451	! or soil carbon pool is used.
452	IF (is_tree(start_index(ivma))) THEN
453	DO iagec = nagec_pft(ivma),1,-1
454
455	!start from oldest age class and then move to younger age classes.
456	ipft = start_index(ivma)+iagec-1
457
458	! Check whether woodmass exceeds boundaries of
459	! the age class.
460	IF(ld_agec)THEN
461	WRITE(numout,*) 'Checking to merge for: '
462	WRITE(numout,*) 'ipft,iagec,ipts: ',ipft,iagec,ipts
463	WRITE(numout,*) 'nagec_pft,woodmass,age_class_bound: ',nagec_pft(ivma),&
464	woodmass(ipts,ipft),bound_spa(ipts,ipft)
465	ENDIF
466
467	IF ( (iagec .EQ. nagec_pft(ivma)) .AND. &
468	woodmass(ipts,ipft) .GT. bound_spa(ipts,ipft) ) THEN
469
470	! If these conditions are satisfied our woodmass is
471	! very unrealist
472	WRITE(numout,*) 'WARNING: age class indicator exceeds: ', &
473	bound_spa(ipts,ipft)
474
475	ELSEIF ( (iagec .NE. nagec_pft(ivma)) .AND. &
476	woodmass(ipts,ipft) .GT. bound_spa(ipts,ipft)) THEN
477
478	IF(ld_agec)THEN
479	WRITE(numout,*) 'Merging biomass'
480	WRITE(numout,*) 'ipts,ipft,iagec: ',ipts,ipft,iagec
481	WRITE(numout,*) 'age_class_bound: ',bound_spa(ipts,ipft)
482	WRITE(numout,*) 'woodmass: ',woodmass(ipts,ipft)
483
484	ENDIF
485
486	!! 2 Merge biomass
487	! Biomass of two age classes needs to be merged. The established
488	! vegetation is stored in ipft+1, the new vegetation is stored in
489	! ipft
490	share_expanded = veget_max(ipts,ipft+1) / &
491	( veget_max(ipts,ipft+1) + veget_max(ipts,ipft) )
492	! We also need a scaling factor which includes the litter
493	DO ilev=1,nlevs
494	DO ilit=1,nlitt
495	IF(litter(ipts,ilit,ipft,ilev,icarbon) .GE. min_stomate)THEN
496	litter_weight_expanded(ilit,ilev)=litter(ipts,ilit,ipft+1,ilev,icarbon) * veget_max(ipts,ipft+1)/ &
497	(litter(ipts,ilit,ipft+1,ilev,icarbon) * veget_max(ipts,ipft+1) + &
498	litter(ipts,ilit,ipft,ilev,icarbon) * veget_max(ipts,ipft))
499	ELSE
500	litter_weight_expanded(ilit,ilev)=zero
501	ENDIF
502	END DO
503	ENDDO
504
505
506
507	! Merge the biomass and ind of the two age classes
508	biomass(ipts,ipft+1,:,:) = share_expanded * biomass(ipts,ipft+1,:,:) + &
509	(un - share_expanded) * biomass(ipts,ipft,:,:)
510	ind(ipts,ipft+1) = share_expanded * ind(ipts,ipft+1) + &
511	(un - share_expanded) * ind(ipts,ipft)
512
513	!! 3 Empty the age class that was merged and update veget_max
514	ind(ipts,ipft) = zero
515	biomass(ipts,ipft,:,:) = zero
516	veget_max(ipts,ipft+1) = veget_max(ipts,ipft+1) + veget_max(ipts,ipft)
517	veget_max(ipts,ipft) = zero
518
519	!! 4 Calculate the PFT characteristics of the merged PFT
520	! Take the weighted mean of the existing vegetation and the new
521	! vegetation joining this PFT.
522	! Note that co2_to_bm is in gC. m-2 dt-1 ,
523	! so we should also take the weighted mean (rather than sum if
524	! this where absolute values).
525	lm_lastyearmax(ipts,ipft+1) = share_expanded * lm_lastyearmax(ipts,ipft+1) + &
526	(un - share_expanded) * lm_lastyearmax(ipts,ipft)
527	lm_lastyearmax(ipts,ipft) = zero
528	!age(ipts,ipft+1) = share_expanded * age(ipts,ipft+1) + &
529	! (un - share_expanded) * age(ipts,ipft)
530	!age(ipts,ipft) = zero
531
532	!CHECK: more strictly this should be considered together with leaf mass
533	leaf_frac(ipts,ipft+1,:) = share_expanded * leaf_frac(ipts,ipft+1,:) + &
534	(un - share_expanded) * leaf_frac(ipts,ipft,:)
535	leaf_frac(ipts,ipft,:) = zero
536	leaf_age(ipts,ipft+1,:) = share_expanded * leaf_age(ipts,ipft+1,:) + &
537	(un - share_expanded) * leaf_age(ipts,ipft,:)
538	leaf_age(ipts,ipft,:) = zero
539	co2_to_bm(ipts,ipft+1) = share_expanded * co2_to_bm(ipts,ipft+1) + &
540	(un - share_expanded) * co2_to_bm(ipts,ipft)
541	co2_to_bm(ipts,ipft) = zero
542
543	! Everywhere deals with the migration of vegetation. Copy the
544	! status of the most migrated vegetation for the whole PFT
545	everywhere(ipts,ipft+1) = MAX(everywhere(ipts,ipft), everywhere(ipts,ipft+1))
546	everywhere(ipts,ipft) = zero
547
548	! The new soil&litter pools are the weighted mean of the newly
549	! established vegetation for that PFT and the soil&litter pools
550	! of the original vegetation that already exists in that PFT.
551	! Since it is not only the amount of vegetation present (veget_max) but also
552	! the amount of structural litter (litter) that is important, we have to
553	! weight by both items here.
554	DO ilev=1,nlevs
555	lignin_struc(ipts,ipft+1,ilev) = litter_weight_expanded(istructural,ilev) * lignin_struc(ipts,ipft+1,ilev) + &
556	(un - litter_weight_expanded(istructural,ilev)) * lignin_struc(ipts,ipft,ilev)
557	lignin_struc(ipts,ipft,ilev) = zero
558	ENDDO
559	litter(ipts,:,ipft+1,:,:) = share_expanded * litter(ipts,:,ipft+1,:,:) + &
560	(un - share_expanded) * litter(ipts,:,ipft,:,:)
561	litter(ipts,:,ipft,:,:) = zero
562
563	fuel_1hr(ipts,ipft+1,:,:) = share_expanded * fuel_1hr(ipts,ipft+1,:,:) + &
564	(un - share_expanded) * fuel_1hr(ipts,ipft,:,:)
565	fuel_1hr(ipts,ipft,:,:) = zero
566
567	fuel_10hr(ipts,ipft+1,:,:) = share_expanded * fuel_10hr(ipts,ipft+1,:,:) + &
568	(un - share_expanded) * fuel_10hr(ipts,ipft,:,:)
569	fuel_10hr(ipts,ipft,:,:) = zero
570
571	fuel_100hr(ipts,ipft+1,:,:) = share_expanded * fuel_100hr(ipts,ipft+1,:,:) + &
572	(un - share_expanded) * fuel_100hr(ipts,ipft,:,:)
573	fuel_100hr(ipts,ipft,:,:) = zero
574
575	fuel_1000hr(ipts,ipft+1,:,:) = share_expanded * fuel_1000hr(ipts,ipft+1,:,:) + &
576	(un - share_expanded) * fuel_1000hr(ipts,ipft,:,:)
577	fuel_1000hr(ipts,ipft,:,:) = zero
578
579	carbon(ipts,:,ipft+1) = share_expanded * carbon(ipts,:,ipft+1) + &
580	(un - share_expanded) * carbon(ipts,:,ipft)
581	carbon(ipts,:,ipft) = zero
582
583	deepC_a(ipts,:,ipft+1) = share_expanded * deepC_a(ipts,:,ipft+1) + &
584	(un - share_expanded) * deepC_a(ipts,:,ipft)
585	deepC_a(ipts,:,ipft) = zero
586
587	deepC_s(ipts,:,ipft+1) = share_expanded * deepC_s(ipts,:,ipft+1) + &
588	(un - share_expanded) * deepC_s(ipts,:,ipft)
589	deepC_s(ipts,:,ipft) = zero
590
591	deepC_p(ipts,:,ipft+1) = share_expanded * deepC_p(ipts,:,ipft+1) + &
592	(un - share_expanded) * deepC_p(ipts,:,ipft)
593	deepC_p(ipts,:,ipft) = zero
594
595	bm_to_litter(ipts,ipft+1,:,:) = share_expanded * bm_to_litter(ipts,ipft+1,:,:) + &
596	(un - share_expanded) * bm_to_litter(ipts,ipft,:,:)
597	bm_to_litter(ipts,ipft,:,:) = zero
598
599	! Copy variables that depend on veget_max
600	when_growthinit(ipts,ipft+1) = share_expanded * when_growthinit(ipts,ipft+1) + &
601	(un - share_expanded) * when_growthinit(ipts,ipft)
602	when_growthinit(ipts,ipft) = zero
603	gdd_from_growthinit(ipts,ipft+1) = share_expanded * &
604	gdd_from_growthinit(ipts,ipft+1) + &
605	(un - share_expanded) * gdd_from_growthinit(ipts,ipft)
606	gdd_from_growthinit(ipts,ipft) = zero
607	gdd_midwinter(ipts,ipft+1) = share_expanded * gdd_midwinter(ipts,ipft+1) + &
608	(un - share_expanded) * gdd_midwinter(ipts,ipft)
609	gdd_midwinter(ipts,ipft) = zero
610	time_hum_min(ipts,ipft+1) = share_expanded * time_hum_min(ipts,ipft+1) + &
611	(un - share_expanded) * time_hum_min(ipts,ipft)
612	time_hum_min(ipts,ipft) = zero
613	gdd_m5_dormance(ipts,ipft+1) = share_expanded * gdd_m5_dormance(ipts,ipft+1) + &
614	(un - share_expanded) * gdd_m5_dormance(ipts,ipft)
615	gdd_m5_dormance(ipts,ipft) = zero
616	ncd_dormance(ipts,ipft+1) = share_expanded * ncd_dormance(ipts,ipft+1) + &
617	(un - share_expanded) * ncd_dormance(ipts,ipft)
618	ncd_dormance(ipts,ipft) = zero
619	moiavail_month(ipts,ipft+1) = share_expanded * moiavail_month(ipts,ipft+1) + &
620	(un - share_expanded) * moiavail_month(ipts,ipft)
621	moiavail_month(ipts,ipft) = zero
622	moiavail_week(ipts,ipft+1) = share_expanded * moiavail_week(ipts,ipft+1) + &
623	(un - share_expanded) * moiavail_week(ipts,ipft)
624	moiavail_week(ipts,ipft) = zero
625	ngd_minus5(ipts,ipft+1) = share_expanded * ngd_minus5(ipts,ipft+1) + &
626	(un - share_expanded) * ngd_minus5(ipts,ipft)
627	ngd_minus5(ipts,ipft) = zero
628
629	! Copy remaining properties
630	PFTpresent(ipts,ipft+1) = PFTpresent(ipts,ipft)
631	PFTpresent(ipts,ipft) = .FALSE.
632	senescence(ipts,ipft+1) = senescence(ipts,ipft)
633	senescence(ipts,ipft) = .FALSE.
634	npp_longterm(ipts,ipft+1) = share_expanded * npp_longterm(ipts,ipft+1) + &
635	(un - share_expanded) * npp_longterm(ipts,ipft)
636	npp_longterm(ipts,ipft) = zero
637	gpp_daily(ipts,ipft+1) = share_expanded * gpp_daily(ipts,ipft+1) + &
638	(un - share_expanded) * gpp_daily(ipts,ipft)
639	gpp_daily(ipts,ipft) = zero
640	gpp_week(ipts,ipft+1) = share_expanded * gpp_week(ipts,ipft+1) + &
641	(un - share_expanded) * gpp_week(ipts,ipft)
642	gpp_week(ipts,ipft) = zero
643	resp_maint(ipts,ipft+1) = share_expanded * resp_maint(ipts,ipft+1) + &
644	(un - share_expanded) * resp_maint(ipts,ipft)
645	resp_maint(ipts,ipft) = zero
646	resp_growth(ipts,ipft+1) = share_expanded * resp_growth(ipts,ipft+1) + &
647	(un - share_expanded) * resp_growth(ipts,ipft)
648	resp_growth(ipts,ipft) = zero
649	npp_daily(ipts,ipft+1) = share_expanded * npp_daily(ipts,ipft+1) + &
650	(un - share_expanded) * npp_daily(ipts,ipft)
651	npp_daily(ipts,ipft) = zero
652
653	ENDIF
654	ENDDO
655	! concerned MTC is grass/pasture/crop
656	ELSE
657	DO iagec = 1,nagec_pft(ivma),1
658
659	! As the soil C gets smaller when forest-generating crop gets older,
660	! we start from young age class and then move to older age classes.
661	! If the soil C of ipft is smaller than the threshold, then it should
662	! go to the next age class.
663	ipft = start_index(ivma)+iagec-1
664
665	! Check whether woodmass exceeds boundaries of
666	! the age class.
667	IF(ld_agec)THEN
668	WRITE(numout,*) 'Checking to merge for: '
669	WRITE(numout,*) 'ipft,iagec,ipts: ',ipft,iagec,ipts
670	WRITE(numout,*) 'nagec_pft,woodmass,age_class_bound: ',nagec_pft(ivma),&
671	woodmass(ipts,ipft),bound_spa(ipts,ipft)
672	ENDIF
673
674	!IF ( (iagec .EQ. 1) .AND. &
675	! woodmass(ipts,ipft) .GT. bound_spa(ipts,ipft) ) THEN
676	!
677	! ! If this is satisfied than we're having a quite large
678	! ! soil C in the newly initiated crop
679	! WRITE(numout,*) 'WARNING: age class indicator exceeds: ', &
680	! bound_spa(ipts,ipft)
681
682	!ELSEIF ( (iagec .NE. nagec_pft(ivma)) .AND. &
683	! woodmass(ipts,ipft) .LT. bound_spa(ipts,ipft)) THEN
684
685	! If the soil C is smaller than the threshold and the concerned
686	! ipft is not the oldest age class, then it should move to the
687	! next (older) age class. So we have to set the soil C threshold
688	! for crop as:
689
690	! youngest: 0.9 of maximum end-spinup forest soil C
691	! 2nd young: 0.75 of maximum end-spniup forest soil C
692	! old: 0.55 of maximum end-spniup forest soil C
693	! oldest: the oldest one should not be less than zero.
694	IF ( (iagec .NE. nagec_pft(ivma)) .AND. &
695	woodmass(ipts,ipft) .LT. bound_spa(ipts,ipft) .AND. veget_max(ipts,ipft) .GT. min_stomate) THEN
696	IF(ld_agec)THEN
697	WRITE(numout,*) 'Merging biomass'
698	WRITE(numout,*) 'ipts,ipft,iagec: ',ipts,ipft,iagec
699	WRITE(numout,*) 'age_class_bound: ',bound_spa(ipts,ipft)
700	WRITE(numout,*) 'woodmass: ',woodmass(ipts,ipft)
701
702	ENDIF
703
704	!! 2 Merge biomass
705	! Biomass of two age classes needs to be merged. The established
706	! vegetation is stored in ipft+1, the new vegetation is stored in
707	! ipft
708	share_expanded = veget_max(ipts,ipft+1) / &
709	( veget_max(ipts,ipft+1) + veget_max(ipts,ipft) )
710	! We also need a scaling factor which includes the litter
711	DO ilev=1,nlevs
712	DO ilit=1,nlitt
713	IF(litter(ipts,ilit,ipft,ilev,icarbon) .GE. min_stomate)THEN
714	litter_weight_expanded(ilit,ilev)=litter(ipts,ilit,ipft+1,ilev,icarbon) * veget_max(ipts,ipft+1)/ &
715	(litter(ipts,ilit,ipft+1,ilev,icarbon) * veget_max(ipts,ipft+1) + &
716	litter(ipts,ilit,ipft,ilev,icarbon) * veget_max(ipts,ipft))
717	ELSE
718	litter_weight_expanded(ilit,ilev)=zero
719	ENDIF
720	END DO
721	ENDDO
722
723	! Merge the biomass and ind of the two age classes
724	biomass(ipts,ipft+1,:,:) = share_expanded * biomass(ipts,ipft+1,:,:) + &
725	(un - share_expanded) * biomass(ipts,ipft,:,:)
726	ind(ipts,ipft+1) = share_expanded * ind(ipts,ipft+1) + &
727	(un - share_expanded) * ind(ipts,ipft)
728
729	!! 3 Empty the age class that was merged and update veget_max
730	ind(ipts,ipft) = zero
731	biomass(ipts,ipft,:,:) = zero
732	veget_max(ipts,ipft+1) = veget_max(ipts,ipft+1) + veget_max(ipts,ipft)
733	veget_max(ipts,ipft) = zero
734
735	!! 4 Calculate the PFT characteristics of the merged PFT
736	! Take the weighted mean of the existing vegetation and the new
737	! vegetation joining this PFT.
738	! Note that co2_to_bm is in gC. m-2 dt-1 ,
739	! so we should also take the weighted mean (rather than sum if
740	! this where absolute values).
741	lm_lastyearmax(ipts,ipft+1) = share_expanded * lm_lastyearmax(ipts,ipft+1) + &
742	(un - share_expanded) * lm_lastyearmax(ipts,ipft)
743	lm_lastyearmax(ipts,ipft) = zero
744	!age(ipts,ipft+1) = share_expanded * age(ipts,ipft+1) + &
745	! (un - share_expanded) * age(ipts,ipft)
746	!age(ipts,ipft) = zero
747
748	!CHECK: more strictly this should be considered together with leaf mass
749	leaf_frac(ipts,ipft+1,:) = share_expanded * leaf_frac(ipts,ipft+1,:) + &
750	(un - share_expanded) * leaf_frac(ipts,ipft,:)
751	leaf_frac(ipts,ipft,:) = zero
752	leaf_age(ipts,ipft+1,:) = share_expanded * leaf_age(ipts,ipft+1,:) + &
753	(un - share_expanded) * leaf_age(ipts,ipft,:)
754	leaf_age(ipts,ipft,:) = zero
755	co2_to_bm(ipts,ipft+1) = share_expanded * co2_to_bm(ipts,ipft+1) + &
756	(un - share_expanded) * co2_to_bm(ipts,ipft)
757	co2_to_bm(ipts,ipft) = zero
758
759	! Everywhere deals with the migration of vegetation. Copy the
760	! status of the most migrated vegetation for the whole PFT
761	everywhere(ipts,ipft+1) = MAX(everywhere(ipts,ipft), everywhere(ipts,ipft+1))
762	everywhere(ipts,ipft) = zero
763
764	! The new soil&litter pools are the weighted mean of the newly
765	! established vegetation for that PFT and the soil&litter pools
766	! of the original vegetation that already exists in that PFT.
767	! Since it is not only the amount of vegetation present (veget_max) but also
768	! the amount of structural litter (litter) that is important, we have to
769	! weight by both items here.
770	DO ilev=1,nlevs
771	lignin_struc(ipts,ipft+1,ilev) = litter_weight_expanded(istructural,ilev) * lignin_struc(ipts,ipft+1,ilev) + &
772	(un - litter_weight_expanded(istructural,ilev)) * lignin_struc(ipts,ipft,ilev)
773	lignin_struc(ipts,ipft,ilev) = zero
774	ENDDO
775	litter(ipts,:,ipft+1,:,:) = share_expanded * litter(ipts,:,ipft+1,:,:) + &
776	(un - share_expanded) * litter(ipts,:,ipft,:,:)
777	litter(ipts,:,ipft,:,:) = zero
778
779	fuel_1hr(ipts,ipft+1,:,:) = share_expanded * fuel_1hr(ipts,ipft+1,:,:) + &
780	(un - share_expanded) * fuel_1hr(ipts,ipft,:,:)
781	fuel_1hr(ipts,ipft,:,:) = zero
782
783	fuel_10hr(ipts,ipft+1,:,:) = share_expanded * fuel_10hr(ipts,ipft+1,:,:) + &
784	(un - share_expanded) * fuel_10hr(ipts,ipft,:,:)
785	fuel_10hr(ipts,ipft,:,:) = zero
786
787	fuel_100hr(ipts,ipft+1,:,:) = share_expanded * fuel_100hr(ipts,ipft+1,:,:) + &
788	(un - share_expanded) * fuel_100hr(ipts,ipft,:,:)
789	fuel_100hr(ipts,ipft,:,:) = zero
790
791	fuel_1000hr(ipts,ipft+1,:,:) = share_expanded * fuel_1000hr(ipts,ipft+1,:,:) + &
792	(un - share_expanded) * fuel_1000hr(ipts,ipft,:,:)
793	fuel_1000hr(ipts,ipft,:,:) = zero
794
795	carbon(ipts,:,ipft+1) = share_expanded * carbon(ipts,:,ipft+1) + &
796	(un - share_expanded) * carbon(ipts,:,ipft)
797	carbon(ipts,:,ipft) = zero
798
799	deepC_a(ipts,:,ipft+1) = share_expanded * deepC_a(ipts,:,ipft+1) + &
800	(un - share_expanded) * deepC_a(ipts,:,ipft)
801	deepC_a(ipts,:,ipft) = zero
802
803	deepC_s(ipts,:,ipft+1) = share_expanded * deepC_s(ipts,:,ipft+1) + &
804	(un - share_expanded) * deepC_s(ipts,:,ipft)
805	deepC_s(ipts,:,ipft) = zero
806
807	deepC_p(ipts,:,ipft+1) = share_expanded * deepC_p(ipts,:,ipft+1) + &
808	(un - share_expanded) * deepC_p(ipts,:,ipft)
809	deepC_p(ipts,:,ipft) = zero
810
811	bm_to_litter(ipts,ipft+1,:,:) = share_expanded * bm_to_litter(ipts,ipft+1,:,:) + &
812	(un - share_expanded) * bm_to_litter(ipts,ipft,:,:)
813	bm_to_litter(ipts,ipft,:,:) = zero
814
815	! Copy variables that depend on veget_max
816	when_growthinit(ipts,ipft+1) = share_expanded * when_growthinit(ipts,ipft+1) + &
817	(un - share_expanded) * when_growthinit(ipts,ipft)
818	when_growthinit(ipts,ipft) = zero
819	gdd_from_growthinit(ipts,ipft+1) = share_expanded * &
820	gdd_from_growthinit(ipts,ipft+1) + &
821	(un - share_expanded) * gdd_from_growthinit(ipts,ipft)
822	gdd_from_growthinit(ipts,ipft) = zero
823	gdd_midwinter(ipts,ipft+1) = share_expanded * gdd_midwinter(ipts,ipft+1) + &
824	(un - share_expanded) * gdd_midwinter(ipts,ipft)
825	gdd_midwinter(ipts,ipft) = zero
826	time_hum_min(ipts,ipft+1) = share_expanded * time_hum_min(ipts,ipft+1) + &
827	(un - share_expanded) * time_hum_min(ipts,ipft)
828	time_hum_min(ipts,ipft) = zero
829	gdd_m5_dormance(ipts,ipft+1) = share_expanded * gdd_m5_dormance(ipts,ipft+1) + &
830	(un - share_expanded) * gdd_m5_dormance(ipts,ipft)
831	gdd_m5_dormance(ipts,ipft) = zero
832	ncd_dormance(ipts,ipft+1) = share_expanded * ncd_dormance(ipts,ipft+1) + &
833	(un - share_expanded) * ncd_dormance(ipts,ipft)
834	ncd_dormance(ipts,ipft) = zero
835	moiavail_month(ipts,ipft+1) = share_expanded * moiavail_month(ipts,ipft+1) + &
836	(un - share_expanded) * moiavail_month(ipts,ipft)
837	moiavail_month(ipts,ipft) = zero
838	moiavail_week(ipts,ipft+1) = share_expanded * moiavail_week(ipts,ipft+1) + &
839	(un - share_expanded) * moiavail_week(ipts,ipft)
840	moiavail_week(ipts,ipft) = zero
841	ngd_minus5(ipts,ipft+1) = share_expanded * ngd_minus5(ipts,ipft+1) + &
842	(un - share_expanded) * ngd_minus5(ipts,ipft)
843	ngd_minus5(ipts,ipft) = zero
844
845	! Copy remaining properties
846	PFTpresent(ipts,ipft+1) = PFTpresent(ipts,ipft)
847	PFTpresent(ipts,ipft) = .FALSE.
848	senescence(ipts,ipft+1) = senescence(ipts,ipft)
849	senescence(ipts,ipft) = .FALSE.
850	npp_longterm(ipts,ipft+1) = share_expanded * npp_longterm(ipts,ipft+1) + &
851	(un - share_expanded) * npp_longterm(ipts,ipft)
852	npp_longterm(ipts,ipft) = zero
853	gpp_daily(ipts,ipft+1) = share_expanded * gpp_daily(ipts,ipft+1) + &
854	(un - share_expanded) * gpp_daily(ipts,ipft)
855	gpp_daily(ipts,ipft) = zero
856	gpp_week(ipts,ipft+1) = share_expanded * gpp_week(ipts,ipft+1) + &
857	(un - share_expanded) * gpp_week(ipts,ipft)
858	gpp_week(ipts,ipft) = zero
859	resp_maint(ipts,ipft+1) = share_expanded * resp_maint(ipts,ipft+1) + &
860	(un - share_expanded) * resp_maint(ipts,ipft)
861	resp_maint(ipts,ipft) = zero
862	resp_growth(ipts,ipft+1) = share_expanded * resp_growth(ipts,ipft+1) + &
863	(un - share_expanded) * resp_growth(ipts,ipft)
864	resp_growth(ipts,ipft) = zero
865	npp_daily(ipts,ipft+1) = share_expanded * npp_daily(ipts,ipft+1) + &
866	(un - share_expanded) * npp_daily(ipts,ipft)
867	npp_daily(ipts,ipft) = zero
868
869	ENDIF
870	ENDDO
871
872	ENDIF
873
874	END SUBROUTINE check_merge_same_MTC
875
876	! ================================================================================================================================
877	!! SUBROUTINE : harvest_forest
878	!!
879	!>\BRIEF : Handle forest harvest before its legacy is transferred to
880	! newly initialized youngest-age-class PFT.
881	!!
882	!>\DESCRIPTION
883	!_ ================================================================================================================================
884	!!++TEMP++ biomass,veget_frac are not used because the remaining biomass to be
885	!! harvested is calculated within the deforestation fire module.
886	SUBROUTINE harvest_forest (npts,ipts,ivm,biomass,frac, &
887	litter, deforest_biomass_remain,&
888	fuel_1hr,fuel_10hr,&
889	fuel_100hr,fuel_1000hr,&
890	lignin_struc,&
891	bm_to_litter_pro,convflux,prod10,prod100,&
892	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
893	fuel_1000hr_pro, lignin_content_pro)
894
895
896	IMPLICIT NONE
897
898	!! 0.1 Input variables
899	INTEGER, INTENT(in) :: npts
900	INTEGER, INTENT(in) :: ipts
901	INTEGER, INTENT(in) :: ivm
902	REAL(r_std), INTENT(in) :: frac !! the fraction of land covered by forest to be deforested
903	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
904	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
905	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
906	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
907	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
908	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! Vegetmax-weighted remaining litter on the ground for
909	!! deforestation region.
910	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
911	!! deforestation region.
912	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
913	!! above and below ground
914
915	!! 0.2 Modified variables
916	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro !! conversion of biomass to litter
917	!! @tex ($gC m^{-2} day^{-1}$) @endtex
918	REAL(r_std), DIMENSION(:), INTENT(inout) :: convflux !! release during first year following land cover
919	!! change
920
921	REAL(r_std), DIMENSION(npts,0:10), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
922	!! pool after the annual release for each
923	!! compartment (10 + 1 : input from year of land
924	!! cover change)
925	REAL(r_std), DIMENSION(npts,0:100), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
926	!! pool after the annual release for each
927	!! compartment (100 + 1 : input from year of land
928	!! cover change)
929
930	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: litter_pro
931	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1hr_pro
932	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_10hr_pro
933	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_100hr_pro
934	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1000hr_pro
935	REAL(r_std), DIMENSION(:),INTENT(inout) :: lignin_content_pro
936
937
938
939	!! 0.4 Local variables
940	REAL(r_std) :: above
941
942	! harvest of aboveground sap- and heartwood biomass after taking into
943	! account of deforestation fire
944	IF (allow_deforest_fire) THEN
945	above = deforest_biomass_remain(ipts,ivm,isapabove,icarbon)+ &
946	deforest_biomass_remain(ipts,ivm,iheartabove,icarbon)
947	convflux(ipts) = convflux(ipts) + 0
948	prod10(ipts,0) = prod10(ipts,0) + 0.4*above
949	prod100(ipts,0) = prod100(ipts,0) + 0.6*above
950	ELSE
951	above = (biomass(ipts,ivm,isapabove,icarbon)+ &
952	biomass(ipts,ivm,iheartabove,icarbon))*frac
953	convflux(ipts) = convflux(ipts) + coeff_lcchange_1(ivm) * above
954	prod10(ipts,0) = prod10(ipts,0) + coeff_lcchange_10(ivm) * above
955	prod100(ipts,0) = prod100(ipts,0) + coeff_lcchange_100(ivm) * above
956	ENDIF
957
958	! the transfer of dead biomass to litter
959	bm_to_litter_pro(isapbelow,:) = bm_to_litter_pro(isapbelow,:) + &
960	biomass(ipts,ivm,isapbelow,:)*frac
961	bm_to_litter_pro(iheartbelow,:) = bm_to_litter_pro(iheartbelow,:) + &
962	biomass(ipts,ivm,iheartbelow,:)*frac
963	bm_to_litter_pro(iroot,:) = bm_to_litter_pro(iroot,:) + &
964	biomass(ipts,ivm,iroot,:)*frac
965	bm_to_litter_pro(ifruit,:) = bm_to_litter_pro(ifruit,:) + &
966	biomass(ipts,ivm,ifruit,:)*frac
967	bm_to_litter_pro(icarbres,:) = bm_to_litter_pro(icarbres,:) + &
968	biomass(ipts,ivm,icarbres,:)*frac
969	bm_to_litter_pro(ileaf,:) = bm_to_litter_pro(ileaf,:) + &
970	biomass(ipts,ivm,ileaf,:)*frac
971
972	!update litter_pro
973	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
974	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
975	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
976	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
977	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
978	!don't forget to hanle litter lignin content
979	lignin_content_pro(:)= lignin_content_pro(:) + &
980	litter(ipts,istructural,ivm,:,icarbon)fraclignin_struc(ipts,ivm,:)
981
982	END SUBROUTINE harvest_forest
983
984	! ================================================================================================================================
985	!! SUBROUTINE : harvest_herb
986	!!
987	!>\BRIEF : Handle herbaceous PFT clearing before its legacy is transferred to
988	! newly initialized youngest-age-class PFT.
989	!!
990	!>\DESCRIPTION
991	!_ ================================================================================================================================
992	SUBROUTINE harvest_herb (ipts,ivm,biomass,veget_frac,bm_to_litter_pro)
993
994	IMPLICIT NONE
995
996	!! 0.1 Input variables
997	INTEGER, INTENT(in) :: ipts
998	INTEGER, INTENT(in) :: ivm
999	REAL(r_std), INTENT(in) :: veget_frac !! the fraction of land covered by herbaceous PFT to be cleared
1000	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
1001
1002	!! 0.2 Modified variables
1003	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro
1004
1005
1006
1007	! the transfer of dead biomass to litter
1008	bm_to_litter_pro(:,:) = bm_to_litter_pro(:,:) + biomass(ipts,ivm,:,:)*veget_frac
1009
1010	END SUBROUTINE harvest_herb
1011
1012
1013	! ================================================================================================================================
1014	!! SUBROUTINE : initialize_proxy_pft
1015	!!
1016	!>\BRIEF Initialize a proxy new youngest age class PFT.
1017	!!
1018	!>\DESCRIPTION Initialize a proxy new youngest age class PFT that will be
1019	!! merged with existing yongest age class, or fill the empty
1020	!! niche of the youngest age class PFT.
1021	!_ ================================================================================================================================
1022	SUBROUTINE initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro, &
1023	biomass_pro, co2_to_bm_pro, ind_pro, age_pro, &
1024	senescence_pro, PFTpresent_pro, &
1025	lm_lastyearmax_pro, everywhere_pro, npp_longterm_pro, &
1026	leaf_frac_pro,leaf_age_pro)
1027
1028	IMPLICIT NONE
1029
1030	!! 0.1 Input variables
1031	INTEGER, INTENT(in) :: ipts !!
1032	INTEGER, INTENT(in) :: ipft_young_agec !! index of the concerned youngest-age-class PFT
1033	REAL(r_std), INTENT(in) :: veget_max_pro !! fraction of grid cell land area that's to be occupied
1034
1035	!! 0.2 Modified variables
1036	REAL(r_std), INTENT(inout) :: co2_to_bm_pro
1037
1038	!! 0.3 Output variables
1039	REAL(r_std), DIMENSION(:,:), INTENT(out) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
1040	REAL(r_std), DIMENSION(:), INTENT(out) :: leaf_frac_pro !! fraction of leaves in leaf age class
1041	REAL(r_std), DIMENSION(:), INTENT(out) :: leaf_age_pro !! fraction of leaves in leaf age class
1042	REAL(r_std), INTENT(out) :: age_pro, ind_pro, lm_lastyearmax_pro
1043	REAL(r_std), INTENT(out) :: npp_longterm_pro
1044	REAL(r_std), INTENT(out) :: everywhere_pro !! is the PFT everywhere in the grid box or very
1045	LOGICAL, INTENT(out) :: senescence_pro !! plant senescent (only for deciduous trees) Set
1046	!! to .FALSE. if PFT is introduced or killed
1047	LOGICAL, INTENT(out) :: PFTpresent_pro !! Is pft there (unitless)
1048
1049	!! 0.4 Local variables
1050	!REAL(r_std), DIMENSION(npts,nvm) :: when_growthinit !! how many days ago was the beginning of the
1051	! !! growing season (days)
1052
1053	REAL(r_std), DIMENSION(nparts,nelements) :: bm_new !! biomass increase @tex ($gC m^{-2}$) @endtex
1054	REAL(r_std) :: cn_ind,ind
1055	INTEGER :: i,j,k,l
1056
1057	! -Note-
1058	! This part of codes are copied from the original lcchange_main subroutine
1059	! that initialize a new PFT.
1060
1061	i=ipts
1062	j=ipft_young_agec
1063
1064	!! Initialization of some variables
1065	leaf_frac_pro(:) = zero
1066	leaf_age_pro(:) = zero
1067
1068	!! Initial setting of new establishment
1069	IF (is_tree(j)) THEN
1070	! cn_sapl(j)=0.5; stomate_data.f90
1071	cn_ind = cn_sapl(j)
1072	ELSE
1073	cn_ind = un
1074	ENDIF
1075	ind = veget_max_pro / cn_ind
1076	ind_pro = ind*veget_max_pro
1077	PFTpresent_pro = .TRUE.
1078	senescence_pro = .FALSE.
1079	everywhere_pro = 1.*veget_max_pro
1080	age_pro = zero
1081
1082	! large_value = 1.E33_r_std
1083	! when_growthinit(i,j) = large_value
1084	leaf_frac_pro(1) = 1.0 * veget_max_pro
1085	leaf_age_pro(1) = 1.0 * veget_max_pro !This was not included in original lcchange_main subroutine
1086	npp_longterm_pro = npp_longterm_init * veget_max_pro
1087	lm_lastyearmax_pro = bm_sapl(j,ileaf,icarbon) * ind * veget_max_pro
1088
1089	!! Update of biomass in each each carbon stock component (leaf, sapabove, sapbelow,
1090	!> heartabove, heartbelow, root, fruit, and carbres)\n
1091	DO k = 1, nparts ! loop over # carbon stock components, nparts = 8; stomate_constant.f90
1092	DO l = 1,nelements ! loop over # elements
1093	biomass_pro(k,l) = ind * bm_sapl(j,k,l)
1094	END DO ! loop over # elements
1095	co2_to_bm_pro = co2_to_bm_pro + ind * bm_sapl(j,k,icarbon)
1096	ENDDO ! loop over # carbon stock components
1097
1098	END SUBROUTINE initialize_proxy_pft
1099
1100	! ================================================================================================================================
1101	!! SUBROUTINE sap_take
1102	!!
1103	!>\BRIEF : Take the sapling biomass of the new PFTs from the existing biomass, otherwise
1104	! take from co2_to_bm
1105	!!
1106	!>\DESCRIPTION
1107	!_ ================================================================================================================================
1108	SUBROUTINE sap_take (ipts,ivma,veget_max,biomass_pro,biomass,co2_to_bm_pro)
1109
1110	INTEGER, INTENT(in) :: ipts !!
1111	INTEGER, INTENT(in) :: ivma
1112	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT (LAI ->
1113	REAL(r_std), DIMENSION(:,:), INTENT(in) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
1114
1115	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
1116	REAL(r_std), INTENT(inout) :: co2_to_bm_pro
1117
1118
1119	REAL(r_std), DIMENSION(nparts,nelements) :: biomass_total !! biomass @tex ($gC m^{-2}$) @endtex
1120	REAL(r_std) :: bm_org,bmpro_share
1121	INTEGER :: i,ivm,ipart
1122
1123	biomass_total(:,:) = zero
1124	bm_org = zero
1125	bmpro_share = zero
1126
1127	DO i = 1,nagec_pft(ivma)
1128	ivm = start_index(ivma)+i-1
1129	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
1130	biomass_total = biomass_total + biomass(ipts,ivm,:,:)*veget_max(ipts,ivm)
1131	ENDIF
1132	ENDDO
1133
1134	DO ipart = 1, nparts
1135	IF (biomass_total(ipart,icarbon) .GT. biomass_pro(ipart,icarbon)) THEN
1136	co2_to_bm_pro = co2_to_bm_pro - biomass_pro(ipart,icarbon)
1137	!treat each PFT of the MTC
1138	DO i = 1,nagec_pft(ivma)
1139	ivm = start_index(ivma)+i-1
1140	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
1141	bm_org = biomass(ipts,ivm,ipart,icarbon) * veget_max(ipts,ivm)
1142	bmpro_share = bm_org/biomass_total(ipart,icarbon) * biomass_pro(ipart,icarbon)
1143	biomass(ipts,ivm,ipart,icarbon) = (bm_org - bmpro_share)/veget_max(ipts,ivm)
1144	ENDIF
1145	ENDDO
1146	ENDIF
1147	ENDDO
1148
1149	END SUBROUTINE sap_take
1150
1151	! ================================================================================================================================
1152	!! SUBROUTINE collect_legacy_pft
1153	!!
1154	!>\BRIEF : Collect the legacy variables that are going to be included
1155	! in the newly initialized PFT.
1156	!!
1157	!>\DESCRIPTION
1158	!_ ================================================================================================================================
1159	SUBROUTINE collect_legacy_pft(npts, ipts, ivma, glcc_pftmtc, &
1160	biomass, bm_to_litter, carbon, litter, &
1161	deepC_a, deepC_s, deepC_p, &
1162	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1163	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
1164	resp_maint, resp_growth, resp_hetero, co2_fire, &
1165	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1166	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1167	deforest_litter_remain, deforest_biomass_remain, &
1168	veget_max_pro, carbon_pro, lignin_struc_pro, litter_pro, &
1169	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
1170	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
1171	bm_to_litter_pro, co2_to_bm_pro, gpp_daily_pro, &
1172	npp_daily_pro, resp_maint_pro, resp_growth_pro, &
1173	resp_hetero_pro, co2_fire_pro, &
1174	convflux,prod10,prod100)
1175
1176	IMPLICIT NONE
1177
1178	!! 0.1 Input variables
1179	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
1180	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
1181	INTEGER, INTENT(in) :: ivma !! Index for metaclass
1182	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
1183	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
1184	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: bm_to_litter !! Transfer of biomass to litter
1185	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1186	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: carbon !! carbon pool: active, slow, or passive
1187	!! @tex ($gC m^{-2}$) @endtex
1188	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_a !! Permafrost soil carbon (g/m**3) active
1189	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1190	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1191	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! metabolic and structural litter, above and
1192	!! below ground @tex ($gC m^{-2}$) @endtex
1193	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
1194	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
1195	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
1196	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
1197	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1198	!! above and below ground
1199	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_to_bm !! biomass uptaken
1200	!! @tex ($gC m^{-2} day^{-1}$) @endtex
1201	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gpp_daily !! Daily gross primary productivity
1202	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1203	REAL(r_std), DIMENSION(:,:), INTENT(in) :: npp_daily !! Net primary productivity
1204	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1205	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_maint !! Maintenance respiration
1206	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1207	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_growth !! Growth respiration
1208	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1209	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_hetero !! Heterotrophic respiration
1210	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1211	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_fire !! Heterotrophic respiration
1212	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1213	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_1hr_remain
1214	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_10hr_remain
1215	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_100hr_remain
1216	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_1000hr_remain
1217	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
1218	!! deforestation region.
1219	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
1220	!! deforestation region.
1221
1222	!! 0.2 Output variables
1223	REAL(r_std), DIMENSION(:), INTENT(out) :: carbon_pro
1224	REAL(r_std), DIMENSION(:), INTENT(out) :: deepC_a_pro
1225	REAL(r_std), DIMENSION(:), INTENT(out) :: deepC_s_pro
1226	REAL(r_std), DIMENSION(:), INTENT(out) :: deepC_p_pro
1227	REAL(r_std), DIMENSION(:), INTENT(out) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
1228	!! above and below ground
1229	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: litter_pro
1230	REAL(r_std), DIMENSION(:,:), INTENT(out) :: fuel_1hr_pro
1231	REAL(r_std), DIMENSION(:,:), INTENT(out) :: fuel_10hr_pro
1232	REAL(r_std), DIMENSION(:,:), INTENT(out) :: fuel_100hr_pro
1233	REAL(r_std), DIMENSION(:,:), INTENT(out) :: fuel_1000hr_pro
1234	REAL(r_std), DIMENSION(:,:), INTENT(out) :: bm_to_litter_pro
1235	REAL(r_std), INTENT(out) :: veget_max_pro, co2_to_bm_pro
1236	REAL(r_std), INTENT(out) :: gpp_daily_pro, npp_daily_pro
1237	REAL(r_std), INTENT(out) :: resp_maint_pro, resp_growth_pro
1238	REAL(r_std), INTENT(out) :: resp_hetero_pro, co2_fire_pro
1239
1240	!! 0.3 Modified variables
1241	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: convflux !! release during first year following land cover
1242	!! change
1243
1244	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
1245	!! pool after the annual release for each
1246	!! compartment (10 + 1 : input from year of land
1247	!! cover change)
1248	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
1249	!! pool after the annual release for each
1250	!! compartment (100 + 1 : input from year of land
1251	!! cover change)
1252
1253	!! 0.4 Local variables
1254	REAL(r_std), DIMENSION(nlevs) :: lignin_content_pro
1255	REAL(r_std) :: frac
1256	INTEGER :: ivm
1257
1258
1259	! All *_pro variables collect the legacy pools/fluxes of the ancestor
1260	! PFTs for the receiving youngest age class. All *_pro variables
1261	! represent the quantity weighted by the fraction of ancestor contributing
1262	! PFTs.
1263	! Exceptions:
1264	! lignin_struc_pro:: the ratio of lignin content in structural litter.
1265
1266	veget_max_pro=zero
1267	carbon_pro(:)=zero
1268	deepC_a_pro(:)=zero
1269	deepC_s_pro(:)=zero
1270	deepC_p_pro(:)=zero
1271	lignin_struc_pro(:)=zero
1272	lignin_content_pro(:)=zero
1273	litter_pro(:,:,:)=zero
1274	fuel_1hr_pro(:,:)=zero
1275	fuel_10hr_pro(:,:)=zero
1276	fuel_100hr_pro(:,:)=zero
1277	fuel_1000hr_pro(:,:)=zero
1278	bm_to_litter_pro(:,:)=zero
1279	co2_to_bm_pro=zero
1280	gpp_daily_pro=zero
1281	npp_daily_pro=zero
1282	resp_maint_pro=zero
1283	resp_growth_pro=zero
1284	resp_hetero_pro=zero
1285	co2_fire_pro=zero
1286
1287	DO ivm = 1,nvm
1288	frac = glcc_pftmtc(ipts,ivm,ivma)
1289	IF (frac>zero) THEN
1290	veget_max_pro = veget_max_pro+frac
1291
1292	IF (is_tree(ivm)) THEN
1293	IF (is_tree(start_index(ivma))) THEN
1294	CALL harvest_forest (npts,ipts,ivm,biomass,frac, &
1295	litter, deforest_biomass_remain,&
1296	fuel_1hr,fuel_10hr,&
1297	fuel_100hr,fuel_1000hr,&
1298	lignin_struc,&
1299	bm_to_litter_pro,convflux(:,iwphar),prod10(:,:,iwphar),prod100(:,:,iwphar),&
1300	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
1301	fuel_1000hr_pro, lignin_content_pro)
1302	ELSE
1303	CALL harvest_forest (npts,ipts,ivm,biomass,frac, &
1304	litter, deforest_biomass_remain,&
1305	fuel_1hr,fuel_10hr,&
1306	fuel_100hr,fuel_1000hr,&
1307	lignin_struc,&
1308	bm_to_litter_pro,convflux(:,iwplcc),prod10(:,:,iwplcc),prod100(:,:,iwplcc),&
1309	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
1310	fuel_1000hr_pro, lignin_content_pro)
1311	ENDIF
1312	ELSE
1313	CALL harvest_herb(ipts,ivm,biomass,frac, &
1314	bm_to_litter_pro)
1315	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
1316	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
1317	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
1318	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
1319	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
1320	!don't forget to hanle litter lignin content
1321	lignin_content_pro(:)= lignin_content_pro(:) + &
1322	litter(ipts,istructural,ivm,:,icarbon)lignin_struc(ipts,ivm,:)frac
1323	ENDIF
1324
1325	!! scalar variables to be accumulated and inherited
1326	!! by the destination PFT
1327	bm_to_litter_pro(:,:) = bm_to_litter_pro(:,:) + &
1328	bm_to_litter(ipts,ivm,:,:)*frac
1329	carbon_pro(:) = carbon_pro(:)+carbon(ipts,:,ivm)*frac
1330	deepC_a_pro(:) = deepC_a_pro(:)+deepC_a(ipts,:,ivm)*frac
1331	deepC_s_pro(:) = deepC_s_pro(:)+deepC_s(ipts,:,ivm)*frac
1332	deepC_p_pro(:) = deepC_p_pro(:)+deepC_p(ipts,:,ivm)*frac
1333	co2_to_bm_pro = co2_to_bm_pro + co2_to_bm(ipts,ivm)*frac
1334
1335	gpp_daily_pro = gpp_daily_pro + gpp_daily(ipts,ivm)*frac
1336	npp_daily_pro = npp_daily_pro + npp_daily(ipts,ivm)*frac
1337	resp_maint_pro = resp_maint_pro + resp_maint(ipts,ivm)*frac
1338	resp_growth_pro = resp_growth_pro + resp_growth(ipts,ivm)*frac
1339	resp_hetero_pro = resp_hetero_pro + resp_hetero(ipts,ivm)*frac
1340	co2_fire_pro = co2_fire_pro + co2_fire(ipts,ivm)*frac
1341	ENDIF
1342	ENDDO
1343
1344	WHERE (litter_pro(istructural,:,icarbon) .GT. min_stomate)
1345	lignin_struc_pro(:) = lignin_content_pro(:)/litter_pro(istructural,:,icarbon)
1346	ENDWHERE
1347
1348	END SUBROUTINE collect_legacy_pft
1349
1350
1351	! ================================================================================================================================
1352	!! SUBROUTINE gross_lcchange
1353	!!
1354	!>\BRIEF : Apply gross land cover change.
1355	!!
1356	!>\DESCRIPTION
1357	!_ ================================================================================================================================
1358	SUBROUTINE gross_glcchange_fh (npts, dt_days, harvest_matrix, &
1359	glccSecondShift,glccPrimaryShift,glccNetLCC,&
1360	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1361	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1362	deforest_litter_remain, deforest_biomass_remain, &
1363	convflux, cflux_prod10, cflux_prod100, &
1364	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
1365	veget_max, prod10, prod100, flux10, flux100, &
1366	PFTpresent, senescence, moiavail_month, moiavail_week, &
1367	gpp_week, ngd_minus5, resp_maint, resp_growth, &
1368	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
1369	ind, lm_lastyearmax, everywhere, age, &
1370	co2_to_bm, gpp_daily, co2_fire, &
1371	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
1372	gdd_m5_dormance, ncd_dormance, &
1373	lignin_struc, carbon, leaf_frac, &
1374	deepC_a, deepC_s, deepC_p, &
1375	leaf_age, bm_to_litter, biomass, litter, &
1376	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
1377
1378	IMPLICIT NONE
1379
1380	!! 0.1 Input variables
1381
1382	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
1383	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
1384	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccSecondShift !! the land-cover-change (LCC) matrix in case a gross LCC is
1385	!! used.
1386	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccPrimaryShift !! the land-cover-change (LCC) matrix in case a gross LCC is
1387	!! used.
1388	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
1389	!! used.
1390	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: harvest_matrix !!
1391	!!
1392
1393	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1hr_remain
1394	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_10hr_remain
1395	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_100hr_remain
1396	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1000hr_remain
1397	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
1398	!! deforestation region.
1399	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
1400	!! deforestation region.
1401
1402
1403	!! 0.2 Output variables
1404	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: convflux !! release during first year following land cover
1405	!! change
1406	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: cflux_prod10 !! total annual release from the 10 year-turnover
1407	!! pool @tex ($gC m^{-2}$) @endtex
1408	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: cflux_prod100 !! total annual release from the 100 year-
1409	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
1410	!! after considering the consistency between presribed
1411	!! glcc matrix and existing vegetation fractions.
1412	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
1413	!! there are not enough fractions in the source PFTs
1414	!! /vegetations to target PFTs/vegetations. I.e., these
1415	!! fraction transfers are presribed in LCC matrix but
1416	!! not realized.
1417	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
1418	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout):: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
1419	!! i.e., the contribution of each PFT to the youngest age-class of MTC
1420
1421	!! 0.3 Modified variables
1422	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT (LAI ->
1423	!! infinity) on ground (unitless)
1424	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
1425	!! pool after the annual release for each
1426	!! compartment (10 + 1 : input from year of land
1427	!! cover change)
1428	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
1429	!! pool after the annual release for each
1430	!! compartment (100 + 1 : input from year of land
1431	!! cover change)
1432	REAL(r_std), DIMENSION(npts,10,nwp), INTENT(inout) :: flux10 !! annual release from the 10/100 year-turnover
1433	!! pool compartments
1434	REAL(r_std), DIMENSION(npts,100,nwp), INTENT(inout) :: flux100 !! annual release from the 10/100 year-turnover
1435	!! pool compartments
1436	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
1437	!! each pixel
1438	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
1439	!! for deciduous trees)
1440	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
1441	!! unitless)
1442	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
1443	!! (0 to 1, unitless)
1444	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
1445	!! @tex $(gC m^{-2} day^{-1})$ @endtex
1446	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
1447	!! -5 deg C (for phenology)
1448	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
1449	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1450	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
1451	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1452	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
1453	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1454	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
1455	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1456	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
1457	!! the growing season (days)
1458	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
1459	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
1460	!! @tex $(m^{-2})$ @endtex
1461	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
1462	!! @tex ($gC m^{-2}$) @endtex
1463	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
1464	!! very localized (after its introduction) (?)
1465	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
1466	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
1467	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
1468	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
1469	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1470	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_fire !! Fire carbon emissions
1471	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1472	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
1473	!! availability (days)
1474	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
1475	!! (for phenology) - this is written to the
1476	!! history files
1477	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
1478	!! for crops
1479	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
1480	!! C (for phenology)
1481	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
1482	!! leaves were lost (for phenology)
1483	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1484	!! above and below ground
1485	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
1486	!! @tex ($gC m^{-2}$) @endtex
1487	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
1488	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1489	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1490	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
1491	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
1492	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
1493	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1494	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
1495	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
1496	!! below ground @tex ($gC m^{-2}$) @endtex
1497	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
1498	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
1499	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
1500	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
1501
1502	!! 0.4 Local variables
1503	REAL(r_std), DIMENSION(nparts,nelements) :: bm_to_litter_pro !! conversion of biomass to litter
1504	!! @tex ($gC m^{-2} day^{-1}$) @endtex
1505	REAL(r_std), DIMENSION(nparts,nelements) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
1506	REAL(r_std) :: veget_max_pro !! "maximal" coverage fraction of a PFT (LAI ->
1507	!! infinity) on ground (unitless)
1508	REAL(r_std), DIMENSION(ncarb) :: carbon_pro !! carbon pool: active, slow, or passive
1509	!! @tex ($gC m^{-2}$) @endtex
1510	REAL(r_std), DIMENSION(ndeep) :: deepC_a_pro !! Permafrost carbon pool: active, slow, or passive
1511	!! @tex ($gC m^{-3}$) @endtex
1512	REAL(r_std), DIMENSION(ndeep) :: deepC_s_pro !! Permafrost carbon pool: active, slow, or passive
1513	!! @tex ($gC m^{-3}$) @endtex
1514	REAL(r_std), DIMENSION(ndeep) :: deepC_p_pro !! Permafrost carbon pool: active, slow, or passive
1515	!! @tex ($gC m^{-3}$) @endtex
1516	REAL(r_std), DIMENSION(nlitt,nlevs,nelements) :: litter_pro !! metabolic and structural litter, above and
1517	!! below ground @tex ($gC m^{-2}$) @endtex
1518	REAL(r_std), DIMENSION(nlitt,nelements) :: fuel_1hr_pro
1519	REAL(r_std), DIMENSION(nlitt,nelements) :: fuel_10hr_pro
1520	REAL(r_std), DIMENSION(nlitt,nelements) :: fuel_100hr_pro
1521	REAL(r_std), DIMENSION(nlitt,nelements) :: fuel_1000hr_pro
1522	REAL(r_std), DIMENSION(nlevs) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
1523	!! above and below ground
1524	REAL(r_std), DIMENSION(nleafages) :: leaf_frac_pro !! fraction of leaves in leaf age class
1525	REAL(r_std), DIMENSION(nleafages) :: leaf_age_pro !! fraction of leaves in leaf age class
1526	LOGICAL :: PFTpresent_pro, senescence_pro !! Is pft there (unitless)
1527	REAL(r_std) :: ind_pro, age_pro, lm_lastyearmax_pro, npp_longterm_pro
1528	REAL(r_std) :: everywhere_pro
1529	REAL(r_std) :: gpp_daily_pro, npp_daily_pro, co2_to_bm_pro
1530	REAL(r_std) :: resp_maint_pro, resp_growth_pro
1531	REAL(r_std) :: resp_hetero_pro, co2_fire_pro
1532
1533	INTEGER :: ipts,ivm,ivma,l,m,ipft_young_agec
1534	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
1535
1536	REAL(r_std), DIMENSION(npts,nvmap) :: glcc_mtc !! Increase in fraction of each MTC in its youngest age-class
1537	REAL(r_std), DIMENSION(npts,nvm) :: glccReal_tmp !! A temporary variable to hold glccReal
1538	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf_final !!
1539	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf_final !!
1540	REAL(r_std), DIMENSION(npts) :: pf2yf_compen_sf2yf !!
1541	REAL(r_std), DIMENSION(npts) :: sf2yf_compen_pf2yf !!
1542	REAL(r_std), DIMENSION(npts,nvm) :: glcc_harvest !! Loss of fraction due to forestry harvest
1543
1544	WRITE(numout,*) 'Entering gross_glcchange_fh'
1545	glcc_harvest(:,:) = zero
1546	glccReal_tmp(:,:) = zero
1547
1548	!! Some initialization
1549	convflux(:,:)=zero
1550	prod10(:,0,:) = zero
1551	prod100(:,0,:) = zero
1552	cflux_prod10(:,:) = zero
1553	cflux_prod100(:,:) = zero
1554
1555	CALL gross_glcc_firstday_fh(npts,veget_max,harvest_matrix, &
1556	glccSecondShift,glccPrimaryShift,glccNetLCC,&
1557	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
1558	Deficit_pf2yf_final, Deficit_sf2yf_final, &
1559	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
1560
1561	glcc_mtc(:,:) = SUM(glcc_pftmtc,DIM=2)
1562	DO ipts=1,npts
1563	! Note that we assume people don't intentionally change baresoil to
1564	! vegetated land.
1565	DO ivma = 2,nvmap
1566	! we assume only the youngest age class receives the incoming PFT
1567	! [chaoyuejoy@gmail.com 2015-08-04] This line is commented to allow
1568	! the case of only single age class being handled.
1569	IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
1570	ipft_young_agec = start_index(ivma)
1571
1572	! 1. we accumulate the scalar variables that will be inherited
1573	! note we don't handle the case of harvesting forest because
1574	! we assume glcc_pftmtc(forest->forest) would be zero and this
1575	! case won't occur as it's filtered by the condition of
1576	! (frac>min_stomate)
1577	CALL collect_legacy_pft(npts, ipts, ivma, glcc_pftmtc, &
1578	biomass, bm_to_litter, carbon, litter, &
1579	deepC_a, deepC_s, deepC_p, &
1580	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1581	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
1582	resp_maint, resp_growth, resp_hetero, co2_fire, &
1583	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1584	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1585	deforest_litter_remain, deforest_biomass_remain, &
1586	veget_max_pro, carbon_pro, lignin_struc_pro, litter_pro, &
1587	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
1588	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
1589	bm_to_litter_pro, co2_to_bm_pro, gpp_daily_pro, &
1590	npp_daily_pro, resp_maint_pro, resp_growth_pro, &
1591	resp_hetero_pro, co2_fire_pro, &
1592	convflux,prod10,prod100)
1593
1594	!++TEMP++
1595	! Here we substract the outgoing fraction from the source PFT.
1596	! If a too small fraction remains in this source PFT, then it is
1597	! exhausted, we empty it. The subroutine 'empty_pft' might be
1598	! combined with 'collect_legacy_pft', but now we just put it here.
1599	DO ivm = 1,nvm
1600	IF( glcc_pftmtc(ipts,ivm,ivma)>min_stomate ) THEN
1601	veget_max(ipts,ivm) = veget_max(ipts,ivm)-glcc_pftmtc(ipts,ivm,ivma)
1602	IF ( veget_max(ipts,ivm)<min_stomate ) THEN
1603	CALL empty_pft(ipts, ivm, veget_max, biomass, ind, &
1604	carbon, litter, lignin_struc, bm_to_litter, &
1605	deepC_a, deepC_s, deepC_p, &
1606	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1607	gpp_daily, npp_daily, gpp_week, npp_longterm, &
1608	co2_to_bm, resp_maint, resp_growth, resp_hetero, &
1609	lm_lastyearmax, leaf_frac, leaf_age, age, &
1610	everywhere, PFTpresent, when_growthinit, &
1611	senescence, gdd_from_growthinit, gdd_midwinter, &
1612	time_hum_min, gdd_m5_dormance, ncd_dormance, &
1613	moiavail_month, moiavail_week, ngd_minus5)
1614	ENDIF
1615	ENDIF
1616	ENDDO
1617
1618	! 2. we establish a proxy PFT with the fraction of veget_max_pro,
1619	! which is going to be either merged with existing target
1620	! `ipft_young_agec` PFT, or fill the place if no existing target PFT
1621	! exits.
1622	CALL initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro, &
1623	biomass_pro, co2_to_bm_pro, ind_pro, age_pro, &
1624	senescence_pro, PFTpresent_pro, &
1625	lm_lastyearmax_pro, everywhere_pro, npp_longterm_pro, &
1626	leaf_frac_pro,leaf_age_pro)
1627
1628	CALL sap_take (ipts,ivma,veget_max,biomass_pro,biomass,co2_to_bm_pro)
1629
1630	! 3. we merge the newly initiazlized proxy PFT into existing one
1631	! or use it to fill an empty PFT slot.
1632	CALL add_incoming_proxy_pft(npts, ipts, ipft_young_agec, veget_max_pro,&
1633	carbon_pro, litter_pro, lignin_struc_pro, bm_to_litter_pro, &
1634	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
1635	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
1636	biomass_pro, co2_to_bm_pro, npp_longterm_pro, ind_pro, &
1637	lm_lastyearmax_pro, age_pro, everywhere_pro, &
1638	leaf_frac_pro, leaf_age_pro, PFTpresent_pro, senescence_pro, &
1639	gpp_daily_pro, npp_daily_pro, resp_maint_pro, resp_growth_pro, &
1640	resp_hetero_pro, co2_fire_pro, &
1641	veget_max, carbon, litter, lignin_struc, bm_to_litter, &
1642	deepC_a, deepC_s, deepC_p, &
1643	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1644	biomass, co2_to_bm, npp_longterm, ind, &
1645	lm_lastyearmax, age, everywhere, &
1646	leaf_frac, leaf_age, PFTpresent, senescence, &
1647	gpp_daily, npp_daily, resp_maint, resp_growth, &
1648	resp_hetero, co2_fire)
1649
1650	ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
1651
1652	ENDDO
1653	ENDDO
1654
1655	!! Update 10 year-turnover pool content following flux emission
1656	!! (linear decay (10%) of the initial carbon input)
1657	DO l = 0, 8
1658	m = 10 - l
1659	cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,m,:)
1660	prod10(:,m,:) = prod10(:,m-1,:) - flux10(:,m-1,:)
1661	flux10(:,m,:) = flux10(:,m-1,:)
1662	WHERE (prod10(:,m,:) .LT. 1.0) prod10(:,m,:) = zero
1663	ENDDO
1664
1665	cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,1,:)
1666	flux10(:,1,:) = 0.1 * prod10(:,0,:)
1667	prod10(:,1,:) = prod10(:,0,:)
1668
1669	!! 2.4.3 update 100 year-turnover pool content following flux emission\n
1670	DO l = 0, 98
1671	m = 100 - l
1672	cflux_prod100(:,:) = cflux_prod100(:,:) + flux100(:,m,:)
1673	prod100(:,m,:) = prod100(:,m-1,:) - flux100(:,m-1,:)
1674	flux100(:,m,:) = flux100(:,m-1,:)
1675
1676	WHERE (prod100(:,m,:).LT.1.0) prod100(:,m,:) = zero
1677	ENDDO
1678
1679	cflux_prod100(:,:) = cflux_prod100(:,:) + flux100(:,1,:)
1680	flux100(:,1,:) = 0.01 * prod100(:,0,:)
1681	prod100(:,1,:) = prod100(:,0,:)
1682	prod10(:,0,:) = zero
1683	prod100(:,0,:) = zero
1684
1685	convflux = convflux/one_year*dt_days
1686	cflux_prod10 = cflux_prod10/one_year*dt_days
1687	cflux_prod100 = cflux_prod100/one_year*dt_days
1688
1689	! Write out history files
1690	CALL histwrite_p (hist_id_stomate, 'glcc_pft', itime, &
1691	glcc_pft, npts*nvm, horipft_index)
1692
1693	glccReal_tmp(:,1:12) = glccReal
1694	CALL histwrite_p (hist_id_stomate, 'glccReal', itime, &
1695	glccReal_tmp, npts*nvm, horipft_index)
1696
1697	! Write out forestry harvest variables
1698	DO ipts = 1,npts
1699	DO ivm = 1,nvm
1700	DO ivma = 1,nvmap
1701	IF (is_tree(ivm) .AND. is_tree(start_index(ivma))) THEN
1702	glcc_harvest(ipts,ivm) = glcc_harvest(ipts,ivm) + glcc_pftmtc(ipts,ivm,ivma)
1703	ENDIF
1704	ENDDO
1705	ENDDO
1706	ENDDO
1707	CALL histwrite_p (hist_id_stomate, 'glcc_harvest', itime, &
1708	glcc_harvest, npts*nvm, horipft_index)
1709
1710	glccReal_tmp(:,:) = zero
1711	glccReal_tmp(:,1:12) = IncreDeficit
1712	CALL histwrite_p (hist_id_stomate, 'IncreDeficit', itime, &
1713	glccReal_tmp, npts*nvm, horipft_index)
1714
1715	glccReal_tmp(:,:) = zero
1716	glccReal_tmp(:,1) = Deficit_pf2yf_final
1717	glccReal_tmp(:,2) = Deficit_sf2yf_final
1718	glccReal_tmp(:,3) = pf2yf_compen_sf2yf
1719	glccReal_tmp(:,4) = sf2yf_compen_pf2yf
1720
1721	CALL histwrite_p (hist_id_stomate, 'DefiComForHarvest', itime, &
1722	glccReal_tmp, npts*nvm, horipft_index)
1723
1724	DO ivma = 1, nvmap
1725	WRITE(part_str,'(I2)') ivma
1726	IF (ivma < 10) part_str(1:1) = '0'
1727	CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_'//part_str(1:LEN_TRIM(part_str)), &
1728	itime, glcc_pftmtc(:,:,ivma), npts*nvm, horipft_index)
1729	ENDDO
1730	END SUBROUTINE gross_glcchange_fh
1731
1732
1733	! ================================================================================================================================
1734	!! SUBROUTINE : add_incoming_proxy_pft
1735	!!
1736	!>\BRIEF : Merge the newly incoming proxy PFT cohort with the exisiting
1737	!! cohort.
1738	!! \n
1739	!
1740	!_ ================================================================================================================================
1741	SUBROUTINE add_incoming_proxy_pft(npts, ipts, ipft, veget_max_pro, &
1742	carbon_pro, litter_pro, lignin_struc_pro, bm_to_litter_pro, &
1743	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
1744	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
1745	biomass_pro, co2_to_bm_pro, npp_longterm_pro, ind_pro, &
1746	lm_lastyearmax_pro, age_pro, everywhere_pro, &
1747	leaf_frac_pro, leaf_age_pro, PFTpresent_pro, senescence_pro, &
1748	gpp_daily_pro, npp_daily_pro, resp_maint_pro, resp_growth_pro, &
1749	resp_hetero_pro, co2_fire_pro, &
1750	veget_max, carbon, litter, lignin_struc, bm_to_litter, &
1751	deepC_a, deepC_s, deepC_p, &
1752	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1753	biomass, co2_to_bm, npp_longterm, ind, &
1754	lm_lastyearmax, age, everywhere, &
1755	leaf_frac, leaf_age, PFTpresent, senescence, &
1756	gpp_daily, npp_daily, resp_maint, resp_growth, &
1757	resp_hetero, co2_fire)
1758
1759	IMPLICIT NONE
1760
1761	!! 0.1 Input variables
1762	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
1763	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
1764	INTEGER, INTENT(in) :: ipft
1765	REAL(r_std), INTENT(in) :: veget_max_pro !! The land fraction of incoming new PFTs that are
1766	!! the sum of all its ancestor PFTs
1767	REAL(r_std), DIMENSION(:), INTENT(in) :: carbon_pro
1768	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_a_pro
1769	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_s_pro
1770	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_p_pro
1771	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: litter_pro
1772	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_1hr_pro
1773	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_10hr_pro
1774	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_100hr_pro
1775	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_1000hr_pro
1776	REAL(r_std), DIMENSION(:,:), INTENT(in) :: bm_to_litter_pro
1777	REAL(r_std), DIMENSION(:), INTENT(in) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
1778	!! above and below ground
1779	REAL(r_std), DIMENSION(:,:), INTENT(in) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
1780	REAL(r_std), DIMENSION(:), INTENT(in) :: leaf_frac_pro !! fraction of leaves in leaf age class
1781	REAL(r_std), DIMENSION(:), INTENT(in) :: leaf_age_pro !! fraction of leaves in leaf age class
1782	REAL(r_std), INTENT(in) :: ind_pro, age_pro, lm_lastyearmax_pro
1783	REAL(r_std), INTENT(in) :: npp_longterm_pro, co2_to_bm_pro
1784	REAL(r_std), INTENT(in) :: everywhere_pro !! is the PFT everywhere in the grid box or very
1785	LOGICAL, INTENT(in) :: PFTpresent_pro, senescence_pro !! Is pft there (unitless)
1786
1787	REAL(r_std), INTENT(in) :: gpp_daily_pro, npp_daily_pro
1788	REAL(r_std), INTENT(in) :: resp_maint_pro, resp_growth_pro
1789	REAL(r_std), INTENT(in) :: resp_hetero_pro, co2_fire_pro
1790
1791	!! 0.2 Output variables
1792
1793	!! 0.3 Modified variables
1794
1795	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1796	!! May sum to
1797	!! less than unity if the pixel has
1798	!! nobio area. (unitless, 0-1)
1799
1800	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
1801	!! @tex ($gC m^{-2}$) @endtex
1802	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
1803	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1804	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1805	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
1806	!! below ground @tex ($gC m^{-2}$) @endtex
1807	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
1808	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
1809	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
1810	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
1811	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1812	!! above and below ground
1813	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
1814	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1815	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
1816	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
1817	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
1818
1819	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
1820	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
1821	!! @tex $(m^{-2})$ @endtex
1822	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
1823	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
1824	!! each pixel
1825	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
1826	!! for deciduous trees)
1827	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
1828	!! @tex ($gC m^{-2}$) @endtex
1829	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
1830	!! very localized (after its introduction) (?)
1831	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
1832	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
1833
1834	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
1835	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1836	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
1837	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1838	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
1839	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1840	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
1841	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1842	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
1843	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1844	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_fire !! Heterotrophic respiration
1845	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1846
1847	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
1848	! !! unitless)
1849	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
1850	! !! (0 to 1, unitless)
1851	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
1852	! !! @tex $(gC m^{-2} day^{-1})$ @endtex
1853	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
1854	! !! -5 deg C (for phenology)
1855	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
1856	! !! the growing season (days)
1857	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
1858	! !! availability (days)
1859	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
1860	! !! (for phenology) - this is written to the
1861	! !! history files
1862	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
1863	! !! for crops
1864	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
1865	! !! C (for phenology)
1866	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
1867	! !! leaves were lost (for phenology)
1868
1869	!! 0.4 Local variables
1870
1871	INTEGER(i_std) :: iele !! Indeces(unitless)
1872	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
1873	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements) :: litter_old !! metabolic and structural litter, above and
1874	!! below ground @tex ($gC m^{-2}$) @endtex
1875	REAL(r_std) :: veget_old,veget_total
1876
1877
1878	! Back up some variables in case they're needed later
1879	litter_old(:,:,:,:,:) = litter(:,:,:,:,:)
1880
1881	!! General idea
1882	! The established proxy vegetation has a fraction of 'veget_max_pro'; the
1883	! existing iPFT has a fraction of veget_max(ipts,ipft).
1884	! Suppose we want to merge a scalar variable B, the value of B after merging
1885	! is (BiVi+BjVj)/(Vi+Vj), where Vi is the original veget_max, Vj is the
1886	! incoming veget_max. Note that in case Vi=0, this equation remains solid,
1887	! i.e. the veget_max after merging is Vj and B after merging is Bj. In other
1888	! words, the proxy vegetation "fills" up the empty niche of iPFT.
1889	! Also note that for many scalar variables our input value is Bj*Vj, which
1890	! is accumulated from multiple ancestor PFTs.
1891	veget_old = veget_max(ipts,ipft)
1892	veget_total = veget_old+veget_max_pro
1893
1894	!! Different ways of handling merging depending on nature of variables:
1895
1896	!! 1. Area-based scalar variables, use the equation above
1897	! biomass,carbon, litter, bm_to_litter, co2_to_bm, ind,
1898	! lm_lastyearmax, npp_longterm, lm_lastyearmax,
1899	! lignin_struc (ratio variable depending on area-based variable)
1900
1901	!! 2. Variables are tentatively handled like area-based variables:
1902	! leaf_frac, leaf_age,
1903
1904	!! 3. Variables that are overwritten by the newly initialized PFT:
1905	! PFTpresent, senescence
1906
1907	!! 4. Variables whose operation is uncertain and are not handled currently:
1908	! when_growthinit :: how many days ago was the beginning of the growing season (days)
1909	! gdd_from_growthinit :: growing degree days, since growthinit
1910	! gdd_midwinter, time_hum_min, gdd_m5_dormance, ncd_dormance,
1911	! moiavail_month, moiavail_week, ngd_minus5
1912
1913	!! 5. Variables that concern with short-term fluxes that do not apply in
1914	! this case:
1915	! gpp_daily, npp_daily etc.
1916
1917	! Add the coming veget_max_pro into existing veget_max
1918	veget_max(ipts,ipft) = veget_total
1919
1920	! Merge scalar variables which are defined on area basis
1921	carbon(ipts,:,ipft) = (veget_old * carbon(ipts,:,ipft) + &
1922	carbon_pro(:))/veget_total
1923	deepC_a(ipts,:,ipft) = (veget_old * deepC_a(ipts,:,ipft) + &
1924	deepC_a_pro(:))/veget_total
1925	deepC_s(ipts,:,ipft) = (veget_old * deepC_s(ipts,:,ipft) + &
1926	deepC_s_pro(:))/veget_total
1927	deepC_p(ipts,:,ipft) = (veget_old * deepC_p(ipts,:,ipft) + &
1928	deepC_p_pro(:))/veget_total
1929	litter(ipts,:,ipft,:,:) = (veget_old * litter(ipts,:,ipft,:,:) + &
1930	litter_pro(:,:,:))/veget_total
1931	fuel_1hr(ipts,ipft,:,:) = (veget_old * fuel_1hr(ipts,ipft,:,:) + &
1932	fuel_1hr_pro(:,:))/veget_total
1933	fuel_10hr(ipts,ipft,:,:) = (veget_old * fuel_10hr(ipts,ipft,:,:) + &
1934	fuel_10hr_pro(:,:))/veget_total
1935	fuel_100hr(ipts,ipft,:,:) = (veget_old * fuel_100hr(ipts,ipft,:,:) + &
1936	fuel_100hr_pro(:,:))/veget_total
1937	fuel_1000hr(ipts,ipft,:,:) = (veget_old * fuel_1000hr(ipts,ipft,:,:) + &
1938	fuel_1000hr_pro(:,:))/veget_total
1939
1940	WHERE (litter(ipts,istructural,ipft,:,icarbon) .GT. min_stomate)
1941	lignin_struc(ipts,ipft,:) = (veget_oldlitter_old(ipts,istructural,ipft,:,icarbon) &
1942	lignin_struc(ipts,ipft,:) + litter_pro(istructural,:,icarbon)* &
1943	lignin_struc_pro(:))/(veget_total*litter(ipts,istructural,ipft,:,icarbon))
1944	ENDWHERE
1945	bm_to_litter(ipts,ipft,:,:) = (veget_old * bm_to_litter(ipts,ipft,:,:) + &
1946	bm_to_litter_pro(:,:))/veget_total
1947
1948	biomass(ipts,ipft,:,:) = (biomass(ipts,ipft,:,:)*veget_old + &
1949	biomass_pro(:,:))/veget_total
1950	co2_to_bm(ipts,ipft) = (veget_old*co2_to_bm(ipts,ipft) + &
1951	co2_to_bm_pro)/veget_total
1952	ind(ipts,ipft) = (ind(ipts,ipft)*veget_old + ind_pro)/veget_total
1953	lm_lastyearmax(ipts,ipft) = (lm_lastyearmax(ipts,ipft)*veget_old + &
1954	lm_lastyearmax_pro)/veget_total
1955	npp_longterm(ipts,ipft) = (veget_old * npp_longterm(ipts,ipft) + &
1956	npp_longterm_pro)/veget_total
1957
1958	!CHECK: Here follows the original idea in DOFOCO, more strictly,
1959	! leas mass should be considered together. The same also applies on
1960	! leaf age.
1961	leaf_frac(ipts,ipft,:) = (leaf_frac(ipts,ipft,:)*veget_old + &
1962	leaf_frac_pro(:))/veget_total
1963	leaf_age(ipts,ipft,:) = (leaf_age(ipts,ipft,:)*veget_old + &
1964	leaf_age_pro(:))/veget_total
1965	age(ipts,ipft) = (veget_old * age(ipts,ipft) + &
1966	age_pro)/veget_total
1967
1968	! Everywhere deals with the migration of vegetation. Copy the
1969	! status of the most migrated vegetation for the whole PFT
1970	everywhere(ipts,ipft) = MAX(everywhere(ipts,ipft), everywhere_pro)
1971
1972	! Overwrite the original variables with that from newly initialized
1973	! proxy PFT
1974	PFTpresent(ipts,ipft) = PFTpresent_pro
1975	senescence(ipts,ipft) = senescence_pro
1976
1977	! This is to close carbon loop when writing history variables.
1978	gpp_daily(ipts,ipft) = (veget_old * gpp_daily(ipts,ipft) + &
1979	gpp_daily_pro)/veget_total
1980	npp_daily(ipts,ipft) = (veget_old * npp_daily(ipts,ipft) + &
1981	npp_daily_pro)/veget_total
1982	resp_maint(ipts,ipft) = (veget_old * resp_maint(ipts,ipft) + &
1983	resp_maint_pro)/veget_total
1984	resp_growth(ipts,ipft) = (veget_old * resp_growth(ipts,ipft) + &
1985	resp_growth_pro)/veget_total
1986	resp_hetero(ipts,ipft) = (veget_old * resp_hetero(ipts,ipft) + &
1987	resp_hetero_pro)/veget_total
1988	co2_fire(ipts,ipft) = (veget_old * co2_fire(ipts,ipft) + &
1989	co2_fire_pro)/veget_total
1990
1991	! Phenology- or time-related variables will be copied from original values if
1992	! there is already youngest-age-class PFT there, otherwise they're left
1993	! untouched, because 1. to initiliaze all new PFTs here is wrong and
1994	! phenology is not explicitly considered, so we cannot assign a value
1995	! to these variables. 2. We assume they will be correctly filled if
1996	! other variables are in place (e.g., non-zero leaf mass will lead to
1997	! onset of growing season). In this case, merging a newly initialized PFT
1998	! to an existing one is not the same as merging PFTs when they grow
1999	! old enough to exceed thresholds.
2000
2001	! gpp_week(ipts,ipft) = (veget_old * gpp_week(ipts,ipft) + &
2002	! gpp_week_pro)/veget_total
2003	! when_growthinit(ipts,ipft) = (veget_old * when_growthinit(ipts,ipft) + &
2004	! when_growthinit_pro)/veget_total
2005	! gdd_from_growthinit(ipts,ipft) = (veget_old * gdd_from_growthinit(ipts,ipft) + &
2006	! gdd_from_growthinit_pro)/veget_total
2007	! gdd_midwinter(ipts,ipft) = (veget_old * gdd_midwinter(ipts,ipft) + &
2008	! gdd_midwinter_pro)/veget_total
2009	! time_hum_min(ipts,ipft) = (veget_old * time_hum_min(ipts,ipft) + &
2010	! time_hum_min_pro)/veget_total
2011	! gdd_m5_dormance(ipts,ipft) = (veget_old * gdd_m5_dormance(ipts,ipft) + &
2012	! gdd_m5_dormance_pro)/veget_total
2013	! ncd_dormance(ipts,ipft) = (veget_old * ncd_dormance(ipts,ipft) + &
2014	! ncd_dormance_pro)/veget_total
2015	! moiavail_month(ipts,ipft) = (veget_old * moiavail_month(ipts,ipft) + &
2016	! moiavail_month_pro)/veget_total
2017	! moiavail_week(ipts,ipft) = (veget_old * moiavail_week(ipts,ipft) + &
2018	! moiavail_week_pro)/veget_total
2019	! ngd_minus5(ipts,ipft) = (veget_old * ngd_minus5(ipts,ipft) + &
2020	! ngd_minus5_pro)/veget_total
2021
2022
2023	END SUBROUTINE add_incoming_proxy_pft
2024
2025
2026	! ================================================================================================================================
2027	!! SUBROUTINE : empty_pft
2028	!!
2029	!>\BRIEF : Empty a PFT when,
2030	!! - it is exhausted because of land cover change.
2031	!! - it moves to the next age class
2032	!! \n
2033	!_ ================================================================================================================================
2034	SUBROUTINE empty_pft(ipts, ivm, veget_max, biomass, ind, &
2035	carbon, litter, lignin_struc, bm_to_litter, &
2036	deepC_a, deepC_s, deepC_p, &
2037	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
2038	gpp_daily, npp_daily, gpp_week, npp_longterm, &
2039	co2_to_bm, resp_maint, resp_growth, resp_hetero, &
2040	lm_lastyearmax, leaf_frac, leaf_age, age, &
2041	everywhere, PFTpresent, when_growthinit, &
2042	senescence, gdd_from_growthinit, gdd_midwinter, &
2043	time_hum_min, gdd_m5_dormance, ncd_dormance, &
2044	moiavail_month, moiavail_week, ngd_minus5)
2045
2046	IMPLICIT NONE
2047
2048	!! 0.1 Input variables
2049	INTEGER, INTENT(in) :: ipts !! index for grid cell
2050	INTEGER, INTENT(in) :: ivm !! index for pft
2051
2052	!! 0.2 Output variables
2053
2054	!! 0.3 Modified variables
2055
2056	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
2057	!! May sum to
2058	!! less than unity if the pixel has
2059	!! nobio area. (unitless, 0-1)
2060	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
2061	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
2062	!! @tex $(m^{-2})$ @endtex
2063	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
2064	!! @tex ($gC m^{-2}$) @endtex
2065	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
2066	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
2067	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
2068	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
2069	!! below ground @tex ($gC m^{-2}$) @endtex
2070	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
2071	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
2072	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
2073	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
2074	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
2075	!! above and below ground
2076	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
2077	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2078	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
2079	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2080	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
2081	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2082	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
2083	!! @tex $(gC m^{-2} day^{-1})$ @endtex
2084	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
2085	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
2086	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
2087	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
2088	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2089	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
2090	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2091	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
2092	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2093	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
2094	!! @tex ($gC m^{-2}$) @endtex
2095	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
2096	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
2097	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
2098	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
2099	!! very localized (after its introduction) (?)
2100	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
2101	!! each pixel
2102	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
2103	!! the growing season (days)
2104	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
2105	!! for deciduous trees)
2106	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
2107	!! for crops
2108	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
2109	!! (for phenology) - this is written to the
2110	!! history files
2111	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
2112	!! availability (days)
2113	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
2114	!! C (for phenology)
2115	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
2116	!! leaves were lost (for phenology)
2117	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
2118	!! unitless)
2119	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
2120	!! (0 to 1, unitless)
2121	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
2122	!! -5 deg C (for phenology)
2123
2124	!! 0.4 Local variables
2125	INTEGER(i_std) :: iele !! Indeces(unitless)
2126	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
2127
2128	veget_max(ipts,ivm) = zero
2129	ind(ipts,ivm) = zero
2130	biomass(ipts,ivm,:,:) = zero
2131	litter(ipts,:,ivm,:,:) = zero
2132	fuel_1hr(ipts,ivm,:,:) = zero
2133	fuel_10hr(ipts,ivm,:,:) = zero
2134	fuel_100hr(ipts,ivm,:,:) = zero
2135	fuel_1000hr(ipts,ivm,:,:) = zero
2136	carbon(ipts,:,ivm) = zero
2137	deepC_a(ipts,:,ivm) = zero
2138	deepC_s(ipts,:,ivm) = zero
2139	deepC_p(ipts,:,ivm) = zero
2140	bm_to_litter(ipts,ivm,:,:) = zero
2141	DO ilev=1,nlevs
2142	lignin_struc(ipts,ivm,ilev) = zero
2143	ENDDO
2144	npp_longterm(ipts,ivm) = zero
2145	gpp_daily(ipts,ivm) = zero
2146	gpp_week(ipts,ivm) = zero
2147	resp_maint(ipts,ivm) = zero
2148	resp_growth(ipts,ivm) = zero
2149	resp_hetero(ipts,ivm) = zero
2150	npp_daily(ipts,ivm) = zero
2151	co2_to_bm(ipts,ivm) = zero
2152	lm_lastyearmax(ipts,ivm) = zero
2153	age(ipts,ivm) = zero
2154	leaf_frac(ipts,ivm,:) = zero
2155	leaf_age(ipts,ivm,:) = zero
2156	everywhere(ipts,ivm) = zero
2157	when_growthinit(ipts,ivm) = zero
2158	gdd_from_growthinit(ipts,ivm) = zero
2159	gdd_midwinter(ipts,ivm) = zero
2160	time_hum_min(ipts,ivm) = zero
2161	gdd_m5_dormance(ipts,ivm) = zero
2162	ncd_dormance(ipts,ivm) = zero
2163	moiavail_month(ipts,ivm) = zero
2164	moiavail_week(ipts,ivm) = zero
2165	ngd_minus5(ipts,ivm) = zero
2166	PFTpresent(ipts,ivm) = .FALSE.
2167	senescence(ipts,ivm) = .FALSE.
2168
2169	END SUBROUTINE empty_pft
2170
2171	! ================================================================================================================================
2172	!! SUBROUTINE : gross_lcc_firstday
2173	!!
2174	!>\BRIEF : When necessary, adjust input glcc matrix, and allocate it
2175	!! into different contributing age classes and receiving
2176	!! youngest age classes.
2177	!! \n
2178	!_ ================================================================================================================================
2179
2180	! Note: it has this name because this subroutine will also be called
2181	! the first day of each year to precalculate the forest loss for the
2182	! deforestation fire module.
2183	SUBROUTINE gross_glcc_firstday_fh(npts,veget_max_org,harvest_matrix, &
2184	glccSecondShift,glccPrimaryShift,glccNetLCC,&
2185	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
2186	Deficit_pf2yf_final, Deficit_sf2yf_final, &
2187	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
2188
2189	IMPLICIT NONE
2190
2191	!! 0.1 Input variables
2192
2193	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
2194	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_org !! "maximal" coverage fraction of a PFT on the ground
2195	!! May sum to
2196	!! less than unity if the pixel has
2197	!! nobio area. (unitless, 0-1)
2198	REAL(r_std), DIMENSION(npts,12),INTENT(in) :: harvest_matrix !!
2199	!!
2200	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccSecondShift !! the land-cover-change (LCC) matrix in case a gross LCC is
2201	!! used.
2202	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccPrimaryShift !! the land-cover-change (LCC) matrix in case a gross LCC is
2203	!! used.
2204	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
2205	!! used.
2206
2207	!! 0.2 Output variables
2208	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
2209	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
2210	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
2211	!! after considering the consistency between presribed
2212	!! glcc matrix and existing vegetation fractions.
2213	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
2214	!! there are not enough fractions in the source PFTs
2215	!! /vegetations to target PFTs/vegetations. I.e., these
2216	!! fraction transfers are presribed in LCC matrix but
2217	!! not realized.
2218	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Deficit_pf2yf_final !!
2219	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Deficit_sf2yf_final !!
2220	REAL(r_std), DIMENSION(npts), INTENT(inout) :: pf2yf_compen_sf2yf !!
2221	REAL(r_std), DIMENSION(npts), INTENT(inout) :: sf2yf_compen_pf2yf !!
2222
2223
2224	!! 0.3 Modified variables
2225
2226	!! 0.4 Local variables
2227	REAL(r_std), DIMENSION(npts,nvmap) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
2228	REAL(r_std), DIMENSION(npts,nagec_tree) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
2229	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_grass !! fraction of grass age-class groups, in sequence of old->young
2230	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_pasture !! fraction of pasture age-class groups, in sequence of old->young
2231	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_crop !! fraction of crop age-class groups, in sequence of old->young
2232
2233
2234	REAL(r_std), DIMENSION(npts,4) :: veget_4veg !! "maximal" coverage fraction of a PFT on the ground
2235	REAL(r_std), DIMENSION(npts) :: veget_tree !! "maximal" coverage fraction of a PFT on the ground
2236	REAL(r_std), DIMENSION(npts) :: veget_grass !! "maximal" coverage fraction of a PFT on the ground
2237	REAL(r_std), DIMENSION(npts) :: veget_pasture !! "maximal" coverage fraction of a PFT on the ground
2238	REAL(r_std), DIMENSION(npts) :: veget_crop !! "maximal" coverage fraction of a PFT on the ground
2239
2240	REAL(r_std), DIMENSION(npts,nvm) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
2241	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_tmp !! "maximal" coverage fraction of a PFT on the ground
2242	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_old !! "maximal" coverage fraction of a PFT on the ground
2243	REAL(r_std), DIMENSION(npts,nvm) :: glcc_pft_tmp !! Loss of fraction in each PFT
2244
2245	! Different indexes for convenient local uses
2246	! We define the rules for gross land cover change matrix:
2247	! 1 forest->grass
2248	! 2 forest->pasture
2249	! 3 forest->crop
2250	! 4 grass->forest
2251	! 5 grass->pasture
2252	! 6 grass->crop
2253	! 7 pasture->forest
2254	! 8 pasture->grass
2255	! 9 pasture->crop
2256	! 10 crop->forest
2257	! 11 crop->grass
2258	! 12 crop->pasture
2259	INTEGER :: f2g=1, f2p=2, f2c=3
2260	INTEGER :: g2f=4, g2p=5, g2c=6, p2f=7, p2g=8, p2c=9, c2f=10, c2g=11, c2p=12
2261
2262	INTEGER, ALLOCATABLE :: indall_tree(:) !! Indices for all tree PFTs
2263	INTEGER, ALLOCATABLE :: indold_tree(:) !! Indices for old tree cohort only
2264	INTEGER, ALLOCATABLE :: indagec_tree(:,:) !! Indices for secondary tree cohorts,
2265	!! note the sequence is old->young.
2266	INTEGER, ALLOCATABLE :: indall_grass(:) !! Indices for all grass PFTs
2267	INTEGER, ALLOCATABLE :: indold_grass(:) !! Indices for old grasses only
2268	INTEGER, ALLOCATABLE :: indagec_grass(:,:) !! Indices for secondary grass cohorts
2269	!! note the sequence is old->young.
2270	INTEGER, ALLOCATABLE :: indall_pasture(:) !! Indices for all pasture PFTs
2271	INTEGER, ALLOCATABLE :: indold_pasture(:) !! Indices for old pasture only
2272	INTEGER, ALLOCATABLE :: indagec_pasture(:,:) !! Indices for secondary pasture cohorts
2273	!! note the sequence is old->young.
2274	INTEGER, ALLOCATABLE :: indall_crop(:) !! Indices for all crop PFTs
2275	INTEGER, ALLOCATABLE :: indold_crop(:) !! Indices for old crops only
2276	INTEGER, ALLOCATABLE :: indagec_crop(:,:) !! Indices for secondary crop cohorts
2277	!! note the sequence is old->young.
2278	INTEGER :: num_tree_sinagec,num_tree_mulagec,num_grass_sinagec,num_grass_mulagec, &
2279	num_pasture_sinagec,num_pasture_mulagec,num_crop_sinagec,num_crop_mulagec, &
2280	itree,itree2,igrass,igrass2,ipasture,ipasture2,icrop,icrop2,pf2yf,sf2yf
2281	INTEGER :: i,j,ivma,staind,endind,ivm
2282
2283
2284	REAL(r_std), DIMENSION(npts,12) :: glccDef !! Gross LCC deficit, negative values mean that there
2285	!! are not enough fractions in the source vegetations
2286	!! to the target ones as presribed by the LCC matrix.
2287	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf !!
2288	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf !!
2289	REAL(r_std), DIMENSION(npts) :: Surplus_pf2yf !!
2290	REAL(r_std), DIMENSION(npts) :: Surplus_sf2yf !!
2291	REAL(r_std), DIMENSION(npts,12) :: FHmatrix_remainA !!
2292	REAL(r_std), DIMENSION(npts,12) :: FHmatrix_remainB !!
2293	REAL(r_std), DIMENSION(npts,12) :: glccRemain !!
2294	REAL(r_std), DIMENSION(npts,12) :: glccSecondShift_remain !!
2295	REAL(r_std), DIMENSION(npts,2) :: vegagec_tree_twocl !! Forest fraction in two big classes: the oldest and other
2296	!! age classes.
2297
2298	INTEGER :: ipts,IndStart_f,IndEnd_f
2299
2300
2301	!! 1. We first build all different indices that we are going to use
2302	!! in handling the PFT exchanges, three types of indices are built:
2303	!! - for all age classes
2304	!! - include only oldest age classes
2305	!! - include all age classes excpet the oldest ones
2306	! We have to build these indices because we would like to extract from
2307	! donating PFTs in the sequnce of old->young age classes, and add in the
2308	! receving PFTs only in the youngest-age-class PFTs. These indicies allow
2309	! us to know where the different age classes are.
2310
2311	num_tree_sinagec=0 ! number of tree PFTs with only one single age class
2312	! considered as the oldest age class
2313	num_tree_mulagec=0 ! number of tree PFTs having multiple age classes
2314	num_grass_sinagec=0
2315	num_grass_mulagec=0
2316	num_pasture_sinagec=0
2317	num_pasture_mulagec=0
2318	num_crop_sinagec=0
2319	num_crop_mulagec=0
2320
2321	!! 1.1 Calculate the number of PFTs for different MTCs and allocate
2322	!! the old and all indices arrays.
2323
2324	! [Note here the sequence to identify tree,pasture,grass,crop] is
2325	! critical. The similar sequence is used in the subroutine "calc_cover".
2326	! Do not forget to change the sequence there if you modify here.
2327	DO ivma =2,nvmap
2328	staind=start_index(ivma)
2329	IF (nagec_pft(ivma)==1) THEN
2330	IF (is_tree(staind)) THEN
2331	num_tree_sinagec = num_tree_sinagec+1
2332	ELSE IF (is_grassland_manag(staind)) THEN
2333	num_pasture_sinagec = num_pasture_sinagec+1
2334	ELSE IF (natural(staind)) THEN
2335	num_grass_sinagec = num_grass_sinagec+1
2336	ELSE
2337	num_crop_sinagec = num_crop_sinagec+1
2338	ENDIF
2339
2340	ELSE
2341	IF (is_tree(staind)) THEN
2342	num_tree_mulagec = num_tree_mulagec+1
2343	ELSE IF (is_grassland_manag(staind)) THEN
2344	num_pasture_mulagec = num_pasture_mulagec+1
2345	ELSE IF (natural(staind)) THEN
2346	num_grass_mulagec = num_grass_mulagec+1
2347	ELSE
2348	num_crop_mulagec = num_crop_mulagec+1
2349	ENDIF
2350	ENDIF
2351	ENDDO
2352
2353	!! Allocate index array
2354	! allocate all index
2355	ALLOCATE(indall_tree(num_tree_sinagec+num_tree_mulagec*nagec_tree))
2356	ALLOCATE(indall_grass(num_grass_sinagec+num_grass_mulagec*nagec_herb))
2357	ALLOCATE(indall_pasture(num_pasture_sinagec+num_pasture_mulagec*nagec_herb))
2358	ALLOCATE(indall_crop(num_crop_sinagec+num_crop_mulagec*nagec_herb))
2359
2360	! allocate old-ageclass index
2361	ALLOCATE(indold_tree(num_tree_sinagec+num_tree_mulagec))
2362	ALLOCATE(indold_grass(num_grass_sinagec+num_grass_mulagec))
2363	ALLOCATE(indold_pasture(num_pasture_sinagec+num_pasture_mulagec))
2364	ALLOCATE(indold_crop(num_crop_sinagec+num_crop_mulagec))
2365
2366	!! 1.2 Fill the oldest-age-class and all index arrays
2367	itree=0
2368	igrass=0
2369	ipasture=0
2370	icrop=0
2371	itree2=1
2372	igrass2=1
2373	ipasture2=1
2374	icrop2=1
2375	DO ivma =2,nvmap
2376	staind=start_index(ivma)
2377	IF (is_tree(staind)) THEN
2378	itree=itree+1
2379	indold_tree(itree) = staind+nagec_pft(ivma)-1
2380	DO j = 0,nagec_pft(ivma)-1
2381	indall_tree(itree2+j) = staind+j
2382	ENDDO
2383	itree2=itree2+nagec_pft(ivma)
2384	ELSE IF (natural(staind) .AND. .NOT. is_grassland_manag(staind)) THEN
2385	igrass=igrass+1
2386	indold_grass(igrass) = staind+nagec_pft(ivma)-1
2387	DO j = 0,nagec_pft(ivma)-1
2388	indall_grass(igrass2+j) = staind+j
2389	ENDDO
2390	igrass2=igrass2+nagec_pft(ivma)
2391	ELSE IF (is_grassland_manag(staind)) THEN
2392	ipasture = ipasture+1
2393	indold_pasture(ipasture) = staind+nagec_pft(ivma)-1
2394	DO j = 0,nagec_pft(ivma)-1
2395	indall_pasture(ipasture2+j) = staind+j
2396	ENDDO
2397	ipasture2=ipasture2+nagec_pft(ivma)
2398	ELSE
2399	icrop = icrop+1
2400	indold_crop(icrop) = staind+nagec_pft(ivma)-1
2401	DO j = 0,nagec_pft(ivma)-1
2402	indall_crop(icrop2+j) = staind+j
2403	ENDDO
2404	icrop2=icrop2+nagec_pft(ivma)
2405	ENDIF
2406	ENDDO
2407
2408	!! 1.3 Allocate and fill other age class index
2409
2410	! [chaoyuejoy@gmail.com 2015-08-05]
2411	! note that we treat the case of (num_tree_mulagec==0) differently. In this
2412	! case there is no distinction of age groups among tree PFTs. But we still
2413	! we want to use the "gross_lcchange" subroutine. In this case we consider
2414	! them as having a single age group. In the subroutines
2415	! of "type_conversion" and "cross_give_receive", only the youngest-age-group
2416	! PFTs of a given MTC or vegetation type could receive the incoming fractions.
2417	! To be able to handle this case with least amount of code change, we assign the index
2418	! of PFT between youngest and second-oldes (i.e., indagec_tree etc) the same as
2419	! those of oldest tree PFTs (or all tree PFTs because in this cases these two indices
2420	! are identical) . So that this case could be correctly handled in the subrountines
2421	! of "type_conversion" and "cross_give_receive". This treatment allows use
2422	! of gross land cover change subroutine with only one single age class. This single
2423	! age class is "simultanously the oldest and youngest age class". At the same
2424	! time, we also change the num_tree_mulagec as the same of num_crop_sinagec.
2425	! The similar case also applies in grass,pasture and crop.
2426
2427	IF (num_tree_mulagec .EQ. 0) THEN
2428	ALLOCATE(indagec_tree(num_tree_sinagec,1))
2429	indagec_tree(:,1) = indall_tree(:)
2430	num_tree_mulagec = num_tree_sinagec
2431	ELSE
2432	ALLOCATE(indagec_tree(num_tree_mulagec,nagec_tree-1))
2433	END IF
2434
2435	IF (num_grass_mulagec .EQ. 0) THEN
2436	ALLOCATE(indagec_grass(num_grass_sinagec,1))
2437	indagec_grass(:,1) = indall_grass(:)
2438	num_grass_mulagec = num_grass_sinagec
2439	ELSE
2440	ALLOCATE(indagec_grass(num_grass_mulagec,nagec_herb-1))
2441	END IF
2442
2443	IF (num_pasture_mulagec .EQ. 0) THEN
2444	ALLOCATE(indagec_pasture(num_pasture_sinagec,1))
2445	indagec_pasture(:,1) = indall_pasture(:)
2446	num_pasture_mulagec = num_pasture_sinagec
2447	ELSE
2448	ALLOCATE(indagec_pasture(num_pasture_mulagec,nagec_herb-1))
2449	END IF
2450
2451	IF (num_crop_mulagec .EQ. 0) THEN
2452	ALLOCATE(indagec_crop(num_crop_sinagec,1))
2453	indagec_crop(:,1) = indall_crop(:)
2454	num_crop_mulagec = num_crop_sinagec
2455	ELSE
2456	ALLOCATE(indagec_crop(num_crop_mulagec,nagec_herb-1))
2457	END IF
2458
2459	! fill the non-oldest age class index arrays when number of age classes
2460	! is more than 1.
2461	! [chaoyuejoy@gmail.com, 2015-08-05]
2462	! Note the corresponding part of code will be automatically skipped
2463	! when nagec_tree ==1 and/or nagec_herb ==1, i.e., the assginment
2464	! in above codes when original num_*_mulagec variables are zero will be retained.
2465	itree=0
2466	igrass=0
2467	ipasture=0
2468	icrop=0
2469	DO ivma = 2,nvmap
2470	staind=start_index(ivma)
2471	IF (nagec_pft(ivma) > 1) THEN
2472	IF (is_tree(staind)) THEN
2473	itree=itree+1
2474	DO j = 1,nagec_tree-1
2475	indagec_tree(itree,j) = staind+nagec_tree-j-1
2476	ENDDO
2477	ELSE IF (natural(staind) .AND. .NOT. is_grassland_manag(staind)) THEN
2478	igrass=igrass+1
2479	DO j = 1,nagec_herb-1
2480	indagec_grass(igrass,j) = staind+nagec_herb-j-1
2481	ENDDO
2482	ELSE IF (is_grassland_manag(staind)) THEN
2483	ipasture=ipasture+1
2484	DO j = 1,nagec_herb-1
2485	indagec_pasture(ipasture,j) = staind+nagec_herb-j-1
2486	ENDDO
2487	ELSE
2488	icrop=icrop+1
2489	DO j = 1,nagec_herb-1
2490	indagec_crop(icrop,j) = staind+nagec_herb-j-1
2491	ENDDO
2492	ENDIF
2493	ENDIF
2494	ENDDO
2495
2496	!!! Land cover change processes start here !!!
2497	! we make copies of original input veget_max (which is veget_max_org
2498	! in the subroutine parameter list).
2499	! veget_max will be modified through different operations in order to
2500	! check various purposes, e.g., whether input harvest and glcc matrix
2501	! is compatible with existing veget_max and how to allocate it etc.
2502	! veget_max_old will not be modified
2503	veget_max(:,:) = veget_max_org(:,:)
2504	veget_max_old(:,:) = veget_max_org(:,:)
2505
2506	!******************** block to handle forestry harvest **************
2507	!! 2. Handle the forestry harvest process
2508
2509	!! 2.0 Some preparation
2510
2511	pf2yf=1 !primary to young forest conversion because of harvest
2512	sf2yf=2 !old secondary to young forest conversion because of harvest
2513
2514	! Note that Deficit_pf2yf and Deficit_sf2yf are temporary, intermediate
2515	! variables. The final deficits after mutual compensation are stored in
2516	! Deficit_pf2yf_final and Deficit_sf2yf_final.
2517	Deficit_pf2yf(:) = zero
2518	Deficit_sf2yf(:) = zero
2519	Deficit_pf2yf_final(:) = zero
2520	Deficit_sf2yf_final(:) = zero
2521
2522	! Note that both Surplus_pf2yf and Surplus_sf2yf and temporary intermediate
2523	! variables, the final surplus after mutual compensation are not outputed.
2524	Surplus_pf2yf(:) = zero
2525	Surplus_sf2yf(:) = zero
2526
2527	! Note in the naming of pf2yf_compen_sf2yf and sf2yf_compen_pf2yf, active
2528	! tense is used. I.e., pf2yf_compen_sf2yf means the fraction which pf2yf
2529	! compenstates for sf2yf
2530	pf2yf_compen_sf2yf(:) = zero !primary->young conversion that compensates
2531	!the secondary->young conversion because of deficit
2532	!in the latter
2533	sf2yf_compen_pf2yf(:) = zero !seondary->young conversion that compensates
2534	!the primary->young conversion because of the deficit
2535	!in the latter
2536
2537	! we now have to fill the transtion of forest->forest because of harvest
2538	! into our target matrix glcc_pftmtc. Thus we will initiliaze them first.
2539	glcc_pft(:,:) = 0.
2540	glcc_pft_tmp(:,:) = 0.
2541	glcc_pftmtc(:,:,:) = 0.
2542	glccRemain(:,:) = harvest_matrix(:,:)
2543
2544	!! 2.1 Handle secondary forest harvest
2545
2546	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
2547	vegagec_pasture,vegagec_crop)
2548
2549	! Allocate harvest-caused out-going primary and secondary forest fraction
2550	! into different primary and secondary (all other younger age classes) forest PFTs.
2551	! [Note: below we used the tempelate of type_conversion but in fact we need
2552	! only glcc_pft, which means the fraction loss in each PFT. We then need to
2553	! use glcc_pft to fill glcc_pftmtc (our final target matrix), assuming that
2554	! the loss of forest PFT will go to the youngest age class of its forest MTC.
2555	! Thought glcc_pftmtc and glcc_pft_tmp will be automatically filled when
2556	! we use the tempelate type_conversion by calling it as below, however they
2557	! will be re-set to zero when handling shifting LCC in and net LCC in later
2558	! sections.]
2559
2560	!! 2.1.1 Secondary forest harvest within modeled secondary forest age classes.
2561
2562	! We first handle within the secondary forest age classes, in the sequence
2563	! of old->young
2564
2565	IndStart_f = 2 ! note the indecies and vegetfrac for tree age class
2566	! is from old to young, thus index=2 means the
2567	! 2nd oldest age class.
2568	IndEnd_f = nagec_tree-1 ! the 2nd youngest age class.
2569
2570	DO ipts=1,npts
2571	!sf2yf
2572	CALL type_conversion(ipts,sf2yf,harvest_matrix,veget_mtc, &
2573	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec,&
2574	IndEnd_f,nagec_herb, &
2575	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2576	glccRemain, &
2577	.TRUE., iagec_start=IndStart_f)
2578	ENDDO
2579	FHmatrix_remainA(:,:) = glccRemain
2580
2581	!! 2.1.2 Use primary forest harvest to compensate the deficit in secondary
2582	!! forest harvest within secondary forest in the model.
2583
2584	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
2585	vegagec_pasture,vegagec_crop)
2586
2587	! we check whether the required harvest of secondary forest
2588	! is met by the existing secondary forest fractions. Otherwise
2589	! we use the oldest-age-class forest to compenstate it.
2590	DO ipts=1,npts
2591	IF (FHmatrix_remainA(ipts,sf2yf) .GT. zero) THEN
2592	! in this case, the existing secondary forest fraction
2593	! is not enough for secondary forest harvest, we have to
2594	! use primary (oldest age class) foret to compensate it.
2595
2596	IndStart_f = 1 ! Oldest age class
2597	IndEnd_f = 1 ! Oldest age class
2598
2599	!sf2yf
2600	CALL type_conversion(ipts,sf2yf,FHmatrix_remainA,veget_mtc, &
2601	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec,&
2602	IndEnd_f,nagec_herb, &
2603	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2604	glccRemain, &
2605	.TRUE., iagec_start=IndStart_f)
2606
2607	ENDIF
2608	ENDDO
2609	FHmatrix_remainB(:,:) = glccRemain
2610
2611	!! 2.2 Handle primary forest harvest
2612
2613	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
2614	vegagec_pasture,vegagec_crop)
2615
2616	! we check first if there is still deficit in the required secondary
2617	! harvest. If yes, that means all existing forest (except the youngest
2618	! age class) is depleted, thus required primary harvest will be suppressed.
2619	! Otherwise we will treat primary forest harvest starting from modeled
2620	! oldest-age-class forest
2621
2622	DO ipts=1,npts
2623	IF (FHmatrix_remainB(ipts,sf2yf) .GT. min_stomate) THEN
2624	! in this case, all forest fraction is depleted in handling
2625	! required secondary forest harvest. We thus suppress the
2626	! the required primary forest harvest.
2627	Deficit_sf2yf_final(ipts) = -1 * FHmatrix_remainB(ipts,sf2yf)
2628	Deficit_pf2yf_final(ipts) = -1 * FHmatrix_remainB(ipts,pf2yf)
2629
2630
2631	ELSE
2632	! there are still forest can be used for required primary forest harvest.
2633	! we treat primary harvest wihtin the modeled oldest age class.
2634
2635	IndStart_f = 1 ! Oldest age class
2636	IndEnd_f = nagec_tree-1 ! 2nd youngest age class
2637
2638	!pf2yf
2639	CALL type_conversion(ipts,pf2yf,FHmatrix_remainB,veget_mtc, &
2640	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec,&
2641	IndEnd_f,nagec_herb, &
2642	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2643	glccRemain, &
2644	.TRUE., iagec_start=IndStart_f)
2645	ENDIF
2646
2647	IF (glccRemain(ipts,pf2yf) .GT. min_stomate) THEN
2648	Deficit_pf2yf_final(ipts) = -1 * glccRemain(ipts,pf2yf)
2649	ENDIF
2650	ENDDO
2651
2652	! Because we use the container of type_conversion, now the glcc_pft_tmp
2653	! and glcc_pftmtc have wrong information (because harvest loss is assigned
2654	! on the newly created youngest-age-class pasture/crop MTCs). So they have
2655	! to be re-initialized to zero. Only the information in glcc_pft is what
2656	! we need, as explained above.
2657	glcc_pft_tmp(:,:) = 0.
2658	glcc_pftmtc(:,:,:) = 0.
2659	!Here we need to put glcc_pft into glcc_pftmtc for forestry harvest.
2660	!The same MTC will be maintained when forest is harvested.
2661	DO ivm =1,nvm
2662	IF (is_tree(ivm)) THEN
2663	glcc_pftmtc(:,ivm,pft_to_mtc(ivm)) = glcc_pft(:,ivm)
2664	ENDIF
2665	ENDDO
2666	!**************** end block to handle forestry harvest **************
2667
2668	!! 3. Treat secondary-agriculture shifting cultivation transition matrix
2669	!! [The primary-agriculture shifting cultivation will be treated together
2670	!! with the netLCC transitions, with the conversion sequence of oldest->
2671	!! youngest is applied.]
2672	! When we prepare the driving data, secondary-agriculture shifting cultivation
2673	! is intended to include the "constant transitions" over time. Ideally, we
2674	! should start applying this secondary-agriculture shifting cultivation with
2675	! the "secondary forest" in the model. Here we tentatively start with the 3rd
2676	! youngest age class and move to the 2ne youngest age class. But if the prescribed
2677	! transition fraction is not met, we then move further to 4th youngest age class
2678	! and then move to the oldest age class sequentially.
2679
2680	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
2681	vegagec_pasture,vegagec_crop)
2682
2683	!! 3.1 We start treating secondary-agriculture cultivation from the 3rd youngest
2684	!! age class and then move to the younger age class.
2685	! Because it's rather complicated to calculate which transtion fraction between
2686	! which vegetation types should stay in here in case there is deficit occuring
2687	! for the overall donation vegetation type, we will just start from some
2688	! priority and leave the unrealized parts into the latter section.
2689
2690	! For this purpose, we should first make a copy of glccSecondShift into
2691	! glccRemain. glccRemain will tell us the transition fractions that have to
2692	! be treated starting from 3rd oldest age class and moving torward older
2693	! age class.
2694	glccRemain(:,:) = glccSecondShift(:,:)
2695
2696	! Now we will call type_conversion for each of the 12 transitions, starting
2697	! from 2nd age class moving to the youngest age class. We use glccRemain
2698	! to track the transtion fractions we should leave for the second case.
2699	! To make the code more flexible, we will store the start and end indecies
2700	! in variables.
2701
2702	!*[Note: we do above process only for forest now, as we assume the conversion
2703	! of crop/pasture/grass to other types will start always from the oldest
2704	! age class]
2705
2706	IndStart_f = nagec_tree-2 ! note the indecies and vegetfrac for tree age class
2707	! is from old to young, thus nagec_tree-1 means the
2708	! 3rd youngest age class.
2709	IndEnd_f = nagec_tree-2 ! nagec_tree-2: The 3rd youngest age class
2710	! nagec_tree-1: The 2nd youngest age class
2711	! nagec_tree: The youngest age class
2712
2713
2714	DO ipts=1,npts
2715	!f2c
2716	CALL type_conversion(ipts,f2c,glccSecondShift,veget_mtc, &
2717	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
2718	IndEnd_f,nagec_herb, &
2719	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
2720	glccRemain, &
2721	.TRUE., iagec_start=IndStart_f)
2722	!f2p
2723	CALL type_conversion(ipts,f2p,glccSecondShift,veget_mtc, &
2724	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
2725	IndEnd_f,nagec_herb, &
2726	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
2727	glccRemain, &
2728	.TRUE., iagec_start=IndStart_f)
2729	!f2g
2730	CALL type_conversion(ipts,f2g,glccSecondShift,veget_mtc, &
2731	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
2732	IndEnd_f,nagec_herb, &
2733	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2734	glccRemain, &
2735	.TRUE., iagec_start=IndStart_f)
2736	!g2c
2737	CALL type_conversion(ipts,g2c,glccSecondShift,veget_mtc, &
2738	indold_grass,indagec_grass,indagec_crop,num_crop_mulagec, &
2739	nagec_herb,nagec_herb, &
2740	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2741	glccRemain, &
2742	.TRUE.)
2743	!g2p
2744	CALL type_conversion(ipts,g2p,glccSecondShift,veget_mtc, &
2745	indold_grass,indagec_grass,indagec_pasture,num_pasture_mulagec, &
2746	nagec_herb,nagec_herb, &
2747	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2748	glccRemain, &
2749	.TRUE.)
2750	!g2f
2751	CALL type_conversion(ipts,g2f,glccSecondShift,veget_mtc, &
2752	indold_grass,indagec_grass,indagec_tree,num_tree_mulagec, &
2753	nagec_herb,nagec_tree, &
2754	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2755	glccRemain, &
2756	.TRUE.)
2757	!p2c
2758	CALL type_conversion(ipts,p2c,glccSecondShift,veget_mtc, &
2759	indold_pasture,indagec_pasture,indagec_crop,num_crop_mulagec, &
2760	nagec_herb,nagec_herb, &
2761	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2762	glccRemain, &
2763	.TRUE.)
2764	!p2g
2765	CALL type_conversion(ipts,p2g,glccSecondShift,veget_mtc, &
2766	indold_pasture,indagec_pasture,indagec_grass,num_grass_mulagec, &
2767	nagec_herb,nagec_herb, &
2768	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2769	glccRemain, &
2770	.TRUE.)
2771	!p2f
2772	CALL type_conversion(ipts,p2f,glccSecondShift,veget_mtc, &
2773	indold_pasture,indagec_pasture,indagec_tree,num_tree_mulagec, &
2774	nagec_herb,nagec_tree, &
2775	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2776	glccRemain, &
2777	.TRUE.)
2778	!c2p
2779	CALL type_conversion(ipts,c2p,glccSecondShift,veget_mtc, &
2780	indold_crop,indagec_crop,indagec_pasture,num_pasture_mulagec, &
2781	nagec_herb,nagec_herb, &
2782	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2783	glccRemain, &
2784	.TRUE.)
2785	!c2g
2786	CALL type_conversion(ipts,c2g,glccSecondShift,veget_mtc, &
2787	indold_crop,indagec_crop,indagec_grass,num_grass_mulagec, &
2788	nagec_herb,nagec_herb, &
2789	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2790	glccRemain, &
2791	.TRUE.)
2792	!c2f
2793	CALL type_conversion(ipts,c2f,glccSecondShift,veget_mtc, &
2794	indold_crop,indagec_crop,indagec_tree,num_tree_mulagec, &
2795	nagec_herb,nagec_tree, &
2796	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2797	glccRemain, &
2798	.TRUE.)
2799	ENDDO
2800	glccSecondShift_remain(:,:) = glccRemain(:,:)
2801
2802	!! 3.2 We treat the remaing unrealized transtions from forest. Now we will
2803	!! start with the 3rd oldest age class and then move to the oldest age class.
2804
2805	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
2806	vegagec_pasture,vegagec_crop)
2807
2808	IndStart_f = nagec_tree-3 ! note the indecies and vegetfrac for tree age class
2809	! is from old to young, thus nagec_tree-2 means the
2810	! 3rd oldest age class.
2811	IndEnd_f = 1
2812
2813	! we start with the 3rd youngest age class and move up to the oldest age
2814	! class in the sequence of young->old, as indicated by the .FALSE. parameter
2815	! when calling the subroutine type_conversion.
2816	DO ipts=1,npts
2817	!f2c
2818	CALL type_conversion(ipts,f2c,glccSecondShift_remain,veget_mtc, &
2819	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
2820	IndEnd_f,nagec_herb, &
2821	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
2822	glccRemain, &
2823	.FALSE., iagec_start=IndStart_f)
2824	!f2p
2825	CALL type_conversion(ipts,f2p,glccSecondShift_remain,veget_mtc, &
2826	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
2827	IndEnd_f,nagec_herb, &
2828	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
2829	glccRemain, &
2830	.FALSE., iagec_start=IndStart_f)
2831	!f2g
2832	CALL type_conversion(ipts,f2g,glccSecondShift_remain,veget_mtc, &
2833	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
2834	IndEnd_f,nagec_herb, &
2835	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2836	glccRemain, &
2837	.FALSE., iagec_start=IndStart_f)
2838	ENDDO
2839
2840	! we put the remaining glccRemain into the deficit
2841	IncreDeficit(:,:) = -1*glccRemain
2842	!***end block to handle secondary-agriculture shifting cultivation *****
2843
2844
2845	!+++ Code freezing: Compensation among different transition fractions +++
2846	!+++ Description: This block of code and associated subroutines are originally
2847	!+++ developed to make the LCC module compatible with DGVM.
2848	!
2849	!! we copy updated veget_max to veget_max_tmp.
2850	!! The latter will be used to retrieve the values of veget_max after checking
2851	!! the consistency of input glcc with existing vegetation fractions.
2852	!veget_max_tmp(:,:) = veget_max(:,:)
2853
2854	!!************************************************************************!
2855	!!**block to calculate fractions for basic veg types and age classes *!
2856	!! Note:
2857	!! 1. "calc_cover" subroutine does not depend on how many age classes
2858	!! there are in each MTC.
2859	!! 2. Fraction of baresoil is excluded here. This means transformation
2860	!! of baresoil to a vegetated PFT is excluded in gross land cover change.
2861	!veget_mtc(:,:) = 0.
2862	!vegagec_tree(:,:) = 0.
2863	!vegagec_grass(:,:) = 0.
2864	!vegagec_pasture(:,:) = 0.
2865	!vegagec_crop(:,:) = 0.
2866
2867
2868	!CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
2869	! vegagec_pasture,vegagec_crop)
2870
2871	!veget_tree(:) = SUM(vegagec_tree(:,:),DIM=2)
2872	!veget_grass(:) = SUM(vegagec_grass(:,:),DIM=2)
2873	!veget_pasture(:) = SUM(vegagec_pasture(:,:),DIM=2)
2874	!veget_crop(:) = SUM(vegagec_crop(:,:),DIM=2)
2875	!itree=1
2876	!igrass=2
2877	!ipasture=3
2878	!icrop=4
2879	!veget_4veg(:,itree) = veget_tree(:)
2880	!veget_4veg(:,igrass) = veget_grass(:)
2881	!veget_4veg(:,ipasture) = veget_pasture(:)
2882	!veget_4veg(:,icrop) = veget_crop(:)
2883	!!**end block to calculate fractions for basic veg types and age classes *!
2884	!!****************************************************************************!
2885
2886	!!! 3. Decompose the LCC matrix to different PFTs
2887	!!! We do this through several steps:
2888	!! 3.1 Check whether input LCC matrix is feasible with current PFT fractions
2889	!! (i.e., the fractions of forest,grass,pasture and crops)
2890	!! and if not, adjust the transfer matrix by compensating the deficits
2891	!! using the surpluses.
2892	!! 3.2 Allocate the decreasing fractions of tree/grass/pasture/crop to their
2893	!! respective age classes, in the sequences of old->young.
2894	!! 3.3 Allocate the incoming fractions of tree/grass/pasture/crop to their
2895	!! respective youngest age classes. The incoming fractions are distributed
2896	!! according to the existing fractions of youngest-age-class PFTs of the
2897	!! same receiving vegetation type. If none of them exists, the incoming
2898	!! fraction is distributed equally.
2899
2900	!!! 3.1 Adjust LCC matrix if it's not feasible with current PFT fractions
2901
2902	!IncreDeficit(:,:) = 0.
2903	!glccReal(:,:) = 0.
2904	glccDef(:,:) = 0.
2905
2906	!!to crop - sequence: p2c,g2c,f2c
2907	!CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
2908	! p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
2909	! IncreDeficit)
2910
2911	!!to pasture - sequence: g2p,c2p,f2p
2912	!CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
2913	! g2p,igrass,c2p,icrop,f2p,itree,ipasture, &
2914	! IncreDeficit)
2915
2916	!!to grass - sequence: p2g,c2g,f2g
2917	!CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
2918	! p2g,ipasture,c2g,icrop,f2g,itree,igrass, &
2919	! IncreDeficit)
2920
2921	!!to forest - sequence: c2f,p2f,g2f
2922	!CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
2923	! c2f,icrop,p2f,ipasture,g2f,igrass,itree, &
2924	! IncreDeficit)
2925
2926	!!! 3.2 & 3.3 Allocate LCC matrix to different PFTs/age-classes
2927
2928	!! because we use veget_max as a proxy variable and it has been changed
2929	!! when we derive the glccReal, so here we have to recover its original
2930	!! values, which is veget_max_tmp after the forestry harvest.
2931	!veget_max(:,:) = veget_max_tmp(:,:)
2932	!
2933	!+++ end freezing block of code +++
2934
2935
2936	!! 4. Treat the transtions involving the oldest age classes, which include
2937	!! the first-time primary-agriculture cultivation and the net land cover
2938	!! transtions
2939
2940	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
2941	vegagec_pasture,vegagec_crop)
2942
2943
2944	! the variable "glccReal" is originally for storing the realized maxtrix
2945	! after considering the constraining and compensation of existing vegetation
2946	! fractions. But as this case is not allowed at the moment, we will just
2947	! simply put it as the sum of glccPrimaryShift and glccNetLCC
2948	glccReal(:,:) = glccPrimaryShift+glccNetLCC
2949
2950	! We copy the glccReal to glccRemain in order to track the remaining
2951	! prescribed transtion fraction after applying each transition by calling
2952	! the subroutine "type_conversion". For the moment this is mainly to fufill
2953	! the parameter requirement of the type_conversion subroutine.
2954	glccRemain(:,:) = glccReal(:,:)
2955
2956	! We allocate in the sequences of old->young. Within the same age-class
2957	! group, we allocate in proportion with existing PFT fractions.
2958	DO ipts=1,npts
2959	!f2c
2960	CALL type_conversion(ipts,f2c,glccReal,veget_mtc, &
2961	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
2962	nagec_tree,nagec_herb, &
2963	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
2964	glccRemain, &
2965	.TRUE.)
2966	!f2p
2967	CALL type_conversion(ipts,f2p,glccReal,veget_mtc, &
2968	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
2969	nagec_tree,nagec_herb, &
2970	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
2971	glccRemain, &
2972	.TRUE.)
2973	!f2g
2974	CALL type_conversion(ipts,f2g,glccReal,veget_mtc, &
2975	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
2976	nagec_tree,nagec_herb, &
2977	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2978	glccRemain, &
2979	.TRUE.)
2980	!g2c
2981	CALL type_conversion(ipts,g2c,glccReal,veget_mtc, &
2982	indold_grass,indagec_grass,indagec_crop,num_crop_mulagec, &
2983	nagec_herb,nagec_herb, &
2984	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2985	glccRemain, &
2986	.TRUE.)
2987	!g2p
2988	CALL type_conversion(ipts,g2p,glccReal,veget_mtc, &
2989	indold_grass,indagec_grass,indagec_pasture,num_pasture_mulagec, &
2990	nagec_herb,nagec_herb, &
2991	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2992	glccRemain, &
2993	.TRUE.)
2994	!g2f
2995	CALL type_conversion(ipts,g2f,glccReal,veget_mtc, &
2996	indold_grass,indagec_grass,indagec_tree,num_tree_mulagec, &
2997	nagec_herb,nagec_tree, &
2998	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
2999	glccRemain, &
3000	.TRUE.)
3001	!p2c
3002	CALL type_conversion(ipts,p2c,glccReal,veget_mtc, &
3003	indold_pasture,indagec_pasture,indagec_crop,num_crop_mulagec, &
3004	nagec_herb,nagec_herb, &
3005	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
3006	glccRemain, &
3007	.TRUE.)
3008	!p2g
3009	CALL type_conversion(ipts,p2g,glccReal,veget_mtc, &
3010	indold_pasture,indagec_pasture,indagec_grass,num_grass_mulagec, &
3011	nagec_herb,nagec_herb, &
3012	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
3013	glccRemain, &
3014	.TRUE.)
3015	!p2f
3016	CALL type_conversion(ipts,p2f,glccReal,veget_mtc, &
3017	indold_pasture,indagec_pasture,indagec_tree,num_tree_mulagec, &
3018	nagec_herb,nagec_tree, &
3019	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
3020	glccRemain, &
3021	.TRUE.)
3022	!c2p
3023	CALL type_conversion(ipts,c2p,glccReal,veget_mtc, &
3024	indold_crop,indagec_crop,indagec_pasture,num_pasture_mulagec, &
3025	nagec_herb,nagec_herb, &
3026	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
3027	glccRemain, &
3028	.TRUE.)
3029	!c2g
3030	CALL type_conversion(ipts,c2g,glccReal,veget_mtc, &
3031	indold_crop,indagec_crop,indagec_grass,num_grass_mulagec, &
3032	nagec_herb,nagec_herb, &
3033	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
3034	glccRemain, &
3035	.TRUE.)
3036	!c2f
3037	CALL type_conversion(ipts,c2f,glccReal,veget_mtc, &
3038	indold_crop,indagec_crop,indagec_tree,num_tree_mulagec, &
3039	nagec_herb,nagec_tree, &
3040	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
3041	glccRemain, &
3042	.TRUE.)
3043	ENDDO
3044
3045	! Note here IncreDeficit includes the deficit from secondary<->agriculgure shifting
3046	! cultivation and the primary<->agriculture+NetLCC transitions.
3047	IncreDeficit(:,:) = IncreDeficit(:,:) - glccRemain(:,:)
3048
3049	END SUBROUTINE gross_glcc_firstday_fh
3050
3051
3052	! ================================================================================================================================
3053	!! SUBROUTINE : cross_give_receive
3054	!!
3055	!>\BRIEF : Allocate the outgoing and receving fractions in respective
3056	!! PFTs.
3057	!! \n
3058	!! Notes:
3059	!! 1. veget_max is subtracted when fractions are taken out, but newly added
3060	!! fractions in the youngest age class is not added, to avoid this newly
3061	!! created fractions being used again the following transitions. This is
3062	!! is reasonable because the newly created youngest-age-class PFT fractions
3063	!! have nothing but small sapling biomass and it's unreasonable to use it
3064	!! for any further land use conversion activities.
3065	!_ ================================================================================================================================
3066	SUBROUTINE cross_give_receive(ipts,frac_used,veget_mtc, &
3067	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
3068	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
3069
3070
3071	IMPLICIT NONE
3072
3073	!! 0. Input variables
3074	INTEGER, INTENT(in) :: ipts
3075	REAL(r_std), INTENT(in) :: frac_used !! fraction that the giving PFTs are going to collectively give
3076	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
3077	INTEGER, DIMENSION(:), INTENT(in) :: indold_tree !! Indices for PFTs giving out fractions;
3078	!! here use old tree cohort as an example
3079	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_crop !! Indices for secondary basic-vegetation cohorts; The youngest age classes
3080	!! of these vegetations are going to receive fractions.
3081	!! here we use crop cohorts as an example
3082	INTEGER, INTENT(in) :: num_crop_mulagec !! number of crop MTCs with more than one age classes
3083	INTEGER, INTENT(in) :: nagec_receive !! number of age classes in the receiving basic types
3084	!! (i.e., tree, grass, pasture, crop), here we can use crop
3085	!! as an example, nagec_receive=nagec_herb
3086
3087	!! 1. Modified variables
3088	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
3089	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fractions each PFT is going to lose
3090	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
3091	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
3092	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft_tmp !! a temporary variable to hold the fractions each PFT is going to lose
3093
3094	!! Local vriables
3095	INTEGER :: j,ipft, iyoung
3096	REAL(r_std) :: totalveg
3097
3098
3099	! Out final objective is to know glcc_pftmtc, i.e., the fraction from each PFT
3100	! to the youngest age group of each MTC. We separate this task into two steps:
3101	! 1. we allocate the total outgoing fraction into the same age-class PFTs of
3102	! the a basic-vegetation (for example, the same age-calss PFTs of forest);
3103	! 2. we further allocate the outgoing fraction of each age-class PFT to
3104	! the different receiving youngest age-class PFTs of the same basic-vegetation
3105	! type, for example, the youngest age-calss PFTs of cropland.
3106
3107	! glcc_pft_tmp used only as a temporary variable to store the value
3108	glcc_pft_tmp(ipts,indold_tree) = veget_max(ipts,indold_tree)/SUM(veget_max(ipts,indold_tree))*frac_used
3109	glcc_pft(ipts,indold_tree) = glcc_pft(ipts,indold_tree) + glcc_pft_tmp(ipts,indold_tree)
3110	!we have to remove the outgoing fraction from veget_max in order to use this information for next loop
3111	veget_max(ipts,indold_tree) = veget_max(ipts,indold_tree) - glcc_pft_tmp(ipts,indold_tree)
3112
3113	! when receiving basic-vegetation type has a single age group, it will be considered as
3114	! both old and young age group (thus recevie the fraction donation), otherwise the youngest
3115	! age group is always the final element of indagec_crop.
3116	IF (nagec_receive == 1) THEN
3117	iyoung = 1
3118	ELSE
3119	iyoung = nagec_receive - 1
3120	ENDIF
3121
3122	![20160130 note here totalveg is the total fraction of all existing MTCs
3123	! that are going to recieve newly convervted fractions.]
3124	totalveg = 0.
3125	DO j=1,num_crop_mulagec
3126	totalveg = totalveg + veget_mtc(ipts,agec_group(indagec_crop(j,iyoung)))
3127	ENDDO
3128
3129	IF (totalveg>min_stomate) THEN
3130	DO j=1,num_crop_mulagec
3131	ipft = indagec_crop(j,iyoung)
3132	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) = glcc_pft_tmp(ipts,indold_tree) &
3133	*veget_mtc(ipts,agec_group(ipft))/totalveg
3134	ENDDO
3135	ELSE
3136	DO j=1,num_crop_mulagec
3137	ipft = indagec_crop(j,iyoung)
3138	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) = glcc_pft_tmp(ipts,indold_tree)/num_crop_mulagec
3139	ENDDO
3140	ENDIF
3141
3142	END SUBROUTINE cross_give_receive
3143
3144	! ================================================================================================================================
3145	!! SUBROUTINE : type_conversion
3146	!>\BRIEF : Allocate outgoing into different age classes and incoming into
3147	!! yongest age-class of receiving MTCs.
3148	!!
3149	!! REMARK : The current dummy variables give an example of converting forests
3150	!! to crops.
3151	!! \n
3152	!_ ================================================================================================================================
3153	SUBROUTINE type_conversion(ipts,f2c,glccReal,veget_mtc, &
3154	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
3155	iagec_tree_end,nagec_receive, &
3156	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
3157	glccRemain, &
3158	old_to_young, iagec_start)
3159
3160	IMPLICIT NONE
3161
3162	!! Input variables
3163	INTEGER, INTENT(in) :: ipts,f2c
3164	REAL(r_std), DIMENSION(:,:), INTENT(in) :: glccReal !! The "real" glcc matrix that we apply in the model
3165	!! after considering the consistency between presribed
3166	!! glcc matrix and existing vegetation fractions.
3167	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
3168	INTEGER, DIMENSION(:), INTENT(in) :: indold_tree !! Indices for PFTs giving out fractions;
3169	!! here use old tree cohort as an example
3170	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_tree !! Indices for PFTs giving out fractions;
3171	!! here use old tree cohort as an example
3172	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_crop !! Indices for secondary basic-vegetation cohorts; The youngest age classes
3173	!! of these vegetations are going to receive fractions.
3174	!! here we use crop cohorts as an example
3175	INTEGER, INTENT(in) :: num_crop_mulagec !! number of crop MTCs with more than one age classes
3176	INTEGER, INTENT(in) :: iagec_tree_end !! End index of age classes in the giving basic types
3177	!! (i.e., tree, grass, pasture, crop)
3178	INTEGER, INTENT(in) :: nagec_receive !! number of age classes in the receiving basic types
3179	!! (i.e., tree, grass, pasture, crop), here we can use crop
3180	!! as an example, nagec=nagec_herb
3181	LOGICAL, INTENT(in) :: old_to_young !! an logical variable indicating whether we should handle donation
3182	!! vegetation in a sequence of old->young or young->old. TRUE is for
3183	!! old->young.
3184	INTEGER, OPTIONAL, INTENT(in) :: iagec_start !! starting index for iagec, this is added in order to handle
3185	!! the case of secondary forest harvest.
3186
3187	!! 1. Modified variables
3188	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
3189	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
3190	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fractions each PFT is going to lose
3191	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
3192	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft_tmp !! Loss of fraction in each PFT
3193	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glccRemain !! The remaining glcc matrix after applying the conversion. I.e., it will
3194	!! record the remaining unrealized transition fraction in case the donation
3195	!! vegetation is not enough compared with prescribed transition fraction.
3196	!! This variable should be initialized the same as glccReal before it's fed
3197	!! into this function.
3198
3199	!! Local vriables
3200	INTEGER :: j,iagec,iagec_start_proxy
3201	REAL(r_std) :: frac_begin,frac_used
3202	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
3203	IF (.NOT. PRESENT(iagec_start)) THEN
3204	iagec_start_proxy=1
3205	ELSE
3206	iagec_start_proxy=iagec_start
3207	ENDIF
3208
3209	! This subroutine handles the conversion from one basic-vegetation type
3210	! to another, by calling the subroutine cross_give_receive, which handles
3211	! allocation of giving-receiving fraction among the giving age classes
3212	! and receiving basic-vegetation young age classes.
3213	! We allocate in the sequences of old->young. Within the same age-class
3214	! group, we allocate in proportion with existing PFT fractions. The same
3215	! also applies in the receiving youngest-age-class PFTs, i.e., the receiving
3216	! total fraction is allocated according to existing fractions of
3217	! MTCs of the same basic vegetation type, otherwise it will be equally
3218	! distributed.
3219
3220	frac_begin = glccReal(ipts,f2c)
3221	!DO WHILE (frac_begin>min_stomate)
3222	IF (old_to_young) THEN
3223	! note that both indagec_tree and vegagec_tree are in sequence of old->young
3224	! thus iagec_start_proxy must be smaller than iagec_tree_end
3225	DO iagec=iagec_start_proxy,iagec_tree_end,1
3226	IF (vegagec_tree(ipts,iagec)>frac_begin) THEN
3227	frac_used = frac_begin
3228	ELSE IF (vegagec_tree(ipts,iagec)>min_stomate) THEN
3229	frac_used = vegagec_tree(ipts,iagec)
3230	ELSE
3231	frac_used = 0.
3232	ENDIF
3233
3234	IF (frac_used>min_stomate) THEN
3235	IF (iagec==1) THEN
3236	! Note that vegagec_tree is fractions of tree age-class groups in the
3237	! the sequence of old->young, so iagec==1 means that we're handling
3238	! first the oldest-age-group tree PFTs.
3239	CALL cross_give_receive(ipts,frac_used,veget_mtc, &
3240	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
3241	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
3242	ELSE
3243	! Note also the sequence of indagec_tree is from old->young, so by
3244	! increasing iagec, we're handling progressively the old to young
3245	! tree age-class PFTs.
3246	CALL cross_give_receive(ipts,frac_used,veget_mtc, &
3247	indagec_tree(:,iagec-1),indagec_crop,nagec_receive,num_crop_mulagec, &
3248	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
3249	ENDIF
3250	frac_begin = frac_begin-frac_used
3251	vegagec_tree(ipts,iagec)=vegagec_tree(ipts,iagec)-frac_used
3252	glccRemain(ipts,f2c) = glccRemain(ipts,f2c) - frac_used
3253	ENDIF
3254	ENDDO
3255	ELSE ! in the sequence of young->old
3256	DO iagec=iagec_start_proxy,iagec_tree_end,-1
3257	IF (vegagec_tree(ipts,iagec)>frac_begin) THEN
3258	frac_used = frac_begin
3259	ELSE IF (vegagec_tree(ipts,iagec)>min_stomate) THEN
3260	frac_used = vegagec_tree(ipts,iagec)
3261	ELSE
3262	frac_used = 0.
3263	ENDIF
3264
3265	IF (frac_used>min_stomate) THEN
3266	IF (iagec==1) THEN
3267	! Note that vegagec_tree is fractions of tree age-class groups in the
3268	! the sequence of old->young, so iagec==1 means that we're handling
3269	! first the oldest-age-group tree PFTs.
3270	CALL cross_give_receive(ipts,frac_used,veget_mtc, &
3271	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
3272	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
3273	ELSE
3274	! Note also the sequence of indagec_tree is from old->young, so by
3275	! increasing iagec, we're handling progressively the old to young
3276	! tree age-class PFTs.
3277	CALL cross_give_receive(ipts,frac_used,veget_mtc, &
3278	indagec_tree(:,iagec-1),indagec_crop,nagec_receive,num_crop_mulagec, &
3279	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
3280	ENDIF
3281	frac_begin = frac_begin-frac_used
3282	vegagec_tree(ipts,iagec)=vegagec_tree(ipts,iagec)-frac_used
3283	glccRemain(ipts,f2c) = glccRemain(ipts,f2c) - frac_used
3284	ENDIF
3285	ENDDO
3286	ENDIF
3287	!ENDDO
3288
3289	END SUBROUTINE type_conversion
3290
3291	! ================================================================================================================================
3292	!! SUBROUTINE : calc_cover
3293	!!
3294	!>\BRIEF Calculate coverage fraction for different age classes of forest,
3295	!! grass, pasture and crops and also for each metaclass. Note baresoil is excluded.
3296	!!
3297	!! DESCRIPTION :
3298	!! Note:
3299	!! 1. "calc_cover" subroutine does not depend on how many age classes
3300	!! there are in each MTC.
3301	!! 2. Fraction of baresoil is excluded here. This means transformation
3302	!! of baresoil to a vegetated PFT is excluded in gross land cover change.
3303	!!
3304	!!
3305	!! MAIN OUTPUT VARIABLE(S) :
3306	!!
3307	!! \n
3308	!_ ================================================================================================================================
3309	SUBROUTINE calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
3310	vegagec_pasture,vegagec_crop)
3311
3312
3313	IMPLICIT NONE
3314
3315	!! Input variables
3316	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
3317	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
3318
3319	!! Output variables
3320	REAL(r_std), DIMENSION(npts,nvmap), INTENT(inout) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
3321	REAL(r_std), DIMENSION(npts,nagec_tree), INTENT(inout) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
3322	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_grass !! fraction of grass age-class groups, in sequence of old->young
3323	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_pasture !! fraction of pasture age-class groups, in sequence of old->young
3324	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_crop !! fraction of crop age-class groups, in sequence of old->young
3325
3326	!! Local variables
3327	INTEGER(i_std) :: ivma,staind,endind,j !! indices (unitless)
3328
3329	veget_mtc(:,:) = 0.
3330	vegagec_tree(:,:) = 0.
3331	vegagec_grass(:,:) = 0.
3332	vegagec_pasture(:,:) = 0.
3333	vegagec_crop(:,:) = 0.
3334
3335	! Calculate veget_max for MTCs
3336	DO ivma = 1,nvmap
3337	staind = start_index(ivma)
3338	IF (nagec_pft(ivma) == 1) THEN
3339	veget_mtc(:,ivma) = veget_max(:,staind)
3340	ELSE
3341	veget_mtc(:,ivma) = \
3342	SUM(veget_max(:,staind:staind+nagec_pft(ivma)-1),DIM=2)
3343	ENDIF
3344	ENDDO
3345
3346	! Calculate veget_max for each age class
3347	DO ivma = 2,nvmap !here we start with 2 to exclude baresoil (always PFT1)
3348	staind = start_index(ivma)
3349	endind = staind+nagec_pft(ivma)-1
3350
3351	! Single-age-class MTC goest to oldest age class.
3352	IF (nagec_pft(ivma) == 1) THEN
3353	IF (is_tree(staind)) THEN
3354	vegagec_tree(:,1) = vegagec_tree(:,1)+veget_max(:,staind)
3355	ELSE IF (is_grassland_manag(staind)) THEN
3356	vegagec_pasture(:,1) = vegagec_pasture(:,1)+veget_max(:,staind)
3357	ELSE IF (natural(staind)) THEN
3358	vegagec_grass(:,1) = vegagec_grass(:,1)+veget_max(:,staind)
3359	ELSE
3360	vegagec_crop(:,1) = vegagec_crop(:,1)+veget_max(:,staind)
3361	ENDIF
3362
3363	ELSE
3364	IF (is_tree(staind)) THEN
3365	DO j=1,nagec_tree
3366	vegagec_tree(:,j) = vegagec_tree(:,j)+veget_max(:,endind-j+1)
3367	ENDDO
3368	ELSE IF (is_grassland_manag(staind)) THEN
3369	DO j=1,nagec_herb
3370	vegagec_pasture(:,j) = vegagec_pasture(:,j)+veget_max(:,endind-j+1)
3371	ENDDO
3372	ELSE IF (natural(staind)) THEN
3373	DO j=1,nagec_herb
3374	vegagec_grass(:,j) = vegagec_grass(:,j)+veget_max(:,endind-j+1)
3375	ENDDO
3376	ELSE
3377	DO j=1,nagec_herb
3378	vegagec_crop(:,j) = vegagec_crop(:,j)+veget_max(:,endind-j+1)
3379	ENDDO
3380	ENDIF
3381	ENDIF
3382	ENDDO
3383
3384	END SUBROUTINE calc_cover
3385
3386	! Note this subroutine does not depend on how many age classes there are
3387	! in different MTCs.
3388	SUBROUTINE glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
3389	p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
3390	IncreDeficit)
3391
3392	IMPLICIT NONE
3393
3394	!! 0.1 Input variables
3395	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
3396	INTEGER, INTENT(in) :: p2c,ipasture,g2c,igrass,f2c,itree,icrop
3397	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glcc !! the land-cover-change (LCC) matrix in case a gross LCC is
3398	!! used.
3399
3400	!! 0.2 Output variables
3401
3402
3403	!! 0.3 Modified variables
3404	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: veget_4veg !! "maximal" coverage of tree/grass/pasture/crop
3405	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccDef !! Gross LCC deficit, negative values mean that there
3406	!! are not enough fractions in the source vegetations
3407	!! to the target ones as presribed by the LCC matrix.
3408	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
3409	!! after considering the consistency between presribed
3410	!! glcc matrix and existing vegetation fractions.
3411	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
3412	!! there are not enough fractions in the source PFTs
3413	!! /vegetations to target PFTs/vegetations. I.e., these
3414	!! fraction transfers are presribed in LCC matrix but
3415	!! not realized.
3416
3417	!! 0.4 Local variables
3418	REAL(r_std), DIMENSION(npts) :: tmpdef !! LCC deficits by summing up all the deficits to the
3419	!! the same target vegetation.
3420
3421
3422	!! 0. We first handle the cases where veget_4veg might be very small
3423	!tree
3424	WHERE(veget_4veg(:,itree) > min_stomate)
3425	glccDef(:,f2c) = veget_4veg(:,itree)-glcc(:,f2c)
3426	WHERE(veget_4veg(:,itree)>glcc(:,f2c))
3427	glccReal(:,f2c) = glcc(:,f2c)
3428	ELSEWHERE
3429	glccReal(:,f2c) = veget_4veg(:,itree)
3430	ENDWHERE
3431	ELSEWHERE
3432	glccReal(:,f2c) = 0.
3433	glccDef(:,f2c) = -1*glcc(:,f2c)
3434	ENDWHERE
3435
3436	!pasture
3437	WHERE(veget_4veg(:,ipasture) > min_stomate)
3438	glccDef(:,p2c) = veget_4veg(:,ipasture)-glcc(:,p2c)
3439	WHERE(veget_4veg(:,ipasture)>glcc(:,p2c))
3440	glccReal(:,p2c) = glcc(:,p2c)
3441	ELSEWHERE
3442	glccReal(:,p2c) = veget_4veg(:,ipasture)
3443	ENDWHERE
3444	ELSEWHERE
3445	glccReal(:,p2c) = 0.
3446	glccDef(:,p2c) = -1*glcc(:,p2c)
3447	ENDWHERE
3448
3449	!grass
3450	WHERE(veget_4veg(:,igrass) > min_stomate)
3451	glccDef(:,g2c) = veget_4veg(:,igrass)-glcc(:,g2c)
3452	WHERE(veget_4veg(:,igrass)>glcc(:,g2c))
3453	glccReal(:,g2c) = glcc(:,g2c)
3454	ELSEWHERE
3455	glccReal(:,g2c) = veget_4veg(:,igrass)
3456	ENDWHERE
3457	ELSEWHERE
3458	glccReal(:,g2c) = 0.
3459	glccDef(:,g2c) = -1*glcc(:,g2c)
3460	ENDWHERE
3461
3462	!! 1. Compensation sequence: pasture,grass,forest
3463	tmpdef(:) = glccDef(:,f2c)+glccDef(:,g2c)+glccDef(:,p2c)
3464	WHERE(glccDef(:,p2c)<0)
3465	WHERE(glccDef(:,g2c)<0)
3466	WHERE(glccDef(:,f2c)<0) ! 1 (-,-,-)
3467	IncreDeficit(:,icrop) = tmpdef(:)
3468	ELSEWHERE ! 2 (-,-,+)
3469	WHERE(tmpdef(:)>=min_stomate)
3470	glccReal(:,f2c) = glccReal(:,f2c)-glccDef(:,g2c)-glccDef(:,p2c)
3471	ELSEWHERE
3472	glccReal(:,f2c) = veget_4veg(:,itree)
3473	IncreDeficit(:,icrop) = tmpdef(:)
3474	ENDWHERE
3475	ENDWHERE
3476	ELSEWHERE
3477	WHERE(glccDef(:,f2c)<0) ! 3 (-,+,-)
3478	WHERE(tmpdef(:)>=min_stomate)
3479	glccReal(:,g2c) = glccReal(:,g2c)-glccDef(:,p2c)-glccDef(:,f2c)
3480	ELSEWHERE
3481	glccReal(:,g2c) = veget_4veg(:,igrass)
3482	IncreDeficit(:,icrop) = tmpdef(:)
3483	ENDWHERE
3484	ELSEWHERE ! 4 (-,+,+)
3485	WHERE(tmpdef(:)>=min_stomate)
3486	WHERE((glccDef(:,g2c)+glccDef(:,p2c))>=min_stomate)
3487	glccReal(:,g2c) = glccReal(:,g2c)-glccDef(:,p2c)
3488	ELSEWHERE
3489	glccReal(:,g2c) = veget_4veg(:,igrass)
3490	glccReal(:,f2c) = glccReal(:,f2c)-(glccDef(:,p2c)+glccDef(:,g2c))
3491	ENDWHERE
3492	ELSEWHERE
3493	glccReal(:,g2c) = veget_4veg(:,igrass)
3494	glccReal(:,f2c) = veget_4veg(:,itree)
3495	IncreDeficit(:,icrop) = tmpdef(:)
3496	ENDWHERE
3497	ENDWHERE
3498	ENDWHERE
3499	ELSEWHERE
3500	WHERE(glccDef(:,g2c)<0)
3501	WHERE(glccDef(:,f2c)<0) ! 5 (+,-,-)
3502	WHERE(tmpdef(:)>=min_stomate)
3503	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,g2c)-glccDef(:,f2c)
3504	ELSEWHERE
3505	IncreDeficit(:,icrop) = tmpdef(:)
3506	glccReal(:,p2c) = veget_4veg(:,ipasture)
3507	ENDWHERE
3508	ELSEWHERE ! 6 (+,-,+)
3509	WHERE(tmpdef(:)>=min_stomate)
3510	WHERE((glccDef(:,p2c)+glccDef(:,g2c))>=min_stomate)
3511	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,g2c)
3512	ELSEWHERE
3513	glccReal(:,p2c) = veget_4veg(:,ipasture)
3514	glccReal(:,f2c) = glccReal(:,f2c)-(glccDef(:,g2c)+glccDef(:,p2c))
3515	ENDWHERE
3516	ELSEWHERE
3517	IncreDeficit(:,icrop) = tmpdef(:)
3518	glccReal(:,p2c) = veget_4veg(:,ipasture)
3519	glccReal(:,f2c) = veget_4veg(:,itree)
3520	ENDWHERE
3521	ENDWHERE
3522	ELSEWHERE
3523	WHERE(glccDef(:,f2c)<0) ! 7 (+,+,-)
3524	WHERE(tmpdef(:)>=min_stomate)
3525	WHERE((glccDef(:,p2c)+glccDef(:,f2c))>=min_stomate)
3526	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,f2c)
3527	ELSEWHERE
3528	glccReal(:,p2c) = veget_4veg(:,ipasture)
3529	glccReal(:,g2c) = glccReal(:,g2c)-(glccDef(:,f2c)+glccDef(:,p2c))
3530	ENDWHERE
3531	ELSEWHERE
3532	IncreDeficit(:,icrop) = tmpdef(:)
3533	glccReal(:,g2c) = veget_4veg(:,igrass)
3534	glccReal(:,p2c) = veget_4veg(:,ipasture)
3535	ENDWHERE
3536	ELSEWHERE ! 8 (+,+,+)
3537	!do nothing
3538	ENDWHERE
3539	ENDWHERE
3540	ENDWHERE
3541	veget_4veg(:,itree) = veget_4veg(:,itree) - glccReal(:,f2c)
3542	veget_4veg(:,igrass) = veget_4veg(:,igrass) - glccReal(:,g2c)
3543	veget_4veg(:,ipasture) = veget_4veg(:,ipasture) - glccReal(:,p2c)
3544
3545	END SUBROUTINE glcc_compensation_full
3546
3547
3548
3549	!! This subroutine implements non-full compensation, is currently
3550	!! abandoned.
3551	SUBROUTINE glcc_compensation(npts,veget_4veg,glcc,glccDef, &
3552	p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
3553	IncreDeficit)
3554
3555	IMPLICIT NONE
3556
3557	!! 0.1 Input variables
3558	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
3559	REAL(r_std), DIMENSION(npts,4), INTENT(in) :: veget_4veg !! "maximal" coverage fraction of a PFT on the ground
3560	INTEGER, INTENT(in) :: p2c,ipasture,g2c,igrass,f2c,itree,icrop
3561
3562	!! 0.2 Output variables
3563
3564
3565	!! 0.3 Modified variables
3566	REAL(r_std), DIMENSION (npts,12),INTENT(inout) :: glcc !! the land-cover-change (LCC) matrix in case a gross LCC is
3567	!! used.
3568	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccDef !! Gross LCC deficit, negative values mean that there
3569	!! are not enough fractions in the source vegetations
3570	!! to the target ones as presribed by the LCC matrix.
3571	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
3572	!! there are not enough fractions in the source PFTs
3573	!! /vegetations to target PFTs/vegetations. I.e., these
3574	!! fraction transfers are presribed in LCC matrix but
3575	!! not realized.
3576
3577	!! 0.4 Local variables
3578	REAL(r_std), DIMENSION(npts) :: glccDef_all !! LCC deficits by summing up all the deficits to the
3579	!! the same target vegetation.
3580
3581
3582	WHERE(veget_4veg(:,itree) > min_stomate)
3583	glccDef(:,f2c) = veget_4veg(:,itree)-glcc(:,f2c)
3584	ELSEWHERE
3585	glccDef(:,f2c) = -1*glcc(:,f2c)
3586	glcc(:,f2c) = 0.
3587	ENDWHERE
3588
3589	WHERE(veget_4veg(:,ipasture) > min_stomate)
3590	glccDef(:,p2c) = veget_4veg(:,ipasture)-glcc(:,p2c)
3591	ELSEWHERE
3592	glccDef(:,p2c) = -1*glcc(:,p2c)
3593	glcc(:,p2c) = 0.
3594	ENDWHERE
3595
3596	WHERE(veget_4veg(:,igrass) > min_stomate)
3597	glccDef(:,g2c) = veget_4veg(:,igrass)-glcc(:,g2c)
3598	ELSEWHERE
3599	glccDef(:,g2c) = -1*glcc(:,g2c)
3600	glcc(:,g2c) = 0.
3601	ENDWHERE
3602
3603	glccDef_all(:) = glccDef(:,f2c)+glccDef(:,p2c)+glccDef(:,g2c)
3604
3605	! We allow the surpluses/deficits in p2c and g2c mutually compensating
3606	! for each other. If there are still deficits after this compensation,
3607	! they will be further compensated for by the surpluses from f2c (if there are any
3608	! surpluses). The ultimate deficits that cannot be compensated for
3609	! will be recorded and dropped.
3610
3611	! Because we assume the "pasture rule" is used, i.e., the crops
3612	! are supposed to come primarily from pastures and grasses, normally
3613	! we expect the deficits to occur in p2c or g2c rather than in f2c. But
3614	! if it happens that f2c has deficits while p2c or g2c has surpluse,
3615	! the surpluses will not be used to compensate for the f2c-deficits,
3616	! instead, we will just record and drop the f2c-deficits.
3617
3618	! In following codes for convenience we're not going to check
3619	! whether surpluses in f2c are enough to compensate for deficits
3620	! in p2c or g2c or both. Instead, we just add their deficits on top
3621	! of f2c. The issues of not-enough surpluses in f2c will be left for
3622	! the codes after this section to handle.
3623	WHERE (glccDef(:,p2c) < 0.)
3624	glcc(:,p2c) = veget_4veg(:,ipasture)
3625	WHERE (glccDef(:,g2c) < 0.)
3626	glcc(:,g2c) = veget_4veg(:,igrass)
3627	ELSEWHERE
3628	WHERE (glccDef(:,g2c)+glccDef(:,p2c) > min_stomate)
3629	glcc(:,g2c) = glcc(:,g2c)-glccDef(:,p2c)
3630	ELSEWHERE
3631	glcc(:,g2c) = veget_4veg(:,igrass)
3632	! whatever the case, we simply add the dificts to f2c
3633	glcc(:,f2c) = glcc(:,f2c)-glccDef(:,p2c)-glccDef(:,g2c)
3634	ENDWHERE
3635	ENDWHERE
3636
3637	ELSEWHERE
3638	WHERE(glccDef(:,g2c) < 0.)
3639	glcc(:,g2c) = veget_4veg(:,igrass)
3640	WHERE(glccDef(:,p2c)+glccDef(:,g2c) > min_stomate)
3641	glcc(:,p2c) = glcc(:,p2c)-glccDef(:,g2c)
3642	ELSEWHERE
3643	glcc(:,p2c) = veget_4veg(:,ipasture)
3644	! whatever the case, we simply add the dificts to f2c
3645	glcc(:,f2c) = glcc(:,f2c)-glccDef(:,p2c)-glccDef(:,g2c)
3646	ENDWHERE
3647	ELSEWHERE
3648	!Here p2c and g2c both show surplus, we're not going to check whether
3649	!glccDef(:,f2c) has negative values because we assume a "pasture rule"
3650	!is applied when constructing the gross LCC matrix, so deficits in
3651	!f2c will just be dropped but not be compensated for by the surpluses in
3652	!p2c or g2c.
3653	ENDWHERE
3654	ENDWHERE
3655
3656	! 1. We calculate again the f2c-deficit because f2c-glcc is adjusted in the
3657	! codes above as we allocated the deficits of p2c and g2c into f2c.
3658	! In cases where glccDef_all is less than zero, f2c-glcc will be larger
3659	! than available forest veget_max and we therefore limit the f2c-glcc to
3660	! available forest cover.
3661	! 2. There is (probably) a second case where glccDef_all is larger then zero,
3662	! but f2c-glcc is higher than veget_tree, i.e., Originally f2c is given a
3663	! high value that there is deficit in f2c but surpluses exist for p2c and g2c.
3664	! Normally we
3665	! assume this won't happen as explained above, given that a "pasture rule" was
3666	! used in constructing the gross LCC matrix. Nevertheless if this deos
3667	! happen, we will just drop the f2c deficit without being compensated
3668	! for by the surplus in p2c or g2c.
3669
3670	! we handle the 2nd case first
3671	WHERE(veget_4veg(:,itree) > min_stomate )
3672	WHERE(glccDef(:,f2c) < 0.)
3673	glcc(:,f2c) = veget_4veg(:,itree)
3674	WHERE (glccDef(:,p2c)+glccDef(:,g2c) > min_stomate)
3675	IncreDeficit(:,icrop) = glccDef(:,f2c)
3676	ELSEWHERE
3677	IncreDeficit(:,icrop) = glccDef_all(:)
3678	ENDWHERE
3679	ELSEWHERE
3680	WHERE(glccDef_all(:) < 0.) !handle the 1st case
3681	glcc(:,f2c) = veget_4veg(:,itree)
3682	IncreDeficit(:,icrop) = glccDef_all(:)
3683	ENDWHERE
3684	ENDWHERE
3685	ELSEWHERE
3686	WHERE(glccDef(:,p2c)+glccDef(:,g2c)>min_stomate)
3687	IncreDeficit(:,icrop) = glccDef(:,f2c)
3688	ELSEWHERE
3689	IncreDeficit(:,icrop) = glccDef_all(:)
3690	ENDWHERE
3691	ENDWHERE
3692
3693	END SUBROUTINE glcc_compensation
3694
3695
3696
3697	END MODULE stomate_glcchange_fh

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

Context Navigation

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

Download in other formats: