Context Navigation

stomate_glcchange_MulAgeC.f90 @ 6940

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

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

Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/stomate_glcchange_MulAgeC.f90 @ 6940

Download in other formats: