Context Navigation

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

Last change on this file since 6737 was 6545, checked in by chao.yue, 4 years ago
Apply strict sequence of secondary forest-agricultural turnover: 3rd youngest, then older, then 2nd youngest and youngest
File size: 123.6 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_glcchange_MulAgeC
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 implements gross land use change with age classes.
10	!!
11	!!\n DESCRIPTION: None
12	!!
13	!! RECENT CHANGE(S): Including permafrost carbon
14	!!
15	!! REFERENCE(S) : None
16	!!
17	!! SVN :
18	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/perso/albert.jornet/ORCHIDEE-MICT/src_stomate/stomate_lcchange.f90 $
19	!! $Date: 2015-07-30 15:38:45 +0200 (Thu, 30 Jul 2015) $
20	!! $Revision: 2847 $
21	!! \n
22	!_ ================================================================================================================================
23
24
25	MODULE stomate_glcchange_MulAgeC
26
27	! modules used:
28
29	USE ioipsl_para
30	USE stomate_data
31	USE pft_parameters
32	USE constantes
33	USE constantes_soil_var
34	USE stomate_gluc_common
35	USE stomate_gluc_constants
36
37	IMPLICIT NONE
38
39	PRIVATE
40	PUBLIC glcc_MulAgeC_firstday, glcc_MulAgeC, age_class_distr, type_conversion
41
42	CONTAINS
43
44	! ================================================================================================================================
45	!! SUBROUTINE : age_class_distr
46	!!
47	!>\BRIEF Redistribute biomass, litter, soilcarbon and water across
48	!! the age classes
49	!!
50	!! DESCRIPTION : Following growth, the trees from an age class may have become
51	!! too big to belong to this age class. The biomass, litter, soilcarbon and
52	!! soil water then need to be moved from one age class to the next age class.
53	!!
54	!! RECENT CHANGE(S) :
55	!!
56	!! MAIN OUTPUT VARIABLE(S) :
57	!!
58	!! REFERENCES : None
59	!!
60	!! FLOWCHART :
61	!! \n
62	!_ ================================================================================================================================
63
64	SUBROUTINE age_class_distr(npts, lalo, resolution, bound_spa, &
65	biomass, veget_max, ind, &
66	lm_lastyearmax, leaf_frac, co2_to_bm, &
67	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
68	everywhere, litter, carbon, lignin_struc, &
69	deepC_a, deepC_s, deepC_p, &
70	bm_to_litter, PFTpresent, when_growthinit,&
71	senescence, npp_longterm, gpp_daily, leaf_age, age, &
72	gdd_from_growthinit, gdd_midwinter, time_hum_min, hum_min_dormance, &
73	gdd_m5_dormance, &
74	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
75	gpp_week, resp_maint, resp_growth, npp_daily, resp_hetero)
76
77	IMPLICIT NONE
78
79	!! 0. Variable and parameter declaration
80
81	!! 0.1 Input variables
82
83	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
84	REAL(r_std),DIMENSION(npts,2),INTENT(in) :: lalo !! Geographical coordinates (latitude,longitude)
85	!! for pixels (degrees)
86	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution !! Resolution at each grid point (m)
87	!! [1=E-W, 2=N-S]
88
89	!! 0.2 Output variables
90
91
92	!! 0.3 Modified variables
93
94	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
95	!! each pixel
96	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
97	!! for deciduous trees)
98	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
99	!! unitless)
100	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
101	!! (0 to 1, unitless)
102	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
103	!! @tex $(gC m^{-2} day^{-1})$ @endtex
104	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
105	!! -5 deg C (for phenology)
106	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
107	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
108	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
109	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
110	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
111	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
112	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Net primary productivity
113	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
114	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
115	!! the growing season (days)
116	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
117	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
118	!! @tex $(m^{-2})$ @endtex
119	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
120	!! May sum to
121	!! less than unity if the pixel has
122	!! nobio area. (unitless, 0-1)
123	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
124	!! @tex ($gC m^{-2}$) @endtex
125	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
126	!! very localized (after its introduction) (?)
127	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
128	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
129	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
130	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
131	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
132	!! availability (days)
133	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: hum_min_dormance !! minimum moisture during dormance
134	!! (0-1, unitless)
135	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
136	!! (for phenology) - this is written to the
137	!! history files
138	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
139	!! for crops
140	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
141	!! C (for phenology)
142	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
143	!! leaves were lost (for phenology)
144	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
145	!! above and below ground
146	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
147	!! @tex ($gC m^{-2}$) @endtex
148	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
149	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
150	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
151	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! Age (years)
152	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
153	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
154	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
155	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
156	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
157	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
158	!! below ground @tex ($gC m^{-2}$) @endtex
159	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
160	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
161	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
162	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
163	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bound_spa !! Spatial age class boundaries.
164
165	!! 0.4 Local variables
166
167	INTEGER(i_std) :: ipts,ivm,igroup !! Indeces(unitless)
168	INTEGER(i_std) :: iele,ipar,ipft !! Indeces(unitless)
169	INTEGER(i_std) :: iagec,imbc,icirc !! Indeces(unitless)
170	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
171	INTEGER(i_std) :: ivma !! Indeces(unitless)
172	REAL(r_std) :: share_expanded !! Share of the veget_max of the existing vegetation
173	!! within a PFT over the total veget_max following
174	!! expansion of that PFT (unitless, 0-1)
175	!! @tex $(ind m^{-2})$ @endtex
176	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
177	!! mass balance chech for this routine
178	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
179	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
180	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
181	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start and end pool of this routine
182	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
183	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! Start and end pool of this routine
184	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
185	REAL(r_std), DIMENSION(nelements) :: temp_start !! Start and end pool of this routine
186	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
187	REAL(r_std), DIMENSION(nelements) :: temp_end !! Start and end pool of this routine
188	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
189	REAL(r_std), DIMENSION(nlitt,nlevs) :: litter_weight_expanded !! The fraction of litter on the expanded
190	!! PFT.
191	!! @tex $-$ @endtex
192	REAL(r_std), DIMENSION(npts,nvm) :: woodmass !! Woodmass of individuals (gC)
193	REAL(r_std), DIMENSION(npts,nvm) :: reverse_soilcarbon !!
194	REAL(r_std), DIMENSION(npts,nvm) :: agec_indicator !!
195	CHARACTER(LEN=80) :: data_filename
196
197	!_ ================================================================================================================================
198
199	IF (printlev.GE.3) WRITE(numout,*) 'Entering age class distribution'
200
201	!CALL getin_p('AgeC_Threshold_File',data_filename)
202	!CALL slowproc_read_data(npts, lalo, resolution, bound_spa, data_filename, 'matrix')
203
204	IF (.NOT. use_bound_spa) THEN
205	DO ipts = 1,npts
206	bound_spa(ipts,:) = age_class_bound(:)
207	ENDDO
208	ENDIF
209
210	!! 1. Initialize
211
212	woodmass(:,:) = biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
213	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)
214	reverse_soilcarbon(:,:) = -1 *(SUM(carbon(:,:,:),DIM=2) + &
215	SUM(SUM(litter(:,:,:,:,icarbon),DIM=2),DIM=3))
216
217	!! 1.2 Initialize check for mass balance closure
218	! The mass balance is calculated at the end of this routine
219	! in section 3. Initial biomass and wood product pool all other
220	! relevant pools were just set to zero.
221	pool_start(:,:,:) = zero
222	DO iele = 1,nelements
223
224	! co2_to_bm
225	pool_start(:,:,iele) = pool_start(:,:,iele) + co2_to_bm(:,:)
226
227	! Biomass pool + bm_to_litter
228	DO ipar = 1,nparts
229	pool_start(:,:,iele) = pool_start(:,:,iele) + &
230	(biomass(:,:,ipar,iele) + bm_to_litter(:,:,ipar,iele)) * &
231	veget_max(:,:)
232	ENDDO
233
234	! Litter pool (gC m-2) * (m2 m-2)
235	DO ilit = 1,nlitt
236	DO ilev = 1,nlevs
237	pool_start(:,:,iele) = pool_start(:,:,iele) + &
238	litter(:,ilit,:,ilev,iele) * veget_max(:,:)
239	ENDDO
240	ENDDO
241
242	! Soil carbon (gC m-2) * (m2 m-2)
243	DO icarb = 1,ncarb
244	pool_start(:,:,iele) = pool_start(:,:,iele) + &
245	carbon(:,icarb,:) * veget_max(:,:)
246	ENDDO
247
248	ENDDO
249
250
251	!! 2. Handle the merge of PFTs when one age class moves to the next one.
252
253	! Following growth, the value of age-class indicator variable
254	! from an age class may have become too big to stay
255	! in this age class. The biomass, litter, reverse_soilcarbon and soil
256	! water then need to be moved from one age class to the next age class.
257	DO ipts = 1,npts
258	! This loops over all the MTCs that we have ignoring age classes
259	DO ivma=1,nvmap
260	ivm=start_index(ivma)
261
262	! If we only have a single age class for this
263	! PFT, we can skip it.
264	IF(nagec_pft(ivma) .EQ. 1)CYCLE
265
266	IF(is_tree(ivm)) THEN
267	agec_indicator(:,:) = woodmass(:,:)
268	ELSE
269	agec_indicator(:,:) = reverse_soilcarbon(:,:)
270	ENDIF ! is_tree(ivm)
271
272	CALL check_merge_same_MTC(ipts, ivma, agec_indicator, bound_spa, &
273	biomass, veget_max, ind, &
274	lm_lastyearmax, leaf_frac, co2_to_bm, &
275	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
276	everywhere, litter, carbon, lignin_struc, &
277	deepC_a, deepC_s, deepC_p, &
278	bm_to_litter, PFTpresent, when_growthinit,&
279	senescence, npp_longterm, gpp_daily, leaf_age, age, &
280	gdd_from_growthinit, gdd_midwinter, time_hum_min, hum_min_dormance, &
281	gdd_m5_dormance, &
282	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
283	gpp_week, resp_maint, resp_growth, npp_daily, resp_hetero)
284
285	ENDDO ! Looping over MTCs
286	ENDDO ! loop over #pixels - domain size
287
288
289	!! 3. Mass balance closure
290
291	!! 3.1 Calculate components of the mass balance
292	pool_end(:,:,:) = zero
293	DO iele = 1,nelements
294
295	! co2_to_bm
296	pool_end(:,:,iele) = pool_end(:,:,iele) + co2_to_bm(:,:)
297
298	! Biomass pool + bm_to_litter
299	DO ipar = 1,nparts
300	pool_end(:,:,iele) = pool_end(:,:,iele) + &
301	(biomass(:,:,ipar,iele) + bm_to_litter(:,:,ipar,iele)) * &
302	veget_max(:,:)
303	ENDDO
304
305	! Litter pool (gC m-2) * (m2 m-2)
306	DO ilit = 1,nlitt
307	DO ilev = 1,nlevs
308	pool_end(:,:,iele) = pool_end(:,:,iele) + &
309	litter(:,ilit,:,ilev,iele) * veget_max(:,:)
310	ENDDO
311	ENDDO
312
313	! Soil carbon (gC m-2) * (m2 m-2)
314	DO icarb = 1,ncarb
315	pool_end(:,:,iele) = pool_end(:,:,iele) + &
316	carbon(:,icarb,:) * veget_max(:,:)
317	ENDDO
318	ENDDO
319
320	!! 3.2 Calculate mass balance
321	check_intern(:,:,iatm2land,icarbon) = zero
322	check_intern(:,:,iland2atm,icarbon) = -un * zero
323	check_intern(:,:,ilat2out,icarbon) = zero
324	check_intern(:,:,ilat2in,icarbon) = -un * zero
325	check_intern(:,:,ipoolchange,icarbon) = -un * (pool_end(:,:,icarbon) - pool_start(:,:,icarbon))
326	closure_intern = zero
327	DO imbc = 1,nmbcomp
328	closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + check_intern(:,:,imbc,icarbon)
329	ENDDO
330
331	!! 3.3 Write outcome of the check
332	! Sum over ivm because of age class redistribution
333	DO ipts = 1,npts
334	IF (SUM(closure_intern(ipts,:,icarbon)) .LT. min_stomate .AND. &
335	SUM(closure_intern(ipts,:,icarbon)) .GT. -min_stomate) THEN
336	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure: age_class_distr', ipts
337	ELSE
338	WRITE(numout,*) 'Error: mass balance is not closed in age_class_distr'
339	WRITE(numout,*) ' Difference, ipts, ', ipts, SUM(closure_intern(ipts,:,icarbon))
340	ENDIF
341	ENDDO
342
343	IF (printlev.GE.4) WRITE(numout,*) 'Leaving age class distribution'
344
345	END SUBROUTINE age_class_distr
346
347
348	SUBROUTINE check_merge_same_MTC(ipts, ivma, woodmass, bound_spa, &
349	biomass, veget_max, ind, &
350	lm_lastyearmax, leaf_frac, co2_to_bm, &
351	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
352	everywhere, litter, carbon, lignin_struc, &
353	deepC_a, deepC_s, deepC_p, &
354	bm_to_litter, PFTpresent, when_growthinit,&
355	senescence, npp_longterm, gpp_daily, leaf_age, age, &
356	gdd_from_growthinit, gdd_midwinter, time_hum_min,hum_min_dormance, &
357	gdd_m5_dormance, &
358	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
359	gpp_week, resp_maint, resp_growth, npp_daily, resp_hetero)
360
361	IMPLICIT NONE
362
363	!! 0. Variable and parameter declaration
364
365	!! 0.1 Input variables
366
367	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
368	INTEGER, INTENT(in) :: ivma !!
369	REAL(r_std), DIMENSION(:,:), INTENT(in) :: woodmass !! Woodmass of individuals (gC)
370	REAL(r_std), DIMENSION(:,:), INTENT(in) :: bound_spa !!
371
372	!! 0.2 Output variables
373
374
375	!! 0.3 Modified variables
376
377	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
378	!! each pixel
379	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
380	!! for deciduous trees)
381	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
382	!! unitless)
383	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
384	!! (0 to 1, unitless)
385	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
386	!! @tex $(gC m^{-2} day^{-1})$ @endtex
387	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
388	!! -5 deg C (for phenology)
389	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
390	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
391	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
392	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
393	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
394	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
395	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Net primary productivity
396	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
397	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
398	!! the growing season (days)
399	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
400	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
401	!! @tex $(m^{-2})$ @endtex
402	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
403	!! May sum to
404	!! less than unity if the pixel has
405	!! nobio area. (unitless, 0-1)
406	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
407	!! @tex ($gC m^{-2}$) @endtex
408	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
409	!! very localized (after its introduction) (?)
410	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
411	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
412	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
413	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
414	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
415	!! availability (days)
416	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: hum_min_dormance !! minimum moisture during dormance
417	!! (0-1, unitless)
418	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
419	!! (for phenology) - this is written to the
420	!! history files
421	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
422	!! for crops
423	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
424	!! C (for phenology)
425	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
426	!! leaves were lost (for phenology)
427	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
428	!! above and below ground
429	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
430	!! @tex ($gC m^{-2}$) @endtex
431	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
432	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
433	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
434	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! Age (years)
435	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
436	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
437	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
438	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
439	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
440	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
441	!! below ground @tex ($gC m^{-2}$) @endtex
442	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
443	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
444	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
445	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
446
447	!! 0.4 Local variables
448
449	INTEGER(i_std) :: iele,ipar,ipft !! Indeces(unitless)
450	INTEGER(i_std) :: iagec,imbc,icirc !! Indeces(unitless)
451	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
452	REAL(r_std) :: share_expanded !! Share of the veget_max of the existing vegetation
453	!! within a PFT over the total veget_max following
454	!! expansion of that PFT (unitless, 0-1)
455	!! @tex $(ind m^{-2})$ @endtex
456	REAL(r_std), DIMENSION(nlitt,nlevs) :: litter_weight_expanded !! The fraction of litter on the expanded
457	!! PFT.
458
459
460	!_ ================================================================================================================================
461
462	!! 1 Check if the trees still belong to this age class
463	! Note that the term age class is used but that the classes used in the
464	! code are not defined on an age criterion. Instead the biomass or
465	! or soil carbon pool is used.
466	IF (is_tree(start_index(ivma))) THEN
467	DO iagec = nagec_pft(ivma),1,-1
468
469	!start from oldest age class and then move to younger age classes.
470	ipft = start_index(ivma)+iagec-1
471
472	! Check whether woodmass exceeds boundaries of
473	! the age class. For forest PFTs, bound_spa stores the upper boundary
474	! value .
475	IF(ld_agec)THEN
476	WRITE(numout,*) 'Checking to merge for: '
477	WRITE(numout,*) 'ipft,iagec,ipts: ',ipft,iagec,ipts
478	WRITE(numout,*) 'nagec_pft,woodmass,age_class_bound: ',nagec_pft(ivma),&
479	woodmass(ipts,ipft),bound_spa(ipts,ipft)
480	ENDIF
481
482	IF ( (iagec .EQ. nagec_pft(ivma)) .AND. &
483	woodmass(ipts,ipft) .GT. bound_spa(ipts,ipft) ) THEN
484
485	! If these conditions are satisfied our woodmass is
486	! very unrealist
487	WRITE(numout,*) 'WARNING: age class indicator exceeds: ', &
488	bound_spa(ipts,ipft)
489
490	ELSEIF ( iagec .NE. nagec_pft(ivma) .AND. iagec .NE. nagec_pft(ivma)-1 .AND. &
491	woodmass(ipts,ipft) .GT. bound_spa(ipts,ipft) ) THEN
492
493	IF(ld_agec)THEN
494	WRITE(numout,*) 'Merging biomass'
495	WRITE(numout,*) 'ipts,ipft,iagec: ',ipts,ipft,iagec
496	WRITE(numout,*) 'age_class_bound: ',bound_spa(ipts,ipft)
497	WRITE(numout,*) 'woodmass: ',woodmass(ipts,ipft)
498
499	ENDIF
500
501	!! 2 Merge biomass
502	! Biomass of two age classes needs to be merged. The established
503	! vegetation is stored in ipft+1, the new vegetation is stored in
504	! ipft
505	share_expanded = veget_max(ipts,ipft+1) / &
506	( veget_max(ipts,ipft+1) + veget_max(ipts,ipft) )
507	! We also need a scaling factor which includes the litter
508	DO ilev=1,nlevs
509	DO ilit=1,nlitt
510	IF(litter(ipts,ilit,ipft,ilev,icarbon) .GE. min_stomate)THEN
511	litter_weight_expanded(ilit,ilev)=litter(ipts,ilit,ipft+1,ilev,icarbon) * veget_max(ipts,ipft+1)/ &
512	(litter(ipts,ilit,ipft+1,ilev,icarbon) * veget_max(ipts,ipft+1) + &
513	litter(ipts,ilit,ipft,ilev,icarbon) * veget_max(ipts,ipft))
514	ELSE
515	litter_weight_expanded(ilit,ilev)=zero
516	ENDIF
517	END DO
518	ENDDO
519
520
521
522	! Merge the biomass and ind of the two age classes
523	biomass(ipts,ipft+1,:,:) = share_expanded * biomass(ipts,ipft+1,:,:) + &
524	(un - share_expanded) * biomass(ipts,ipft,:,:)
525	ind(ipts,ipft+1) = share_expanded * ind(ipts,ipft+1) + &
526	(un - share_expanded) * ind(ipts,ipft)
527
528	!! 3 Empty the age class that was merged and update veget_max
529	ind(ipts,ipft) = zero
530	biomass(ipts,ipft,:,:) = zero
531	veget_max(ipts,ipft+1) = veget_max(ipts,ipft+1) + veget_max(ipts,ipft)
532	veget_max(ipts,ipft) = zero
533
534	!! 4 Calculate the PFT characteristics of the merged PFT
535	! Take the weighted mean of the existing vegetation and the new
536	! vegetation joining this PFT.
537	! Note that co2_to_bm is in gC. m-2 dt-1 ,
538	! so we should also take the weighted mean (rather than sum if
539	! this where absolute values).
540
541	lm_lastyearmax(ipts,ipft+1) = share_expanded * lm_lastyearmax(ipts,ipft+1) + &
542	(un - share_expanded) * lm_lastyearmax(ipts,ipft)
543	lm_lastyearmax(ipts,ipft) = zero
544	age(ipts,ipft+1) = share_expanded * age(ipts,ipft+1) + &
545	(un - share_expanded) * age(ipts,ipft)
546	age(ipts,ipft) = zero
547
548	!CHECK: more strictly this should be considered together with leaf mass
549	leaf_frac(ipts,ipft+1,:) = share_expanded * leaf_frac(ipts,ipft+1,:) + &
550	(un - share_expanded) * leaf_frac(ipts,ipft,:)
551	leaf_frac(ipts,ipft,:) = zero
552	leaf_age(ipts,ipft+1,:) = share_expanded * leaf_age(ipts,ipft+1,:) + &
553	(un - share_expanded) * leaf_age(ipts,ipft,:)
554	leaf_age(ipts,ipft,:) = zero
555	co2_to_bm(ipts,ipft+1) = share_expanded * co2_to_bm(ipts,ipft+1) + &
556	(un - share_expanded) * co2_to_bm(ipts,ipft)
557	co2_to_bm(ipts,ipft) = zero
558
559	! Everywhere deals with the migration of vegetation. Copy the
560	! status of the most migrated vegetation for the whole PFT
561	everywhere(ipts,ipft+1) = MAX(everywhere(ipts,ipft), everywhere(ipts,ipft+1))
562	everywhere(ipts,ipft) = zero
563
564	! The new soil&litter pools are the weighted mean of the newly
565	! established vegetation for that PFT and the soil&litter pools
566	! of the original vegetation that already exists in that PFT.
567	! Since it is not only the amount of vegetation present (veget_max) but also
568	! the amount of structural litter (litter) that is important, we have to
569	! weight by both items here.
570	DO ilev=1,nlevs
571	lignin_struc(ipts,ipft+1,ilev) = litter_weight_expanded(istructural,ilev) * lignin_struc(ipts,ipft+1,ilev) + &
572	(un - litter_weight_expanded(istructural,ilev)) * lignin_struc(ipts,ipft,ilev)
573	lignin_struc(ipts,ipft,ilev) = zero
574	ENDDO
575	litter(ipts,:,ipft+1,:,:) = share_expanded * litter(ipts,:,ipft+1,:,:) + &
576	(un - share_expanded) * litter(ipts,:,ipft,:,:)
577	litter(ipts,:,ipft,:,:) = zero
578
579	fuel_1hr(ipts,ipft+1,:,:) = share_expanded * fuel_1hr(ipts,ipft+1,:,:) + &
580	(un - share_expanded) * fuel_1hr(ipts,ipft,:,:)
581	fuel_1hr(ipts,ipft,:,:) = zero
582
583	fuel_10hr(ipts,ipft+1,:,:) = share_expanded * fuel_10hr(ipts,ipft+1,:,:) + &
584	(un - share_expanded) * fuel_10hr(ipts,ipft,:,:)
585	fuel_10hr(ipts,ipft,:,:) = zero
586
587	fuel_100hr(ipts,ipft+1,:,:) = share_expanded * fuel_100hr(ipts,ipft+1,:,:) + &
588	(un - share_expanded) * fuel_100hr(ipts,ipft,:,:)
589	fuel_100hr(ipts,ipft,:,:) = zero
590
591	fuel_1000hr(ipts,ipft+1,:,:) = share_expanded * fuel_1000hr(ipts,ipft+1,:,:) + &
592	(un - share_expanded) * fuel_1000hr(ipts,ipft,:,:)
593	fuel_1000hr(ipts,ipft,:,:) = zero
594
595	carbon(ipts,:,ipft+1) = share_expanded * carbon(ipts,:,ipft+1) + &
596	(un - share_expanded) * carbon(ipts,:,ipft)
597	carbon(ipts,:,ipft) = zero
598
599	deepC_a(ipts,:,ipft+1) = share_expanded * deepC_a(ipts,:,ipft+1) + &
600	(un - share_expanded) * deepC_a(ipts,:,ipft)
601	deepC_a(ipts,:,ipft) = zero
602
603	deepC_s(ipts,:,ipft+1) = share_expanded * deepC_s(ipts,:,ipft+1) + &
604	(un - share_expanded) * deepC_s(ipts,:,ipft)
605	deepC_s(ipts,:,ipft) = zero
606
607	deepC_p(ipts,:,ipft+1) = share_expanded * deepC_p(ipts,:,ipft+1) + &
608	(un - share_expanded) * deepC_p(ipts,:,ipft)
609	deepC_p(ipts,:,ipft) = zero
610
611	bm_to_litter(ipts,ipft+1,:,:) = share_expanded * bm_to_litter(ipts,ipft+1,:,:) + &
612	(un - share_expanded) * bm_to_litter(ipts,ipft,:,:)
613	bm_to_litter(ipts,ipft,:,:) = zero
614
615	! Copy variables that depend on veget_max
616	when_growthinit(ipts,ipft+1) = share_expanded * when_growthinit(ipts,ipft+1) + &
617	(un - share_expanded) * when_growthinit(ipts,ipft)
618	when_growthinit(ipts,ipft) = zero
619	gdd_from_growthinit(ipts,ipft+1) = share_expanded * &
620	gdd_from_growthinit(ipts,ipft+1) + &
621	(un - share_expanded) * gdd_from_growthinit(ipts,ipft)
622	gdd_from_growthinit(ipts,ipft) = zero
623	gdd_midwinter(ipts,ipft+1) = share_expanded * gdd_midwinter(ipts,ipft+1) + &
624	(un - share_expanded) * gdd_midwinter(ipts,ipft)
625	gdd_midwinter(ipts,ipft) = zero
626	time_hum_min(ipts,ipft+1) = share_expanded * time_hum_min(ipts,ipft+1) + &
627	(un - share_expanded) * time_hum_min(ipts,ipft)
628	time_hum_min(ipts,ipft) = zero
629	hum_min_dormance(ipts,ipft+1) = share_expanded * hum_min_dormance(ipts,ipft+1) + &
630	(un - share_expanded) * hum_min_dormance(ipts,ipft)
631	hum_min_dormance(ipts,ipft) = zero
632	gdd_m5_dormance(ipts,ipft+1) = share_expanded * gdd_m5_dormance(ipts,ipft+1) + &
633	(un - share_expanded) * gdd_m5_dormance(ipts,ipft)
634	gdd_m5_dormance(ipts,ipft) = zero
635	ncd_dormance(ipts,ipft+1) = share_expanded * ncd_dormance(ipts,ipft+1) + &
636	(un - share_expanded) * ncd_dormance(ipts,ipft)
637	ncd_dormance(ipts,ipft) = zero
638	moiavail_month(ipts,ipft+1) = share_expanded * moiavail_month(ipts,ipft+1) + &
639	(un - share_expanded) * moiavail_month(ipts,ipft)
640	moiavail_month(ipts,ipft) = zero
641	moiavail_week(ipts,ipft+1) = share_expanded * moiavail_week(ipts,ipft+1) + &
642	(un - share_expanded) * moiavail_week(ipts,ipft)
643	moiavail_week(ipts,ipft) = zero
644	ngd_minus5(ipts,ipft+1) = share_expanded * ngd_minus5(ipts,ipft+1) + &
645	(un - share_expanded) * ngd_minus5(ipts,ipft)
646	ngd_minus5(ipts,ipft) = zero
647
648	! Copy remaining properties
649	PFTpresent(ipts,ipft+1) = PFTpresent(ipts,ipft)
650	PFTpresent(ipts,ipft) = .FALSE.
651	senescence(ipts,ipft+1) = senescence(ipts,ipft)
652	senescence(ipts,ipft) = .FALSE.
653	npp_longterm(ipts,ipft+1) = share_expanded * npp_longterm(ipts,ipft+1) + &
654	(un - share_expanded) * npp_longterm(ipts,ipft)
655	npp_longterm(ipts,ipft) = zero
656	gpp_daily(ipts,ipft+1) = share_expanded * gpp_daily(ipts,ipft+1) + &
657	(un - share_expanded) * gpp_daily(ipts,ipft)
658	gpp_daily(ipts,ipft) = zero
659	gpp_week(ipts,ipft+1) = share_expanded * gpp_week(ipts,ipft+1) + &
660	(un - share_expanded) * gpp_week(ipts,ipft)
661	gpp_week(ipts,ipft) = zero
662	resp_maint(ipts,ipft+1) = share_expanded * resp_maint(ipts,ipft+1) + &
663	(un - share_expanded) * resp_maint(ipts,ipft)
664	resp_maint(ipts,ipft) = zero
665	resp_growth(ipts,ipft+1) = share_expanded * resp_growth(ipts,ipft+1) + &
666	(un - share_expanded) * resp_growth(ipts,ipft)
667	resp_growth(ipts,ipft) = zero
668	npp_daily(ipts,ipft+1) = share_expanded * npp_daily(ipts,ipft+1) + &
669	(un - share_expanded) * npp_daily(ipts,ipft)
670	npp_daily(ipts,ipft) = zero
671	resp_hetero(ipts,ipft+1) = share_expanded * resp_hetero(ipts,ipft+1) + &
672	(un - share_expanded) * resp_hetero(ipts,ipft)
673	resp_hetero(ipts,ipft) = zero
674
675	ENDIF
676	ENDDO
677	ENDIF
678	! ! concerned MTC is grass/pasture/crop
679	! ELSE
680	! DO iagec = 1,nagec_pft(ivma),1
681
682	! ! As the soil C gets smaller when forest-generating crop gets older,
683	! ! we start from young age class and then move to older age classes.
684	! ! If the soil C of ipft is smaller than the threshold, then it should
685	! ! go to the next age class.
686	! ipft = start_index(ivma)+iagec-1
687
688	! ! Check whether woodmass exceeds boundaries of
689	! ! the age class.
690	! IF(ld_agec)THEN
691	! WRITE(numout,*) 'Checking to merge for: '
692	! WRITE(numout,*) 'ipft,iagec,ipts: ',ipft,iagec,ipts
693	! WRITE(numout,*) 'nagec_pft,woodmass,age_class_bound: ',nagec_pft(ivma),&
694	! woodmass(ipts,ipft),bound_spa(ipts,ipft)
695	! ENDIF
696
697	! !IF ( (iagec .EQ. 1) .AND. &
698	! ! woodmass(ipts,ipft) .GT. bound_spa(ipts,ipft) ) THEN
699	! !
700	! ! ! If this is satisfied than we're having a quite large
701	! ! ! soil C in the newly initiated crop
702	! ! WRITE(numout,*) 'WARNING: age class indicator exceeds: ', &
703	! ! bound_spa(ipts,ipft)
704	!
705	! !ELSEIF ( (iagec .NE. nagec_pft(ivma)) .AND. &
706	! ! woodmass(ipts,ipft) .LT. bound_spa(ipts,ipft)) THEN
707
708	! ! If the soil C is smaller than the threshold and the concerned
709	! ! ipft is not the oldest age class, then it should move to the
710	! ! next (older) age class. So we have to set the soil C threshold
711	! ! for crop as:
712
713	! ! youngest: 0.9 of maximum end-spinup forest soil C
714	! ! 2nd young: 0.75 of maximum end-spniup forest soil C
715	! ! old: 0.55 of maximum end-spniup forest soil C
716	! ! oldest: the oldest one should not be less than zero.
717	! IF ( (iagec .NE. nagec_pft(ivma)) .AND. &
718	! woodmass(ipts,ipft) .LT. bound_spa(ipts,ipft) .AND. veget_max(ipts,ipft) .GT. min_stomate) THEN
719	! IF(ld_agec)THEN
720	! WRITE(numout,*) 'Merging biomass'
721	! WRITE(numout,*) 'ipts,ipft,iagec: ',ipts,ipft,iagec
722	! WRITE(numout,*) 'age_class_bound: ',bound_spa(ipts,ipft)
723	! WRITE(numout,*) 'woodmass: ',woodmass(ipts,ipft)
724
725	! ENDIF
726
727	! !! 2 Merge biomass
728	! ! Biomass of two age classes needs to be merged. The established
729	! ! vegetation is stored in ipft+1, the new vegetation is stored in
730	! ! ipft
731	! share_expanded = veget_max(ipts,ipft+1) / &
732	! ( veget_max(ipts,ipft+1) + veget_max(ipts,ipft) )
733	! ! We also need a scaling factor which includes the litter
734	! DO ilev=1,nlevs
735	! DO ilit=1,nlitt
736	! IF(litter(ipts,ilit,ipft,ilev,icarbon) .GE. min_stomate)THEN
737	! litter_weight_expanded(ilit,ilev)=litter(ipts,ilit,ipft+1,ilev,icarbon) * veget_max(ipts,ipft+1)/ &
738	! (litter(ipts,ilit,ipft+1,ilev,icarbon) * veget_max(ipts,ipft+1) + &
739	! litter(ipts,ilit,ipft,ilev,icarbon) * veget_max(ipts,ipft))
740	! ELSE
741	! litter_weight_expanded(ilit,ilev)=zero
742	! ENDIF
743	! END DO
744	! ENDDO
745
746	! ! Merge the biomass and ind of the two age classes
747	! biomass(ipts,ipft+1,:,:) = share_expanded * biomass(ipts,ipft+1,:,:) + &
748	! (un - share_expanded) * biomass(ipts,ipft,:,:)
749	! ind(ipts,ipft+1) = share_expanded * ind(ipts,ipft+1) + &
750	! (un - share_expanded) * ind(ipts,ipft)
751	!
752	! !! 3 Empty the age class that was merged and update veget_max
753	! ind(ipts,ipft) = zero
754	! biomass(ipts,ipft,:,:) = zero
755	! veget_max(ipts,ipft+1) = veget_max(ipts,ipft+1) + veget_max(ipts,ipft)
756	! veget_max(ipts,ipft) = zero
757
758	! !! 4 Calculate the PFT characteristics of the merged PFT
759	! ! Take the weighted mean of the existing vegetation and the new
760	! ! vegetation joining this PFT.
761	! ! Note that co2_to_bm is in gC. m-2 dt-1 ,
762	! ! so we should also take the weighted mean (rather than sum if
763	! ! this where absolute values).
764	! lm_lastyearmax(ipts,ipft+1) = share_expanded * lm_lastyearmax(ipts,ipft+1) + &
765	! (un - share_expanded) * lm_lastyearmax(ipts,ipft)
766	! lm_lastyearmax(ipts,ipft) = zero
767	! !age(ipts,ipft+1) = share_expanded * age(ipts,ipft+1) + &
768	! ! (un - share_expanded) * age(ipts,ipft)
769	! !age(ipts,ipft) = zero
770
771	! !CHECK: more strictly this should be considered together with leaf mass
772	! leaf_frac(ipts,ipft+1,:) = share_expanded * leaf_frac(ipts,ipft+1,:) + &
773	! (un - share_expanded) * leaf_frac(ipts,ipft,:)
774	! leaf_frac(ipts,ipft,:) = zero
775	! leaf_age(ipts,ipft+1,:) = share_expanded * leaf_age(ipts,ipft+1,:) + &
776	! (un - share_expanded) * leaf_age(ipts,ipft,:)
777	! leaf_age(ipts,ipft,:) = zero
778	! co2_to_bm(ipts,ipft+1) = share_expanded * co2_to_bm(ipts,ipft+1) + &
779	! (un - share_expanded) * co2_to_bm(ipts,ipft)
780	! co2_to_bm(ipts,ipft) = zero
781
782	! ! Everywhere deals with the migration of vegetation. Copy the
783	! ! status of the most migrated vegetation for the whole PFT
784	! everywhere(ipts,ipft+1) = MAX(everywhere(ipts,ipft), everywhere(ipts,ipft+1))
785	! everywhere(ipts,ipft) = zero
786
787	! ! The new soil&litter pools are the weighted mean of the newly
788	! ! established vegetation for that PFT and the soil&litter pools
789	! ! of the original vegetation that already exists in that PFT.
790	! ! Since it is not only the amount of vegetation present (veget_max) but also
791	! ! the amount of structural litter (litter) that is important, we have to
792	! ! weight by both items here.
793	! DO ilev=1,nlevs
794	! lignin_struc(ipts,ipft+1,ilev) = litter_weight_expanded(istructural,ilev) * lignin_struc(ipts,ipft+1,ilev) + &
795	! (un - litter_weight_expanded(istructural,ilev)) * lignin_struc(ipts,ipft,ilev)
796	! lignin_struc(ipts,ipft,ilev) = zero
797	! ENDDO
798	! litter(ipts,:,ipft+1,:,:) = share_expanded * litter(ipts,:,ipft+1,:,:) + &
799	! (un - share_expanded) * litter(ipts,:,ipft,:,:)
800	! litter(ipts,:,ipft,:,:) = zero
801
802	! fuel_1hr(ipts,ipft+1,:,:) = share_expanded * fuel_1hr(ipts,ipft+1,:,:) + &
803	! (un - share_expanded) * fuel_1hr(ipts,ipft,:,:)
804	! fuel_1hr(ipts,ipft,:,:) = zero
805
806	! fuel_10hr(ipts,ipft+1,:,:) = share_expanded * fuel_10hr(ipts,ipft+1,:,:) + &
807	! (un - share_expanded) * fuel_10hr(ipts,ipft,:,:)
808	! fuel_10hr(ipts,ipft,:,:) = zero
809
810	! fuel_100hr(ipts,ipft+1,:,:) = share_expanded * fuel_100hr(ipts,ipft+1,:,:) + &
811	! (un - share_expanded) * fuel_100hr(ipts,ipft,:,:)
812	! fuel_100hr(ipts,ipft,:,:) = zero
813
814	! fuel_1000hr(ipts,ipft+1,:,:) = share_expanded * fuel_1000hr(ipts,ipft+1,:,:) + &
815	! (un - share_expanded) * fuel_1000hr(ipts,ipft,:,:)
816	! fuel_1000hr(ipts,ipft,:,:) = zero
817
818	! carbon(ipts,:,ipft+1) = share_expanded * carbon(ipts,:,ipft+1) + &
819	! (un - share_expanded) * carbon(ipts,:,ipft)
820	! carbon(ipts,:,ipft) = zero
821
822	! deepC_a(ipts,:,ipft+1) = share_expanded * deepC_a(ipts,:,ipft+1) + &
823	! (un - share_expanded) * deepC_a(ipts,:,ipft)
824	! deepC_a(ipts,:,ipft) = zero
825
826	! deepC_s(ipts,:,ipft+1) = share_expanded * deepC_s(ipts,:,ipft+1) + &
827	! (un - share_expanded) * deepC_s(ipts,:,ipft)
828	! deepC_s(ipts,:,ipft) = zero
829
830	! deepC_p(ipts,:,ipft+1) = share_expanded * deepC_p(ipts,:,ipft+1) + &
831	! (un - share_expanded) * deepC_p(ipts,:,ipft)
832	! deepC_p(ipts,:,ipft) = zero
833
834	! bm_to_litter(ipts,ipft+1,:,:) = share_expanded * bm_to_litter(ipts,ipft+1,:,:) + &
835	! (un - share_expanded) * bm_to_litter(ipts,ipft,:,:)
836	! bm_to_litter(ipts,ipft,:,:) = zero
837
838	! ! Copy variables that depend on veget_max
839	! when_growthinit(ipts,ipft+1) = share_expanded * when_growthinit(ipts,ipft+1) + &
840	! (un - share_expanded) * when_growthinit(ipts,ipft)
841	! when_growthinit(ipts,ipft) = zero
842	! gdd_from_growthinit(ipts,ipft+1) = share_expanded * &
843	! gdd_from_growthinit(ipts,ipft+1) + &
844	! (un - share_expanded) * gdd_from_growthinit(ipts,ipft)
845	! gdd_from_growthinit(ipts,ipft) = zero
846	! gdd_midwinter(ipts,ipft+1) = share_expanded * gdd_midwinter(ipts,ipft+1) + &
847	! (un - share_expanded) * gdd_midwinter(ipts,ipft)
848	! gdd_midwinter(ipts,ipft) = zero
849	! time_hum_min(ipts,ipft+1) = share_expanded * time_hum_min(ipts,ipft+1) + &
850	! (un - share_expanded) * time_hum_min(ipts,ipft)
851	! time_hum_min(ipts,ipft) = zero
852	! hum_min_dormance(ipts,ipft+1) = share_expanded * hum_min_dormance(ipts,ipft+1) + &
853	! (un - share_expanded) * hum_min_dormance(ipts,ipft)
854	! hum_min_dormance(ipts,ipft) = zero
855	! gdd_m5_dormance(ipts,ipft+1) = share_expanded * gdd_m5_dormance(ipts,ipft+1) + &
856	! (un - share_expanded) * gdd_m5_dormance(ipts,ipft)
857	! gdd_m5_dormance(ipts,ipft) = zero
858	! ncd_dormance(ipts,ipft+1) = share_expanded * ncd_dormance(ipts,ipft+1) + &
859	! (un - share_expanded) * ncd_dormance(ipts,ipft)
860	! ncd_dormance(ipts,ipft) = zero
861	! moiavail_month(ipts,ipft+1) = share_expanded * moiavail_month(ipts,ipft+1) + &
862	! (un - share_expanded) * moiavail_month(ipts,ipft)
863	! moiavail_month(ipts,ipft) = zero
864	! moiavail_week(ipts,ipft+1) = share_expanded * moiavail_week(ipts,ipft+1) + &
865	! (un - share_expanded) * moiavail_week(ipts,ipft)
866	! moiavail_week(ipts,ipft) = zero
867	! ngd_minus5(ipts,ipft+1) = share_expanded * ngd_minus5(ipts,ipft+1) + &
868	! (un - share_expanded) * ngd_minus5(ipts,ipft)
869	! ngd_minus5(ipts,ipft) = zero
870	!
871	! ! Copy remaining properties
872	! PFTpresent(ipts,ipft+1) = PFTpresent(ipts,ipft)
873	! PFTpresent(ipts,ipft) = .FALSE.
874	! senescence(ipts,ipft+1) = senescence(ipts,ipft)
875	! senescence(ipts,ipft) = .FALSE.
876	! npp_longterm(ipts,ipft+1) = share_expanded * npp_longterm(ipts,ipft+1) + &
877	! (un - share_expanded) * npp_longterm(ipts,ipft)
878	! npp_longterm(ipts,ipft) = zero
879	! gpp_daily(ipts,ipft+1) = share_expanded * gpp_daily(ipts,ipft+1) + &
880	! (un - share_expanded) * gpp_daily(ipts,ipft)
881	! gpp_daily(ipts,ipft) = zero
882	! gpp_week(ipts,ipft+1) = share_expanded * gpp_week(ipts,ipft+1) + &
883	! (un - share_expanded) * gpp_week(ipts,ipft)
884	! gpp_week(ipts,ipft) = zero
885	! resp_maint(ipts,ipft+1) = share_expanded * resp_maint(ipts,ipft+1) + &
886	! (un - share_expanded) * resp_maint(ipts,ipft)
887	! resp_maint(ipts,ipft) = zero
888	! resp_growth(ipts,ipft+1) = share_expanded * resp_growth(ipts,ipft+1) + &
889	! (un - share_expanded) * resp_growth(ipts,ipft)
890	! resp_growth(ipts,ipft) = zero
891	! npp_daily(ipts,ipft+1) = share_expanded * npp_daily(ipts,ipft+1) + &
892	! (un - share_expanded) * npp_daily(ipts,ipft)
893	! npp_daily(ipts,ipft) = zero
894	! resp_hetero(ipts,ipft+1) = share_expanded * resp_hetero(ipts,ipft+1) + &
895	! (un - share_expanded) * resp_hetero(ipts,ipft)
896	! resp_hetero(ipts,ipft) = zero
897
898	! ENDIF
899	! ENDDO
900
901	! ENDIF
902
903	END SUBROUTINE check_merge_same_MTC
904
905
906
907	! ================================================================================================================================
908	!! SUBROUTINE gross_lcchange
909	!!
910	!>\BRIEF : Apply gross land cover change.
911	!!
912	!>\DESCRIPTION
913	!_ ================================================================================================================================
914	SUBROUTINE glcc_MulAgeC (npts, dt_days, newvegfrac, &
915	glccSecondShift,glccPrimaryShift,glccNetLCC,&
916	def_fuel_1hr_remain, def_fuel_10hr_remain, &
917	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
918	deforest_litter_remain, deforest_biomass_remain, &
919	convflux, &
920	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
921	veget_max, prod10, prod100, &
922	PFTpresent, senescence, moiavail_month, moiavail_week, &
923	gpp_week, ngd_minus5, resp_maint, resp_growth, &
924	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
925	ind, lm_lastyearmax, everywhere, age, &
926	co2_to_bm, gpp_daily, co2_fire, &
927	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
928	gdd_m5_dormance, ncd_dormance, &
929	lignin_struc, carbon, leaf_frac, &
930	deepC_a, deepC_s, deepC_p, &
931	leaf_age, bm_to_litter, biomass, litter, &
932	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
933
934	IMPLICIT NONE
935
936	!! 0.1 Input variables
937
938	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
939	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
940	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccSecondShift !! the land-cover-change (LCC) matrix in case a gross LCC is
941	!! used.
942	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccPrimaryShift !! the land-cover-change (LCC) matrix in case a gross LCC is
943	!! used.
944	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
945	!! used.
946	REAL(r_std), DIMENSION (npts,nvmap),INTENT(in) :: newvegfrac !! veget max fraction matrix to guide the allocation of newly
947	!! created lands of a given vegetation type.
948
949	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1hr_remain
950	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_10hr_remain
951	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_100hr_remain
952	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1000hr_remain
953	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
954	!! deforestation region.
955	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
956	!! deforestation region.
957
958
959	!! 0.2 Output variables
960	REAL(r_std), DIMENSION(npts,nwp), INTENT(inout) :: convflux !! release during first year following land cover
961	!! change
962	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
963	!! after considering the consistency between presribed
964	!! glcc matrix and existing vegetation fractions.
965	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: IncreDeficit !! Originally "Increment" deficits, negative values mean that
966	!! there are not enough fractions in the source PFTs
967	!! /vegetations to target PFTs/vegetations. I.e., these
968	!! fraction transfers are presribed in LCC matrix but
969	!! not realized. Now the glccDeficit for all land cover changes
970	!! except forestry harvest.
971	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
972	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
973	!! i.e., the contribution of each PFT to the youngest age-class of MTC
974
975	!! 0.3 Modified variables
976	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT (LAI ->
977	!! infinity) on ground (unitless)
978	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
979	!! pool after the annual release for each
980	!! compartment (10 + 1 : input from year of land
981	!! cover change)
982	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
983	!! pool after the annual release for each
984	!! compartment (100 + 1 : input from year of land
985	!! cover change)
986	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
987	!! each pixel
988	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
989	!! for deciduous trees)
990	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
991	!! unitless)
992	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
993	!! (0 to 1, unitless)
994	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
995	!! @tex $(gC m^{-2} day^{-1})$ @endtex
996	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
997	!! -5 deg C (for phenology)
998	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
999	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1000	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
1001	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1002	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
1003	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1004	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
1005	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1006	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
1007	!! the growing season (days)
1008	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
1009	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
1010	!! @tex $(m^{-2})$ @endtex
1011	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
1012	!! @tex ($gC m^{-2}$) @endtex
1013	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
1014	!! very localized (after its introduction) (?)
1015	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
1016	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
1017	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
1018	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
1019	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1020	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_fire !! Fire carbon emissions
1021	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1022	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
1023	!! availability (days)
1024	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
1025	!! (for phenology) - this is written to the
1026	!! history files
1027	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
1028	!! for crops
1029	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
1030	!! C (for phenology)
1031	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
1032	!! leaves were lost (for phenology)
1033	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1034	!! above and below ground
1035	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
1036	!! @tex ($gC m^{-2}$) @endtex
1037	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
1038	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1039	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1040	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
1041	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
1042	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
1043	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1044	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
1045	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
1046	!! below ground @tex ($gC m^{-2}$) @endtex
1047	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
1048	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
1049	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
1050	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
1051
1052	!! 0.4 Local variables
1053	REAL(r_std), DIMENSION(nvmap,nparts,nelements) :: bm_to_litter_pro !! conversion of biomass to litter
1054	!! @tex ($gC m^{-2} day^{-1}$) @endtex
1055	REAL(r_std), DIMENSION(nvmap,nparts,nelements) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
1056	REAL(r_std), DIMENSION(nvmap) :: veget_max_pro !! "maximal" coverage fraction of a PFT (LAI ->
1057	!! infinity) on ground (unitless)
1058	REAL(r_std), DIMENSION(nvmap,ncarb) :: carbon_pro !! carbon pool: active, slow, or passive
1059	!! @tex ($gC m^{-2}$) @endtex
1060	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_a_pro !! Permafrost carbon pool: active, slow, or passive
1061	!! @tex ($gC m^{-3}$) @endtex
1062	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_s_pro !! Permafrost carbon pool: active, slow, or passive
1063	!! @tex ($gC m^{-3}$) @endtex
1064	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_p_pro !! Permafrost carbon pool: active, slow, or passive
1065	!! @tex ($gC m^{-3}$) @endtex
1066	REAL(r_std), DIMENSION(nvmap,nlitt,nlevs,nelements) :: litter_pro !! metabolic and structural litter, above and
1067	!! below ground @tex ($gC m^{-2}$) @endtex
1068	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_1hr_pro
1069	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_10hr_pro
1070	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_100hr_pro
1071	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_1000hr_pro
1072	REAL(r_std), DIMENSION(nvmap,nlevs) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
1073	!! above and below ground
1074	REAL(r_std), DIMENSION(nvmap,nleafages) :: leaf_frac_pro !! fraction of leaves in leaf age class
1075	REAL(r_std), DIMENSION(nvmap,nleafages) :: leaf_age_pro !! fraction of leaves in leaf age class
1076	LOGICAL, DIMENSION(nvmap) :: PFTpresent_pro, senescence_pro !! Is pft there (unitless)
1077	REAL(r_std), DIMENSION(nvmap) :: ind_pro, age_pro, lm_lastyearmax_pro, npp_longterm_pro
1078	REAL(r_std), DIMENSION(nvmap) :: everywhere_pro
1079	REAL(r_std), DIMENSION(nvmap) :: gpp_daily_pro, npp_daily_pro, co2_to_bm_pro
1080	REAL(r_std), DIMENSION(nvmap) :: resp_maint_pro, resp_growth_pro
1081	REAL(r_std), DIMENSION(nvmap) :: resp_hetero_pro, co2_fire_pro
1082
1083	INTEGER :: ipts,ivm,ivma,l,m,ipft_young_agec
1084	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
1085
1086	REAL(r_std), DIMENSION(npts,nvmap) :: glcc_mtc !! Increase in fraction of each MTC in its youngest age-class
1087	REAL(r_std), DIMENSION(npts,nvm) :: glccReal_tmp !! A temporary variable to hold glccReal
1088
1089	WRITE(numout,*) 'Entering glcc_MulAgeC'
1090	glccReal_tmp(:,:) = zero
1091
1092	CALL glcc_MulAgeC_firstday(npts,veget_max,newvegfrac, &
1093	glccSecondShift,glccPrimaryShift,glccNetLCC,&
1094	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit)
1095
1096	glcc_mtc(:,:) = SUM(glcc_pftmtc,DIM=2)
1097
1098	DO ipts=1,npts
1099
1100	!! Initialize the _pro variables
1101	bm_to_litter_pro(:,:,:)=zero
1102	biomass_pro(:,:,:)=zero
1103	veget_max_pro(:)=zero
1104	carbon_pro(:,:)=zero
1105	deepC_a_pro(:,:)=zero
1106	deepC_s_pro(:,:)=zero
1107	deepC_p_pro(:,:)=zero
1108	litter_pro(:,:,:,:)=zero
1109	fuel_1hr_pro(:,:,:)=zero
1110	fuel_10hr_pro(:,:,:)=zero
1111	fuel_100hr_pro(:,:,:)=zero
1112	fuel_1000hr_pro(:,:,:)=zero
1113	lignin_struc_pro(:,:)=zero
1114
1115	leaf_frac_pro = zero
1116	leaf_age_pro = zero
1117	PFTpresent_pro(:) = .FALSE.
1118	senescence_pro(:) = .TRUE.
1119	ind_pro = zero
1120	age_pro = zero
1121	lm_lastyearmax_pro = zero
1122	npp_longterm_pro = zero
1123	everywhere_pro = zero
1124
1125	gpp_daily_pro=zero
1126	npp_daily_pro=zero
1127	co2_to_bm_pro=zero
1128	resp_maint_pro=zero
1129	resp_growth_pro=zero
1130	resp_hetero_pro=zero
1131	co2_fire_pro=zero
1132
1133	! Note that we assume people don't intentionally change baresoil to
1134	! vegetated land.
1135	DO ivma = 2,nvmap
1136
1137	! here we set (glcc_mtc(ipts,ivma) GT. min_stomate) as a condition,
1138	! this is necessary because later on in the subroutine of
1139	! `add_incoming_proxy_pft` we have to merge the newly established
1140	! youngest proxy with potentially exisiting youngest receiving MTC,
1141	! thus have to devide a new fraction of (frac_proxy + frac_exist),
1142	! but in case frac_exist = zero, we risk deviding by a very small value
1143	! of frac_proxy and thus we want it to be bigger than min_stomate.
1144	IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
1145
1146	! 1. we accumulate the scalar variables that will be inherited.
1147	! note that in the subroutine of `collect_legacy_pft`, all
1148	! zero transitions will be automatically skipped.
1149	CALL collect_legacy_pft(npts, ipts, ivma, glcc_pftmtc, &
1150	biomass, bm_to_litter, carbon, litter, &
1151	deepC_a, deepC_s, deepC_p, &
1152	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1153	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
1154	resp_maint, resp_growth, resp_hetero, co2_fire, &
1155	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1156	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1157	deforest_litter_remain, deforest_biomass_remain, &
1158	veget_max_pro(ivma), carbon_pro(ivma,:), &
1159	lignin_struc_pro(ivma,:), litter_pro(ivma,:,:,:), &
1160	deepC_a_pro(ivma,:), deepC_s_pro(ivma,:), deepC_p_pro(ivma,:), &
1161	fuel_1hr_pro(ivma,:,:), fuel_10hr_pro(ivma,:,:), &
1162	fuel_100hr_pro(ivma,:,:), fuel_1000hr_pro(ivma,:,:), &
1163	bm_to_litter_pro(ivma,:,:), co2_to_bm_pro(ivma), &
1164	gpp_daily_pro(ivma), npp_daily_pro(ivma), &
1165	resp_maint_pro(ivma), resp_growth_pro(ivma), &
1166	resp_hetero_pro(ivma), co2_fire_pro(ivma), &
1167	convflux,prod10,prod100)
1168
1169	! Here we substract the outgoing fraction from the source PFT.
1170	! If a too small fraction remains in this source PFT, then it is
1171	! considered exhausted and we empty it. The subroutine `empty_pft`
1172	! might be combined with `collect_legacy_pft` later.
1173	DO ivm = 1,nvm
1174	! In the above we limit the collection of legacy pools to
1175	! (glcc_mtc(ipts,ivma) GT. min_stomate). Here we tentatively use
1176	! a lower threshold of `min_stomate*0.1`.
1177	IF( glcc_pftmtc(ipts,ivm,ivma)>min_stomate*0.1 ) THEN
1178	! this is the key line to implement reduction of fraction of source
1179	! PFT.
1180	veget_max(ipts,ivm) = veget_max(ipts,ivm)-glcc_pftmtc(ipts,ivm,ivma)
1181	IF ( veget_max(ipts,ivm)<min_stomate ) THEN
1182	CALL empty_pft(ipts, ivm, veget_max, biomass, ind, &
1183	carbon, litter, lignin_struc, bm_to_litter, &
1184	deepC_a, deepC_s, deepC_p, &
1185	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1186	gpp_daily, npp_daily, gpp_week, npp_longterm, &
1187	co2_to_bm, resp_maint, resp_growth, resp_hetero, &
1188	lm_lastyearmax, leaf_frac, leaf_age, age, &
1189	everywhere, PFTpresent, when_growthinit, &
1190	senescence, gdd_from_growthinit, gdd_midwinter, &
1191	time_hum_min, gdd_m5_dormance, ncd_dormance, &
1192	moiavail_month, moiavail_week, ngd_minus5)
1193	ENDIF
1194	ENDIF
1195	ENDDO
1196
1197	ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
1198	ENDDO ! (DO ivma = 2,nvmap)
1199
1200	! We establish new youngest proxy and add it to the
1201	! existing youngest-age PFT.
1202	DO ivma = 2,nvmap
1203	ipft_young_agec = start_index(ivma)
1204
1205	! 2. we establish a proxy PFT with the fraction of veget_max_pro,
1206	! which is going to be either merged with existing target
1207	! `ipft_young_agec` PFT, or fill the place if no existing target PFT
1208	! exits.
1209	CALL initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro(ivma), &
1210	biomass_pro(ivma,:,:), co2_to_bm_pro(ivma), ind_pro(ivma), &
1211	age_pro(ivma), &
1212	senescence_pro(ivma), PFTpresent_pro(ivma), &
1213	lm_lastyearmax_pro(ivma), everywhere_pro(ivma), &
1214	npp_longterm_pro(ivma), &
1215	leaf_frac_pro(ivma,:),leaf_age_pro(ivma,:))
1216
1217	! we take as a priority from exsiting PFTs of the same meta-class
1218	! the sapling biomass needed to initialize the youngest-age-class
1219	! PFT, to avoid a too much high amount of CO2 dragged down from
1220	! the air.
1221	CALL sap_take (ipts,ivma,veget_max,biomass_pro(ivma,:,:), &
1222	biomass,co2_to_bm_pro(ivma))
1223
1224	! 3. we merge the newly initiazlized proxy PFT into existing one
1225	! or use it to fill an empty PFT slot.
1226	CALL add_incoming_proxy_pft(npts, ipts, ipft_young_agec, veget_max_pro(ivma),&
1227	carbon_pro(ivma,:), litter_pro(ivma,:,:,:), lignin_struc_pro(ivma,:), &
1228	bm_to_litter_pro(ivma,:,:), &
1229	deepC_a_pro(ivma,:), deepC_s_pro(ivma,:), deepC_p_pro(ivma,:), &
1230	fuel_1hr_pro(ivma,:,:), fuel_10hr_pro(ivma,:,:), &
1231	fuel_100hr_pro(ivma,:,:), fuel_1000hr_pro(ivma,:,:), &
1232	biomass_pro(ivma,:,:), co2_to_bm_pro(ivma), &
1233	npp_longterm_pro(ivma), ind_pro(ivma), &
1234	lm_lastyearmax_pro(ivma), age_pro(ivma), everywhere_pro(ivma), &
1235	leaf_frac_pro(ivma,:), leaf_age_pro(ivma,:), &
1236	PFTpresent_pro(ivma), senescence_pro(ivma), &
1237	gpp_daily_pro(ivma), npp_daily_pro(ivma), &
1238	resp_maint_pro(ivma), resp_growth_pro(ivma), &
1239	resp_hetero_pro(ivma), co2_fire_pro(ivma), &
1240	veget_max, carbon, litter, lignin_struc, bm_to_litter, &
1241	deepC_a, deepC_s, deepC_p, &
1242	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1243	biomass, co2_to_bm, npp_longterm, ind, &
1244	lm_lastyearmax, age, everywhere, &
1245	leaf_frac, leaf_age, PFTpresent, senescence, &
1246	gpp_daily, npp_daily, resp_maint, resp_growth, &
1247	resp_hetero, co2_fire)
1248
1249	ENDDO !(DO ivma=1,nvmap)
1250
1251	ENDDO !(DO ipts=1,npts)
1252
1253	! Write out history files
1254	CALL histwrite_p (hist_id_stomate, 'glcc_pft', itime, &
1255	glcc_pft, npts*nvm, horipft_index)
1256
1257	glccReal_tmp(:,1:12) = glccReal
1258	CALL histwrite_p (hist_id_stomate, 'glccReal', itime, &
1259	glccReal_tmp, npts*nvm, horipft_index)
1260
1261	glccReal_tmp(:,:) = zero
1262	glccReal_tmp(:,1:12) = IncreDeficit
1263	CALL histwrite_p (hist_id_stomate, 'IncreDeficit', itime, &
1264	glccReal_tmp, npts*nvm, horipft_index)
1265
1266	DO ivma = 1, nvmap
1267	WRITE(part_str,'(I2)') ivma
1268	IF (ivma < 10) part_str(1:1) = '0'
1269	CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_'//part_str(1:LEN_TRIM(part_str)), &
1270	itime, glcc_pftmtc(:,:,ivma), npts*nvm, horipft_index)
1271	ENDDO
1272
1273	END SUBROUTINE glcc_MulAgeC
1274
1275
1276	! ================================================================================================================================
1277	!! SUBROUTINE : glcc_MulAgeC_firstday
1278	!!
1279	!>\BRIEF : When necessary, adjust input glcc matrix, and allocate it
1280	!! into different contributing age classes and receiving
1281	!! youngest age classes.
1282	!! \n
1283	!_ ================================================================================================================================
1284
1285	! Note: it has this name because this subroutine will also be called
1286	! the first day of each year to precalculate the forest loss for the
1287	! deforestation fire module.
1288	SUBROUTINE glcc_MulAgeC_firstday(npts,veget_max_org,newvegfrac, &
1289	glccSecondShift,glccPrimaryShift,glccNetLCC,&
1290	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit)
1291
1292	IMPLICIT NONE
1293
1294	!! 0.1 Input variables
1295
1296	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
1297	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_org !! "maximal" coverage fraction of a PFT on the ground
1298	!! May sum to
1299	!! less than unity if the pixel has
1300	!! nobio area. (unitless, 0-1)
1301	REAL(r_std), DIMENSION(npts,nvmap), INTENT(in) :: newvegfrac !! used to guid the allocation of new PFTs.
1302	!!
1303	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccSecondShift !! the land-cover-change (LCC) matrix in case a gross LCC is
1304	!! used.
1305	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccPrimaryShift !! the land-cover-change (LCC) matrix in case a gross LCC is
1306	!! used.
1307	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
1308	!! used.
1309
1310	!! 0.2 Output variables
1311	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
1312	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
1313
1314	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! Originally the "real" glcc matrix that we apply in the model
1315	!! after considering the consistency between presribed
1316	!! glcc matrix and existing vegetation fractions. Now the glcc
1317	!! by summing SecShift,NetLCC and PriShift
1318	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: IncreDeficit !! Originally "Increment" deficits, negative values mean that
1319	!! there are not enough fractions in the source PFTs
1320	!! /vegetations to target PFTs/vegetations. I.e., these
1321	!! fraction transfers are presribed in LCC matrix but
1322	!! not realized. Now the glccDeficit for all land cover changes
1323	!! except forestry harvest.
1324
1325	!! 0.3 Modified variables
1326
1327	!! 0.4 Local variables
1328	REAL(r_std), DIMENSION(npts,nvmap) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
1329	REAL(r_std), DIMENSION(npts,nvmap) :: veget_mtc_begin !! "maximal" coverage fraction of a PFT on the ground
1330	REAL(r_std), DIMENSION(npts,nagec_tree) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
1331	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_grass !! fraction of grass age-class groups, in sequence of old->young
1332	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_pasture !! fraction of pasture age-class groups, in sequence of old->young
1333	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_crop !! fraction of crop age-class groups, in sequence of old->young
1334
1335
1336	REAL(r_std), DIMENSION(npts,4) :: veget_4veg !! "maximal" coverage fraction of a PFT on the ground
1337	REAL(r_std), DIMENSION(npts) :: veget_tree !! "maximal" coverage fraction of a PFT on the ground
1338	REAL(r_std), DIMENSION(npts) :: veget_grass !! "maximal" coverage fraction of a PFT on the ground
1339	REAL(r_std), DIMENSION(npts) :: veget_pasture !! "maximal" coverage fraction of a PFT on the ground
1340	REAL(r_std), DIMENSION(npts) :: veget_crop !! "maximal" coverage fraction of a PFT on the ground
1341
1342	REAL(r_std), DIMENSION(npts,nvm) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1343	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_tmp !! "maximal" coverage fraction of a PFT on the ground
1344	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_old !! "maximal" coverage fraction of a PFT on the ground
1345	REAL(r_std), DIMENSION(npts,nvm) :: glcc_pft_tmp !! Loss of fraction in each PFT
1346
1347	REAL(r_std), DIMENSION(npts,nvm,nvmap) :: glcc_pftmtc_SecShift !! a temporary variable to hold the fractions each PFT is going to lose
1348	REAL(r_std), DIMENSION(npts,nvm,nvmap) :: glcc_pftmtc_NetPri !! a temporary variable to hold the fractions each PFT is going to lose
1349
1350	REAL(r_std), DIMENSION(npts,12) :: glccRealSecShift !! real matrix applied for secondary shifting cultivation.
1351	REAL(r_std), DIMENSION(npts,12) :: glccRealNetPriShift !! real matrix applied for NetLCC+primary shifting cultivation.
1352
1353	REAL(r_std), DIMENSION(npts,12) :: glccDefSecShift !! deficit for the glccSecondShift
1354	REAL(r_std), DIMENSION(npts,12) :: glccDefNetPriShift !! deficit for the glccNetLCC + glccPriShift
1355
1356	REAL(r_std), DIMENSION(npts,nvm) :: glccReal_tmp !! A temporary variable to hold glccReal
1357
1358	LOGICAL, SAVE :: glcc_MulAgeC_firstday_done = .FALSE.
1359
1360	! Different indexes for convenient local uses
1361	! We define the rules for gross land cover change matrix:
1362	! 1 forest->grass
1363	! 2 forest->pasture
1364	! 3 forest->crop
1365	! 4 grass->forest
1366	! 5 grass->pasture
1367	! 6 grass->crop
1368	! 7 pasture->forest
1369	! 8 pasture->grass
1370	! 9 pasture->crop
1371	! 10 crop->forest
1372	! 11 crop->grass
1373	! 12 crop->pasture
1374	INTEGER :: f2g=1, f2p=2, f2c=3
1375	INTEGER :: g2f=4, g2p=5, g2c=6, p2f=7, p2g=8, p2c=9, c2f=10, c2g=11, c2p=12
1376
1377	INTEGER :: ivma
1378
1379	REAL(r_std), DIMENSION(npts,12) :: glccRemain
1380	REAL(r_std), DIMENSION(npts,12) :: glccSecondShift_remain
1381
1382	INTEGER :: ipts,IndStart_f,IndEnd_f
1383	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
1384
1385	!Some more local configurations
1386	LOGICAL :: allow_youngest_forest_SecShift = .TRUE.
1387
1388
1389
1390	! check for equal bi-directional transition in glccSecondShift
1391	DO ipts = 1,npts
1392	IF (ABS(glccSecondShift(ipts,f2g)-glccSecondShift(ipts,g2f)) .GE. min_stomate*100) THEN
1393	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1394	WRITE(numout,*) 'transition between f2g and g2f is not equal !'
1395	WRITE(numout,*) 'Grid cell number: ', ipts
1396	STOP
1397	ENDIF
1398
1399	IF (ABS(glccSecondShift(ipts,f2c)-glccSecondShift(ipts,c2f)) .GE. min_stomate*100) THEN
1400	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1401	WRITE(numout,*) 'transition between f2c and c2f is not equal !'
1402	WRITE(numout,*) 'Grid cell number: ', ipts
1403	STOP
1404	ENDIF
1405
1406	IF (ABS(glccSecondShift(ipts,f2p)-glccSecondShift(ipts,p2f)) .GE. min_stomate*100) THEN
1407	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1408	WRITE(numout,*) 'transition between f2p and p2f is not equal !'
1409	WRITE(numout,*) 'Grid cell number: ', ipts
1410	STOP
1411	ENDIF
1412
1413	IF (ABS(glccSecondShift(ipts,g2p)-glccSecondShift(ipts,p2g)) .GE. min_stomate*100) THEN
1414	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1415	WRITE(numout,*) 'transition between g2p and p2g is not equal !'
1416	WRITE(numout,*) 'Grid cell number: ', ipts
1417	STOP
1418	ENDIF
1419
1420	IF (ABS(glccSecondShift(ipts,g2c)-glccSecondShift(ipts,c2g)) .GE. min_stomate*100) THEN
1421	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1422	WRITE(numout,*) 'transition between g2c and c2g is not equal !'
1423	WRITE(numout,*) 'Grid cell number: ', ipts
1424	STOP
1425	ENDIF
1426
1427	IF (ABS(glccSecondShift(ipts,p2c)-glccSecondShift(ipts,c2p)) .GE. min_stomate*100) THEN
1428	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1429	WRITE(numout,*) 'transition between p2c and c2p is not equal !'
1430	WRITE(numout,*) 'Grid cell number: ', ipts
1431	STOP
1432	ENDIF
1433	ENDDO
1434
1435	! check for equal bi-directional transition in glccPrimaryShift
1436	DO ipts = 1,npts
1437	IF (ABS(glccPrimaryShift(ipts,f2g)-glccPrimaryShift(ipts,g2f)) .GE. min_stomate*100) THEN
1438	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1439	WRITE(numout,*) 'transition between f2g and g2f is not equal !'
1440	WRITE(numout,*) 'Grid cell number: ', ipts
1441	STOP
1442	ENDIF
1443
1444	IF (ABS(glccPrimaryShift(ipts,f2c)-glccPrimaryShift(ipts,c2f)) .GE. min_stomate*100) THEN
1445	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1446	WRITE(numout,*) 'transition between f2c and c2f is not equal !'
1447	WRITE(numout,*) 'Grid cell number: ', ipts
1448	STOP
1449	ENDIF
1450
1451	IF (ABS(glccPrimaryShift(ipts,f2p)-glccPrimaryShift(ipts,p2f)) .GE. min_stomate*100) THEN
1452	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1453	WRITE(numout,*) 'transition between f2p and p2f is not equal !'
1454	WRITE(numout,*) 'Grid cell number: ', ipts
1455	STOP
1456	ENDIF
1457
1458	IF (ABS(glccPrimaryShift(ipts,g2p)-glccPrimaryShift(ipts,p2g)) .GE. min_stomate*100) THEN
1459	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1460	WRITE(numout,*) 'transition between g2p and p2g is not equal !'
1461	WRITE(numout,*) 'Grid cell number: ', ipts
1462	STOP
1463	ENDIF
1464
1465	IF (ABS(glccPrimaryShift(ipts,g2c)-glccPrimaryShift(ipts,c2g)) .GE. min_stomate*100) THEN
1466	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1467	WRITE(numout,*) 'transition between g2c and c2g is not equal !'
1468	WRITE(numout,*) 'Grid cell number: ', ipts
1469	STOP
1470	ENDIF
1471
1472	IF (ABS(glccPrimaryShift(ipts,p2c)-glccPrimaryShift(ipts,c2p)) .GE. min_stomate*100) THEN
1473	WRITE(numout,*) 'Error in input transition matrix for shifting cultivation'
1474	WRITE(numout,*) 'transition between p2c and c2p is not equal !'
1475	WRITE(numout,*) 'Grid cell number: ', ipts
1476	STOP
1477	ENDIF
1478	ENDDO
1479
1480	! Initialization
1481	glccReal = zero
1482	glcc_pftmtc = zero
1483	glcc_pft = zero
1484	glcc_pft_tmp = zero
1485
1486	!!! Land cover change processes start here !!!
1487	! we make copies of original input veget_max (which is veget_max_org
1488	! in the subroutine parameter list).
1489	! veget_max will be modified through different operations in order to
1490	! check for various purposes, e.g., whether input glcc matrix
1491	! is compatible with existing veget_max and how to allocate it etc.
1492	! veget_max_old will not be modified
1493	veget_max(:,:) = veget_max_org(:,:)
1494	veget_max_old(:,:) = veget_max_org(:,:)
1495
1496	!! 3. Treat secondary-agriculture shifting cultivation transition matrix.
1497	!! The primary-agriculture shifting cultivation will be treated together
1498	!! with the netLCC transitions, with the conversion sequence of oldest->
1499	!! youngest is applied.
1500
1501	! When we prepare the driving data, secondary-agriculture shifting cultivation
1502	! is intended to include the "constant transitions" over time. Ideally, we
1503	! should start applying this secondary-agriculture shifting cultivation with
1504	! the "secondary forest" in the model. Here we tentatively start with the 3rd
1505	! youngest age class and move to the 2ne youngest age class. But if the prescribed
1506	! transition fraction is not met, we then move further to 4th youngest age class
1507	! and then move to the oldest age class sequentially.
1508
1509	! Note for the first call, we have to pass veget_mtc_begin instead of veget_mtc
1510	! in order to keep the original veget_mtc before any convserion is made. The
1511	! veget_mtc is used to in type_conversion to guide the allocation of newly
1512	! created fraction of a certain mtc to its componenet youngest PFTs.
1513	CALL calc_cover(npts,veget_max,veget_mtc_begin,vegagec_tree,vegagec_grass, &
1514	vegagec_pasture,vegagec_crop)
1515	veget_mtc = veget_mtc_begin
1516
1517	!! 3.1 We start treating secondary-agriculture cultivation from the 3rd youngest
1518	!! age class and then move to the younger age class.
1519	! Because it's rather complicated to calculate which transtion fraction between
1520	! which vegetation types should occur in case there is deficit occuring
1521	! for the overall donation vegetation type, we will just start from some
1522	! priority and leave the unrealized parts into the latter section.
1523
1524	! For this purpose, we should first make a copy of glccSecondShift into
1525	! glccRemain. glccRemain will tell us the transition fractions that have to
1526	! be treated starting from `IndStart_f+1` oldest age class and moving torward older
1527	! age class.
1528	glccRemain(:,:) = glccSecondShift(:,:)
1529
1530	! Now we will call type_conversion for each of the 12 transitions, starting
1531	! from `IndStart_f` age class moving to the 2nd youngest age class. We use glccRemain
1532	! to track the transtion fractions we should leave for the second case.
1533	! To make the code more flexible, we will store the start and end indecies
1534	! in variables.
1535
1536	!*[Note: we do above process only for forest now, as we assume the conversion
1537	! of crop/pasture/grass to other types will start always from the oldest
1538	! age class]
1539
1540	IndStart_f = nagec_tree-2 ! note the indecies and vegetfrac for tree age class
1541	! is from old to young
1542	IndEnd_f = nagec_tree-2 ! nagec_tree-2: The 3rd youngest age class
1543	! nagec_tree-1: The 2nd youngest age class
1544	! nagec_tree: The youngest age class
1545
1546	IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
1547	write(numout,*) 'glcc_MulAgeC: Age class index cannot be negative or zero!'
1548	STOP
1549	ENDIF
1550
1551	DO ipts=1,npts
1552	!f2c
1553	! Note the .TRUE. value in the last line of type_conversion call means
1554	! the forest conversion sequence is old->young.
1555	CALL type_conversion(ipts,f2c,glccSecondShift,veget_mtc,newvegfrac, &
1556	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1557	IndEnd_f,nagec_herb, &
1558	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1559	glccRemain, &
1560	.TRUE., iagec_start=IndStart_f)
1561	!f2p
1562	CALL type_conversion(ipts,f2p,glccSecondShift,veget_mtc,newvegfrac, &
1563	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
1564	IndEnd_f,nagec_herb, &
1565	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1566	glccRemain, &
1567	.TRUE., iagec_start=IndStart_f)
1568	!f2g
1569	CALL type_conversion(ipts,f2g,glccSecondShift,veget_mtc,newvegfrac, &
1570	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
1571	IndEnd_f,nagec_herb, &
1572	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1573	glccRemain, &
1574	.TRUE., iagec_start=IndStart_f)
1575	!g2c
1576	CALL type_conversion(ipts,g2c,glccSecondShift,veget_mtc,newvegfrac, &
1577	indold_grass,indagec_grass,indagec_crop,num_crop_mulagec, &
1578	1,nagec_herb, &
1579	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1580	glccRemain, &
1581	.FALSE., iagec_start=nagec_herb)
1582	!g2p
1583	CALL type_conversion(ipts,g2p,glccSecondShift,veget_mtc,newvegfrac, &
1584	indold_grass,indagec_grass,indagec_pasture,num_pasture_mulagec, &
1585	nagec_herb,nagec_herb, &
1586	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1587	glccRemain, &
1588	.TRUE.)
1589	!g2f
1590	CALL type_conversion(ipts,g2f,glccSecondShift,veget_mtc,newvegfrac, &
1591	indold_grass,indagec_grass,indagec_tree,num_tree_mulagec, &
1592	nagec_herb,nagec_tree, &
1593	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1594	glccRemain, &
1595	.TRUE.)
1596	!p2c
1597	CALL type_conversion(ipts,p2c,glccSecondShift,veget_mtc,newvegfrac, &
1598	indold_pasture,indagec_pasture,indagec_crop,num_crop_mulagec, &
1599	nagec_herb,nagec_herb, &
1600	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1601	glccRemain, &
1602	.TRUE.)
1603	!p2g
1604	CALL type_conversion(ipts,p2g,glccSecondShift,veget_mtc,newvegfrac, &
1605	indold_pasture,indagec_pasture,indagec_grass,num_grass_mulagec, &
1606	nagec_herb,nagec_herb, &
1607	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1608	glccRemain, &
1609	.TRUE.)
1610	!p2f
1611	CALL type_conversion(ipts,p2f,glccSecondShift,veget_mtc,newvegfrac, &
1612	indold_pasture,indagec_pasture,indagec_tree,num_tree_mulagec, &
1613	nagec_herb,nagec_tree, &
1614	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1615	glccRemain, &
1616	.TRUE.)
1617	!c2p
1618	CALL type_conversion(ipts,c2p,glccSecondShift,veget_mtc,newvegfrac, &
1619	indold_crop,indagec_crop,indagec_pasture,num_pasture_mulagec, &
1620	nagec_herb,nagec_herb, &
1621	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1622	glccRemain, &
1623	.TRUE.)
1624	!c2g
1625	CALL type_conversion(ipts,c2g,glccSecondShift,veget_mtc,newvegfrac, &
1626	indold_crop,indagec_crop,indagec_grass,num_grass_mulagec, &
1627	1,nagec_herb, &
1628	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1629	glccRemain, &
1630	.FALSE., iagec_start=nagec_herb)
1631	!c2f
1632	CALL type_conversion(ipts,c2f,glccSecondShift,veget_mtc,newvegfrac, &
1633	indold_crop,indagec_crop,indagec_tree,num_tree_mulagec, &
1634	1,nagec_tree, &
1635	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1636	glccRemain, &
1637	.FALSE., iagec_start=nagec_herb)
1638	ENDDO
1639	glccSecondShift_remain(:,:) = glccRemain(:,:)
1640
1641	!! 3.2 We treat the remaing unrealized transtions from forest. Now we will
1642	!! start with the 3rd oldest age class and then move to the oldest age class.
1643
1644	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
1645	vegagec_pasture,vegagec_crop)
1646	veget_mtc = veget_mtc_begin
1647
1648	IndStart_f = nagec_tree-3 ! note the indecies and vegetfrac for tree age class
1649	! is from old to young.
1650	! nagec_tree -3: The 4th youngest age class.
1651
1652	IndEnd_f = 1 ! oldest-age class forest.
1653
1654	IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
1655	write(numout,*) 'glcc_MulAgeC: Age class index cannot be negative or zero!'
1656	STOP
1657	ENDIF
1658
1659	! we start with the 3rd youngest age class and move up to the oldest age
1660	! class in the sequence of young->old, as indicated by the .FALSE. parameter
1661	! when calling the subroutine type_conversion.
1662	DO ipts=1,npts
1663	!f2c
1664	CALL type_conversion(ipts,f2c,glccSecondShift_remain,veget_mtc,newvegfrac, &
1665	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1666	IndEnd_f,nagec_herb, &
1667	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1668	glccRemain, &
1669	.FALSE., iagec_start=IndStart_f)
1670	!f2p
1671	CALL type_conversion(ipts,f2p,glccSecondShift_remain,veget_mtc,newvegfrac, &
1672	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
1673	IndEnd_f,nagec_herb, &
1674	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1675	glccRemain, &
1676	.FALSE., iagec_start=IndStart_f)
1677	!f2g
1678	CALL type_conversion(ipts,f2g,glccSecondShift_remain,veget_mtc,newvegfrac, &
1679	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
1680	IndEnd_f,nagec_herb, &
1681	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1682	glccRemain, &
1683	.FALSE., iagec_start=IndStart_f)
1684	ENDDO
1685
1686	IF (allow_youngest_forest_SecShift) THEN
1687	!!++Temp++!!
1688	!! this block of 3.3 could be commented to remove this process as desribed
1689	!! below.
1690
1691	! [2016-04-20] This is temporarily added: Normally we assume the youngest
1692	! forest age cohort will not be cut because in a shifting cultivation, they
1693	! are grown to let the land recover from agricultural process. (Or at least)
1694	! we can set the threshold of youngest age cohort to be very small. But there
1695	! are two reasons we allow the youngest forest cohort to be cut for shifting
1696	! cultivation purpose: a). Farmers may decide to harvest the wood of a forest
1697	! and then convert to crop. We don't simulate explicitly this process because
1698	! this will depend on input land change matrix and land use data assumptions.
1699	! However,we can implicitly account for this by assuming "farmers plant young
1700	! trees after harvesting the wood, and immediately convert this young trees
1701	! to crops. b). For the sake of conserving the total sum of veget_max before
1702	! and after the transition as one, we need to allow the youngest forest cohort
1703	! eligible for cutting.
1704
1705	!! 3.3 We treat the remaing unrealized transtions from forest, allowing
1706	!! the youngest forest cohort to be cut. For this purpose, we will
1707	!! start with the 2nd youngest age class and then move to the youngest one.
1708
1709	glccSecondShift_remain(:,:) = glccRemain(:,:)
1710
1711	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
1712	vegagec_pasture,vegagec_crop)
1713	veget_mtc = veget_mtc_begin
1714
1715	! Note: the setting of index here must be consistent with those of 3.1 and 3.2
1716	IndStart_f = nagec_tree-1 ! note the indecies and vegetfrac for tree age class
1717	! is from old to young.
1718	! nagec_tree -1: The 2nd youngest age class.
1719
1720	IndEnd_f = nagec_tree ! youngest class forest.
1721
1722	IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
1723	write(numout,*) 'glcc_MulAgeC: Age class index cannot be negative or zero!'
1724	STOP
1725	ENDIF
1726
1727	! Here again we change the sequence from old to young, as indicated by the
1728	! .TRUE. value in calling the subroutine type_conversion.
1729	DO ipts=1,npts
1730	!f2c
1731	CALL type_conversion(ipts,f2c,glccSecondShift_remain,veget_mtc,newvegfrac, &
1732	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1733	IndEnd_f,nagec_herb, &
1734	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1735	glccRemain, &
1736	.TRUE., iagec_start=IndStart_f)
1737	!f2p
1738	CALL type_conversion(ipts,f2p,glccSecondShift_remain,veget_mtc,newvegfrac, &
1739	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
1740	IndEnd_f,nagec_herb, &
1741	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1742	glccRemain, &
1743	.TRUE., iagec_start=IndStart_f)
1744	!f2g
1745	CALL type_conversion(ipts,f2g,glccSecondShift_remain,veget_mtc,newvegfrac, &
1746	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
1747	IndEnd_f,nagec_herb, &
1748	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1749	glccRemain, &
1750	.TRUE., iagec_start=IndStart_f)
1751	ENDDO
1752	!! End of ++Temp++ Section 3.3
1753	ENDIF
1754
1755	! Final handling of some output variables.
1756	! we separate the glcc_pftmtc_SecShift
1757	glcc_pftmtc_SecShift = glcc_pftmtc
1758
1759	! we put the remaining glccRemain into the deficit
1760	glccDefSecShift = -1 * glccRemain
1761	glccRealSecShift = glccSecondShift - glccRemain
1762
1763	!***end block to handle secondary-agriculture shifting cultivation *****
1764
1765
1766	!! 4. Treat the transtions involving the oldest age classes, which include
1767	!! the first-time primary-agriculture cultivation and the net land cover
1768	!! transtions
1769
1770	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
1771	vegagec_pasture,vegagec_crop)
1772	veget_mtc = veget_mtc_begin
1773
1774
1775	! the variable "glccReal" is originally for storing the realized maxtrix
1776	! after considering the constraining and compensation of existing vegetation
1777	! fractions. But as this case is not allowed at the moment, we will just
1778	! simply put it as the sum of glccPrimaryShift and glccNetLCC
1779	glccReal(:,:) = glccPrimaryShift+glccNetLCC
1780
1781	! We copy the glccReal to glccRemain in order to track the remaining
1782	! prescribed transtion fraction after applying each transition by calling
1783	! the subroutine "type_conversion". For the moment this is mainly to fufill
1784	! the parameter requirement of the type_conversion subroutine.
1785	glccRemain(:,:) = glccReal(:,:)
1786
1787	! We allocate in the sequences of old->young. Within the same age-class
1788	! group, we allocate in proportion with existing PFT fractions.
1789	DO ipts=1,npts
1790	!f2c
1791	CALL type_conversion(ipts,f2c,glccReal,veget_mtc,newvegfrac, &
1792	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1793	nagec_tree,nagec_herb, &
1794	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1795	glccRemain, &
1796	.TRUE.)
1797	!f2p
1798	CALL type_conversion(ipts,f2p,glccReal,veget_mtc,newvegfrac, &
1799	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
1800	nagec_tree,nagec_herb, &
1801	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1802	glccRemain, &
1803	.TRUE.)
1804	!f2g
1805	CALL type_conversion(ipts,f2g,glccReal,veget_mtc,newvegfrac, &
1806	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
1807	nagec_tree,nagec_herb, &
1808	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1809	glccRemain, &
1810	.TRUE.)
1811	!g2c
1812	CALL type_conversion(ipts,g2c,glccReal,veget_mtc,newvegfrac, &
1813	indold_grass,indagec_grass,indagec_crop,num_crop_mulagec, &
1814	1,nagec_herb, &
1815	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1816	glccRemain, &
1817	.FALSE., iagec_start=nagec_herb)
1818	!g2p
1819	CALL type_conversion(ipts,g2p,glccReal,veget_mtc,newvegfrac, &
1820	indold_grass,indagec_grass,indagec_pasture,num_pasture_mulagec, &
1821	1,nagec_herb, &
1822	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1823	glccRemain, &
1824	.FALSE., iagec_start=nagec_herb)
1825	!g2f
1826	CALL type_conversion(ipts,g2f,glccReal,veget_mtc,newvegfrac, &
1827	indold_grass,indagec_grass,indagec_tree,num_tree_mulagec, &
1828	1,nagec_tree, &
1829	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1830	glccRemain, &
1831	.FALSE., iagec_start=nagec_herb)
1832	!p2c
1833	CALL type_conversion(ipts,p2c,glccReal,veget_mtc,newvegfrac, &
1834	indold_pasture,indagec_pasture,indagec_crop,num_crop_mulagec, &
1835	1,nagec_herb, &
1836	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1837	glccRemain, &
1838	.FALSE., iagec_start=nagec_herb)
1839	!p2g
1840	CALL type_conversion(ipts,p2g,glccReal,veget_mtc,newvegfrac, &
1841	indold_pasture,indagec_pasture,indagec_grass,num_grass_mulagec, &
1842	1,nagec_herb, &
1843	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1844	glccRemain, &
1845	.FALSE., iagec_start=nagec_herb)
1846	!p2f
1847	CALL type_conversion(ipts,p2f,glccReal,veget_mtc,newvegfrac, &
1848	indold_pasture,indagec_pasture,indagec_tree,num_tree_mulagec, &
1849	1,nagec_tree, &
1850	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1851	glccRemain, &
1852	.FALSE., iagec_start=nagec_herb)
1853	!c2p
1854	CALL type_conversion(ipts,c2p,glccReal,veget_mtc,newvegfrac, &
1855	indold_crop,indagec_crop,indagec_pasture,num_pasture_mulagec, &
1856	1,nagec_herb, &
1857	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1858	glccRemain, &
1859	.FALSE., iagec_start=nagec_herb)
1860	!c2g
1861	CALL type_conversion(ipts,c2g,glccReal,veget_mtc,newvegfrac, &
1862	indold_crop,indagec_crop,indagec_grass,num_grass_mulagec, &
1863	1,nagec_herb, &
1864	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1865	glccRemain, &
1866	.FALSE., iagec_start=nagec_herb)
1867	!c2f
1868	CALL type_conversion(ipts,c2f,glccReal,veget_mtc,newvegfrac, &
1869	indold_crop,indagec_crop,indagec_tree,num_tree_mulagec, &
1870	1,nagec_tree, &
1871	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1872	glccRemain, &
1873	.FALSE., iagec_start=nagec_herb)
1874	ENDDO
1875
1876	glccDefNetPriShift = -1 * glccRemain
1877	glccRealNetPriShift = glccPrimaryShift + glccNetLCC - glccRemain
1878	glcc_pftmtc_NetPri = glcc_pftmtc - glcc_pftmtc_SecShift
1879	glccReal = glccRealSecShift + glccRealNetPriShift
1880	! Note here IncreDeficit includes the deficit from secondary<->agriculgure shifting
1881	! cultivation and the primary<->agriculture+NetLCC transitions.
1882	IncreDeficit = glccDefSecShift + glccDefNetPriShift
1883
1884	IF (.NOT. glcc_MulAgeC_firstday_done) THEN
1885
1886	glccReal_tmp = zero
1887
1888	glccReal_tmp(:,1:12) = glccRealSecShift
1889	CALL histwrite_p (hist_id_stomate, 'glccRealSecShift', itime, &
1890	glccReal_tmp, npts*nvm, horipft_index)
1891
1892	glccReal_tmp(:,1:12) = glccRealNetPriShift
1893	CALL histwrite_p (hist_id_stomate, 'glccRealNetPriShift', itime, &
1894	glccReal_tmp, npts*nvm, horipft_index)
1895
1896	glccReal_tmp(:,1:12) = glccDefSecShift
1897	CALL histwrite_p (hist_id_stomate, 'glccDefSecShift', itime, &
1898	glccReal_tmp, npts*nvm, horipft_index)
1899
1900	glccReal_tmp(:,1:12) = glccDefNetPriShift
1901	CALL histwrite_p (hist_id_stomate, 'glccDefNetPriShift', itime, &
1902	glccReal_tmp, npts*nvm, horipft_index)
1903
1904	DO ivma = 1, nvmap
1905	WRITE(part_str,'(I2)') ivma
1906	IF (ivma < 10) part_str(1:1) = '0'
1907	CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_SF_'//part_str(1:LEN_TRIM(part_str)), &
1908	itime, glcc_pftmtc_SecShift(:,:,ivma), npts*nvm, horipft_index)
1909	ENDDO
1910
1911	DO ivma = 1, nvmap
1912	WRITE(part_str,'(I2)') ivma
1913	IF (ivma < 10) part_str(1:1) = '0'
1914	CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_NPF_'//part_str(1:LEN_TRIM(part_str)), &
1915	itime, glcc_pftmtc_NetPri(:,:,ivma), npts*nvm, horipft_index)
1916	ENDDO
1917
1918	glcc_MulAgeC_firstday_done = .TRUE.
1919	ENDIF
1920
1921	END SUBROUTINE glcc_MulAgeC_firstday
1922
1923
1924
1925	! ================================================================================================================================
1926	!! SUBROUTINE : type_conversion
1927	!>\BRIEF : Allocate outgoing into different age classes and incoming into
1928	!! yongest age-class of receiving MTCs.
1929	!!
1930	!! REMARK : The current dummy variables give an example of converting forests
1931	!! to crops.
1932	!! \n
1933	!_ ================================================================================================================================
1934	SUBROUTINE type_conversion(ipts,f2c,glccReal,veget_mtc,newvegfrac, &
1935	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1936	iagec_tree_end,nagec_receive, &
1937	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1938	glccRemain, &
1939	old_to_young, iagec_start)
1940
1941	IMPLICIT NONE
1942
1943	!! Input variables
1944	INTEGER, INTENT(in) :: ipts,f2c
1945	REAL(r_std), DIMENSION(:,:), INTENT(in) :: glccReal !! The "real" glcc matrix that we apply in the model
1946	!! after considering the consistency between presribed
1947	!! glcc matrix and existing vegetation fractions.
1948	REAL(r_std), DIMENSION(:,:), INTENT(in) :: newvegfrac
1949	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
1950	INTEGER, DIMENSION(:), INTENT(in) :: indold_tree !! Indices for PFTs giving out fractions;
1951	!! here use old tree cohort as an example. When iagec_start and
1952	!! a 'old_to_young' or 'young_to_old' sequence is prescribed,
1953	!! this index can be possibly skipped.
1954	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_tree !! Indices for other cohorts except the oldest one giving out fractions;
1955	!! here we use an example of other forest chorts except the oldest one.
1956	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_crop !! Indices for secondary basic-vegetation cohorts; The youngest age classes
1957	!! of these vegetations are going to receive fractions.
1958	!! here we use crop cohorts as an example
1959	INTEGER, INTENT(in) :: num_crop_mulagec !! number of crop MTCs with more than one age classes
1960	INTEGER, INTENT(in) :: iagec_tree_end !! End index of age classes in the giving basic types
1961	!! (i.e., tree, grass, pasture, crop)
1962	INTEGER, INTENT(in) :: nagec_receive !! number of age classes in the receiving basic types
1963	!! (i.e., tree, grass, pasture, crop), here we can use crop
1964	!! as an example, nagec=nagec_herb
1965	LOGICAL, INTENT(in) :: old_to_young !! an logical variable indicating whether we should handle donation
1966	!! vegetation in a sequence of old->young or young->old. TRUE is for
1967	!! old->young. If TRUE, the index will be in increasing sequence of
1968	!! (iagec_start,iagec_tree_end).
1969	INTEGER, OPTIONAL, INTENT(in) :: iagec_start !! starting index for iagec, this is added in order to handle
1970	!! the case of secondary forest clearing.
1971
1972	!! 1. Modified variables
1973	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
1974	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1975	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fractions each PFT is going to lose
1976	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
1977	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft_tmp !! Loss of fraction in each PFT
1978	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glccRemain !! The remaining glcc matrix after applying the conversion. I.e., it will
1979	!! record the remaining unrealized transition fraction in case the donation
1980	!! vegetation is not enough compared with prescribed transition fraction.
1981	!! This variable should be initialized the same as glccReal before it's fed
1982	!! into this function.
1983
1984	!! Local vriables
1985	INTEGER :: j,iagec,iagec_start_proxy
1986	REAL(r_std) :: frac_begin,frac_used
1987	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
1988	IF (.NOT. PRESENT(iagec_start)) THEN
1989	iagec_start_proxy=1
1990	ELSE
1991	iagec_start_proxy=iagec_start
1992	ENDIF
1993
1994	! This subroutine handles the conversion from one basic-vegetation type
1995	! to another, by calling the subroutine cross_give_receive, which handles
1996	! allocation of giving-receiving fraction among the giving age classes
1997	! and receiving basic-vegetation young age classes.
1998	! We allocate in the sequences of old->young. Within the same age-class
1999	! group, we allocate in proportion with existing PFT fractions. The same
2000	! also applies in the receiving youngest-age-class PFTs, i.e., the receiving
2001	! total fraction is allocated according to existing fractions of
2002	! MTCs of the same basic vegetation type, otherwise it will be equally
2003	! distributed.
2004
2005	frac_begin = glccReal(ipts,f2c)
2006	!DO WHILE (frac_begin>min_stomate)
2007	IF (old_to_young) THEN
2008	! note that both indagec_tree and vegagec_tree are in sequence of old->young
2009	! thus iagec_start_proxy must be smaller than iagec_tree_end
2010	DO iagec=iagec_start_proxy,iagec_tree_end,1
2011	IF (vegagec_tree(ipts,iagec)>frac_begin) THEN
2012	frac_used = frac_begin
2013	ELSE IF (vegagec_tree(ipts,iagec)>min_stomate) THEN
2014	frac_used = vegagec_tree(ipts,iagec)
2015	ELSE
2016	frac_used = 0.
2017	ENDIF
2018
2019	IF (frac_used>min_stomate) THEN
2020	IF (iagec==1) THEN
2021	! Note that vegagec_tree is fractions of tree age-class groups in the
2022	! the sequence of old->young, so iagec==1 means that we're handling
2023	! first the oldest-age-group tree PFTs.
2024	CALL cross_give_receive(ipts,frac_used,veget_mtc,newvegfrac, &
2025	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
2026	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
2027	ELSE
2028	! Note also the sequence of indagec_tree is from old->young, so by
2029	! increasing iagec, we're handling progressively the old to young
2030	! tree age-class PFTs.
2031	CALL cross_give_receive(ipts,frac_used,veget_mtc,newvegfrac, &
2032	indagec_tree(:,iagec-1),indagec_crop,nagec_receive,num_crop_mulagec, &
2033	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
2034	ENDIF
2035	frac_begin = frac_begin-frac_used
2036	vegagec_tree(ipts,iagec)=vegagec_tree(ipts,iagec)-frac_used
2037	glccRemain(ipts,f2c) = glccRemain(ipts,f2c) - frac_used
2038	ENDIF
2039	ENDDO
2040	ELSE ! in the sequence of young->old
2041	DO iagec=iagec_start_proxy,iagec_tree_end,-1
2042	IF (vegagec_tree(ipts,iagec)>frac_begin) THEN
2043	frac_used = frac_begin
2044	ELSE IF (vegagec_tree(ipts,iagec)>min_stomate) THEN
2045	frac_used = vegagec_tree(ipts,iagec)
2046	ELSE
2047	frac_used = 0.
2048	ENDIF
2049
2050	IF (frac_used>min_stomate) THEN
2051	IF (iagec==1) THEN
2052	! Note that vegagec_tree is fractions of tree age-class groups in the
2053	! the sequence of old->young, so iagec==1 means that we're handling
2054	! first the oldest-age-group tree PFTs.
2055	CALL cross_give_receive(ipts,frac_used,veget_mtc,newvegfrac, &
2056	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
2057	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
2058	ELSE
2059	! Note also the sequence of indagec_tree is from old->young, so by
2060	! increasing iagec, we're handling progressively the old to young
2061	! tree age-class PFTs.
2062	CALL cross_give_receive(ipts,frac_used,veget_mtc,newvegfrac, &
2063	indagec_tree(:,iagec-1),indagec_crop,nagec_receive,num_crop_mulagec, &
2064	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
2065	ENDIF
2066	frac_begin = frac_begin-frac_used
2067	vegagec_tree(ipts,iagec)=vegagec_tree(ipts,iagec)-frac_used
2068	glccRemain(ipts,f2c) = glccRemain(ipts,f2c) - frac_used
2069	ENDIF
2070	ENDDO
2071	ENDIF
2072	!ENDDO
2073
2074	END SUBROUTINE type_conversion
2075
2076	END MODULE stomate_glcchange_MulAgeC

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

Download in other formats: