Context Navigation

stomate_glcchange_SinAgeC.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: 71.1 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_lcchange_fh
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF This module is a copy of stomate_lcchange. It includes the forestry
10	! harvest.
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): Including permafrost carbon
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/perso/albert.jornet/ORCHIDEE-MICT/src_stomate/stomate_lcchange.f90 $
20	!! $Date: 2015-07-30 15:38:45 +0200 (Thu, 30 Jul 2015) $
21	!! $Revision: 2847 $
22	!! \n
23	!_ ================================================================================================================================
24
25
26	MODULE stomate_glcchange_SinAgeC
27
28	! modules used:
29
30	USE ioipsl_para
31	USE stomate_data
32	USE pft_parameters
33	USE constantes
34	USE constantes_soil_var
35	USE stomate_gluc_common
36	USE xios_orchidee
37
38	IMPLICIT NONE
39
40	PRIVATE
41	PUBLIC glcc_SinAgeC_firstday, glcc_SinAgeC, type_conversion
42
43	CONTAINS
44
45	! ================================================================================================================================
46	!! SUBROUTINE gross_lcchange
47	!!
48	!>\BRIEF : Apply gross land cover change.
49	!!
50	!>\DESCRIPTION
51	!_ ================================================================================================================================
52	SUBROUTINE glcc_SinAgeC (npts, dt_days, harvest_matrix,newvegfrac, &
53	glccSecondShift,glccPrimaryShift,glccNetLCC,&
54	def_fuel_1hr_remain, def_fuel_10hr_remain, &
55	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
56	deforest_litter_remain, deforest_biomass_remain, &
57	convflux, cflux_prod10, cflux_prod100, &
58	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
59	veget_max, prod10, prod100, flux10, flux100, &
60	PFTpresent, senescence, moiavail_month, moiavail_week, &
61	gpp_week, ngd_minus5, resp_maint, resp_growth, &
62	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
63	ind, lm_lastyearmax, everywhere, age, &
64	co2_to_bm, gpp_daily, co2_fire, &
65	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
66	gdd_m5_dormance, ncd_dormance, &
67	lignin_struc, carbon, leaf_frac, &
68	deepC_a, deepC_s, deepC_p, &
69	leaf_age, bm_to_litter, biomass, litter, &
70	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
71
72	IMPLICIT NONE
73
74	!! 0.1 Input variables
75
76	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
77	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
78	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccSecondShift !! the land-cover-change (LCC) matrix in case a gross LCC is
79	!! used.
80	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccPrimaryShift !! the land-cover-change (LCC) matrix in case a gross LCC is
81	!! used.
82	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
83	!! used.
84	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: harvest_matrix !!
85	!!
86
87	REAL(r_std), DIMENSION (npts,nvmap),INTENT(in) :: newvegfrac !!
88	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1hr_remain
89	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_10hr_remain
90	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_100hr_remain
91	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1000hr_remain
92	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
93	!! deforestation region.
94	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
95	!! deforestation region.
96
97
98	!! 0.2 Output variables
99	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: convflux !! release during first year following land cover
100	!! change
101	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: cflux_prod10 !! total annual release from the 10 year-turnover
102	!! pool @tex ($gC m^{-2}$) @endtex
103	REAL(r_std), DIMENSION(npts,nwp), INTENT(out) :: cflux_prod100 !! total annual release from the 100 year-
104	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
105	!! after considering the consistency between presribed
106	!! glcc matrix and existing vegetation fractions.
107	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
108	!! there are not enough fractions in the source PFTs
109	!! /vegetations to target PFTs/vegetations. I.e., these
110	!! fraction transfers are presribed in LCC matrix but
111	!! not realized.
112	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
113	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout):: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
114	!! i.e., the contribution of each PFT to the youngest age-class of MTC
115
116	!! 0.3 Modified variables
117	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT (LAI ->
118	!! infinity) on ground (unitless)
119	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
120	!! pool after the annual release for each
121	!! compartment (10 + 1 : input from year of land
122	!! cover change)
123	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
124	!! pool after the annual release for each
125	!! compartment (100 + 1 : input from year of land
126	!! cover change)
127	REAL(r_std), DIMENSION(npts,10,nwp), INTENT(inout) :: flux10 !! annual release from the 10/100 year-turnover
128	!! pool compartments
129	REAL(r_std), DIMENSION(npts,100,nwp), INTENT(inout) :: flux100 !! annual release from the 10/100 year-turnover
130	!! pool compartments
131	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
132	!! each pixel
133	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
134	!! for deciduous trees)
135	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
136	!! unitless)
137	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
138	!! (0 to 1, unitless)
139	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
140	!! @tex $(gC m^{-2} day^{-1})$ @endtex
141	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
142	!! -5 deg C (for phenology)
143	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
144	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
145	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
146	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
147	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
148	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
149	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
150	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
151	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
152	!! the growing season (days)
153	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
154	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
155	!! @tex $(m^{-2})$ @endtex
156	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
157	!! @tex ($gC m^{-2}$) @endtex
158	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
159	!! very localized (after its introduction) (?)
160	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
161	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
162	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
163	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
164	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
165	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_fire !! Fire carbon emissions
166	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
167	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
168	!! availability (days)
169	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
170	!! (for phenology) - this is written to the
171	!! history files
172	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
173	!! for crops
174	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
175	!! C (for phenology)
176	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
177	!! leaves were lost (for phenology)
178	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
179	!! above and below ground
180	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
181	!! @tex ($gC m^{-2}$) @endtex
182	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
183	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
184	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
185	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
186	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
187	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
188	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
189	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
190	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
191	!! below ground @tex ($gC m^{-2}$) @endtex
192	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
193	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
194	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
195	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
196
197	!! 0.4 Local variables
198	REAL(r_std), DIMENSION(nvmap,nparts,nelements) :: bm_to_litter_pro !! conversion of biomass to litter
199	!! @tex ($gC m^{-2} day^{-1}$) @endtex
200	REAL(r_std), DIMENSION(nvmap,nparts,nelements) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
201	REAL(r_std), DIMENSION(nvmap) :: veget_max_pro !! "maximal" coverage fraction of a PFT (LAI ->
202	!! infinity) on ground (unitless)
203	REAL(r_std), DIMENSION(nvmap,ncarb) :: carbon_pro !! carbon pool: active, slow, or passive
204	!! @tex ($gC m^{-2}$) @endtex
205	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_a_pro !! Permafrost carbon pool: active, slow, or passive
206	!! @tex ($gC m^{-3}$) @endtex
207	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_s_pro !! Permafrost carbon pool: active, slow, or passive
208	!! @tex ($gC m^{-3}$) @endtex
209	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_p_pro !! Permafrost carbon pool: active, slow, or passive
210	!! @tex ($gC m^{-3}$) @endtex
211	REAL(r_std), DIMENSION(nvmap,nlitt,nlevs,nelements) :: litter_pro !! metabolic and structural litter, above and
212	!! below ground @tex ($gC m^{-2}$) @endtex
213	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_1hr_pro
214	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_10hr_pro
215	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_100hr_pro
216	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_1000hr_pro
217	REAL(r_std), DIMENSION(nvmap,nlevs) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
218	!! above and below ground
219	REAL(r_std), DIMENSION(nvmap,nleafages) :: leaf_frac_pro !! fraction of leaves in leaf age class
220	REAL(r_std), DIMENSION(nvmap,nleafages) :: leaf_age_pro !! fraction of leaves in leaf age class
221	LOGICAL, DIMENSION(nvmap) :: PFTpresent_pro, senescence_pro !! Is pft there (unitless)
222	REAL(r_std), DIMENSION(nvmap) :: ind_pro, age_pro, lm_lastyearmax_pro, npp_longterm_pro
223	REAL(r_std), DIMENSION(nvmap) :: everywhere_pro
224	REAL(r_std), DIMENSION(nvmap) :: gpp_daily_pro, npp_daily_pro, co2_to_bm_pro
225	REAL(r_std), DIMENSION(nvmap) :: resp_maint_pro, resp_growth_pro
226	REAL(r_std), DIMENSION(nvmap) :: resp_hetero_pro, co2_fire_pro
227
228	INTEGER :: ipts,ivm,ivma,l,m,ipft_young_agec
229	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
230
231	REAL(r_std), DIMENSION(npts,nvmap) :: glcc_mtc !! Increase in fraction of each MTC in its youngest age-class
232	REAL(r_std), DIMENSION(npts,nvm) :: glccReal_tmp !! A temporary variable to hold glccReal
233	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf_final !!
234	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf_final !!
235	REAL(r_std), DIMENSION(npts) :: pf2yf_compen_sf2yf !!
236	REAL(r_std), DIMENSION(npts) :: sf2yf_compen_pf2yf !!
237	REAL(r_std), DIMENSION(npts,nvm) :: glcc_harvest !! Loss of fraction due to forestry harvest
238
239	WRITE(numout,*) 'Entering glcc_SinAgeC'
240	glcc_harvest(:,:) = zero
241	glccReal_tmp(:,:) = zero
242
243	CALL glcc_SinAgeC_firstday(npts,veget_max,newvegfrac,harvest_matrix, &
244	glccSecondShift,glccPrimaryShift,glccNetLCC,&
245	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
246	Deficit_pf2yf_final, Deficit_sf2yf_final, &
247	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
248
249	glcc_mtc(:,:) = SUM(glcc_pftmtc,DIM=2)
250	DO ipts=1,npts
251
252	!! Initialize the _pro variables
253	bm_to_litter_pro(:,:,:)=zero
254	biomass_pro(:,:,:)=zero
255	veget_max_pro(:)=zero
256	carbon_pro(:,:)=zero
257	deepC_a_pro(:,:)=zero
258	deepC_s_pro(:,:)=zero
259	deepC_p_pro(:,:)=zero
260	litter_pro(:,:,:,:)=zero
261	fuel_1hr_pro(:,:,:)=zero
262	fuel_10hr_pro(:,:,:)=zero
263	fuel_100hr_pro(:,:,:)=zero
264	fuel_1000hr_pro(:,:,:)=zero
265	lignin_struc_pro(:,:)=zero
266
267	leaf_frac_pro = zero
268	leaf_age_pro = zero
269	PFTpresent_pro(:) = .FALSE.
270	senescence_pro(:) = .TRUE.
271	ind_pro = zero
272	age_pro = zero
273	lm_lastyearmax_pro = zero
274	npp_longterm_pro = zero
275	everywhere_pro = zero
276
277	gpp_daily_pro=zero
278	npp_daily_pro=zero
279	co2_to_bm_pro=zero
280	resp_maint_pro=zero
281	resp_growth_pro=zero
282	resp_hetero_pro=zero
283	co2_fire_pro=zero
284
285	! Note that we assume people don't intentionally change baresoil to
286	! vegetated land.
287	DO ivma = 2,nvmap
288	! we assume only the youngest age class receives the incoming PFT
289	! [chaoyuejoy@gmail.com 2015-08-04] This line is commented to allow
290	! the case of only single age class being handled.
291
292	! here we set glcc_mtc(ipts,ivma) > min_stomate as a condition,
293	! this is necessary because later on in the subroutine of
294	! add_incoming_proxy_pft we have to merge the newly established
295	! youngest proxy with potentially exisiting youngest receiving MTC,
296	! thus have to devide a new fraction of (frac_proxy + frac_exist),
297	! but in case frac_exist = zero, we risk deviding by a very small value
298	! of frac_proxy and thus we want it to be bigger than min_stomate.
299	IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
300
301	! 1. we accumulate the scalar variables that will be inherited
302	! note we don't handle the case of harvesting forest because
303	! we assume glcc_pftmtc(forest->forest) would be zero and this
304	! case won't occur as it's filtered by the condition of
305	! (frac>min_stomate)
306	CALL collect_legacy_pft(npts, ipts, ivma, glcc_pftmtc, &
307	biomass, bm_to_litter, carbon, litter, &
308	deepC_a, deepC_s, deepC_p, &
309	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
310	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
311	resp_maint, resp_growth, resp_hetero, co2_fire, &
312	def_fuel_1hr_remain, def_fuel_10hr_remain, &
313	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
314	deforest_litter_remain, deforest_biomass_remain, &
315	veget_max_pro(ivma), carbon_pro(ivma,:), &
316	lignin_struc_pro(ivma,:), litter_pro(ivma,:,:,:), &
317	deepC_a_pro(ivma,:), deepC_s_pro(ivma,:), deepC_p_pro(ivma,:), &
318	fuel_1hr_pro(ivma,:,:), fuel_10hr_pro(ivma,:,:), &
319	fuel_100hr_pro(ivma,:,:), fuel_1000hr_pro(ivma,:,:), &
320	bm_to_litter_pro(ivma,:,:), co2_to_bm_pro(ivma), &
321	gpp_daily_pro(ivma), npp_daily_pro(ivma), &
322	resp_maint_pro(ivma), resp_growth_pro(ivma), &
323	resp_hetero_pro(ivma), co2_fire_pro(ivma), &
324	convflux,prod10,prod100)
325
326	!++TEMP++
327	! Here we substract the outgoing fraction from the source PFT.
328	! If a too small fraction remains in this source PFT, then it is
329	! exhausted, we empty it. The subroutine 'empty_pft' might be
330	! combined with 'collect_legacy_pft', but now we just put it here.
331	DO ivm = 1,nvm
332	IF( glcc_pftmtc(ipts,ivm,ivma)>zero ) THEN
333	veget_max(ipts,ivm) = veget_max(ipts,ivm)-glcc_pftmtc(ipts,ivm,ivma)
334	IF ( veget_max(ipts,ivm)<min_stomate ) THEN
335	CALL empty_pft(ipts, ivm, veget_max, biomass, ind, &
336	carbon, litter, lignin_struc, bm_to_litter, &
337	deepC_a, deepC_s, deepC_p, &
338	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
339	gpp_daily, npp_daily, gpp_week, npp_longterm, &
340	co2_to_bm, resp_maint, resp_growth, resp_hetero, &
341	lm_lastyearmax, leaf_frac, leaf_age, age, &
342	everywhere, PFTpresent, when_growthinit, &
343	senescence, gdd_from_growthinit, gdd_midwinter, &
344	time_hum_min, gdd_m5_dormance, ncd_dormance, &
345	moiavail_month, moiavail_week, ngd_minus5)
346	ENDIF
347	ENDIF
348	ENDDO
349
350	ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
351	ENDDO !(DO ivma = 2,nvmap)
352
353	! We can only establish new youngest proxy and add it to the
354	! existing youngest-age PFT after all the harvest is done, to
355	! avoid the dilution of harvestable biomass by the young proxy
356	! and ensure consistency. Therefore now we have to loop again
357	! over nvmap.
358	DO ivma = 2,nvmap
359	IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
360
361	ipft_young_agec = start_index(ivma)
362
363	! 2. we establish a proxy PFT with the fraction of veget_max_pro,
364	! which is going to be either merged with existing target
365	! `ipft_young_agec` PFT, or fill the place if no existing target PFT
366	! exits.
367	CALL initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro(ivma), &
368	biomass_pro(ivma,:,:), co2_to_bm_pro(ivma), ind_pro(ivma), &
369	age_pro(ivma), &
370	senescence_pro(ivma), PFTpresent_pro(ivma), &
371	lm_lastyearmax_pro(ivma), everywhere_pro(ivma), &
372	npp_longterm_pro(ivma), &
373	leaf_frac_pro(ivma,:),leaf_age_pro(ivma,:))
374
375	CALL sap_take (ipts,ivma,veget_max,biomass_pro(ivma,:,:), &
376	biomass,co2_to_bm_pro(ivma))
377
378	! 3. we merge the newly initiazlized proxy PFT into existing one
379	! or use it to fill an empty PFT slot.
380	CALL add_incoming_proxy_pft(npts, ipts, ipft_young_agec, veget_max_pro(ivma),&
381	carbon_pro(ivma,:), litter_pro(ivma,:,:,:), lignin_struc_pro(ivma,:), &
382	bm_to_litter_pro(ivma,:,:), &
383	deepC_a_pro(ivma,:), deepC_s_pro(ivma,:), deepC_p_pro(ivma,:), &
384	fuel_1hr_pro(ivma,:,:), fuel_10hr_pro(ivma,:,:), &
385	fuel_100hr_pro(ivma,:,:), fuel_1000hr_pro(ivma,:,:), &
386	biomass_pro(ivma,:,:), co2_to_bm_pro(ivma), &
387	npp_longterm_pro(ivma), ind_pro(ivma), &
388	lm_lastyearmax_pro(ivma), age_pro(ivma), everywhere_pro(ivma), &
389	leaf_frac_pro(ivma,:), leaf_age_pro(ivma,:), &
390	PFTpresent_pro(ivma), senescence_pro(ivma), &
391	gpp_daily_pro(ivma), npp_daily_pro(ivma), &
392	resp_maint_pro(ivma), resp_growth_pro(ivma), &
393	resp_hetero_pro(ivma), co2_fire_pro(ivma), &
394	veget_max, carbon, litter, lignin_struc, bm_to_litter, &
395	deepC_a, deepC_s, deepC_p, &
396	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
397	biomass, co2_to_bm, npp_longterm, ind, &
398	lm_lastyearmax, age, everywhere, &
399	leaf_frac, leaf_age, PFTpresent, senescence, &
400	gpp_daily, npp_daily, resp_maint, resp_growth, &
401	resp_hetero, co2_fire)
402
403	ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
404	ENDDO !(DO ivma=1,nvmap)
405
406	ENDDO !(DO ipts=1,npts)
407
408	!! Update 10 year-turnover pool content following flux emission
409	!! (linear decay (10%) of the initial carbon input)
410	DO l = 0, 8
411	m = 10 - l
412	cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,m,:)
413	prod10(:,m,:) = prod10(:,m-1,:) - flux10(:,m-1,:)
414	flux10(:,m,:) = flux10(:,m-1,:)
415	ENDDO
416
417	cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,1,:)
418	flux10(:,1,:) = 0.1 * prod10(:,0,:)
419	prod10(:,1,:) = prod10(:,0,:)
420
421	!! 2.4.3 update 100 year-turnover pool content following flux emission\n
422	DO l = 0, 98
423	m = 100 - l
424	cflux_prod100(:,:) = cflux_prod100(:,:) + flux100(:,m,:)
425	prod100(:,m,:) = prod100(:,m-1,:) - flux100(:,m-1,:)
426	flux100(:,m,:) = flux100(:,m-1,:)
427	ENDDO
428
429	cflux_prod100(:,:) = cflux_prod100(:,:) + flux100(:,1,:)
430	flux100(:,1,:) = 0.01 * prod100(:,0,:)
431	prod100(:,1,:) = prod100(:,0,:)
432	prod10(:,0,:) = zero
433	prod100(:,0,:) = zero
434
435	!convflux = convflux/one_year*dt_days
436	!cflux_prod10 = cflux_prod10/one_year*dt_days
437	!cflux_prod100 = cflux_prod100/one_year*dt_days
438
439	! Write out history files
440	CALL histwrite_p (hist_id_stomate, 'glcc_pft', itime, &
441	glcc_pft, npts*nvm, horipft_index)
442	CALL xios_orchidee_send_field ('glcc_pft', glcc_pft) ! kjpindex,nvm
443
444	glccReal_tmp(:,1:12) = glccReal
445	CALL histwrite_p (hist_id_stomate, 'glccReal', itime, &
446	glccReal_tmp, npts*nvm, horipft_index)
447	CALL xios_orchidee_send_field ('glccReal', glccReal_tmp) ! kjpindex,nvm
448
449	! ! Write out forestry harvest variables
450	! DO ipts = 1,npts
451	! DO ivm = 1,nvm
452	! DO ivma = 1,nvmap
453	! IF (is_tree(ivm) .AND. is_tree(start_index(ivma))) THEN
454	! glcc_harvest(ipts,ivm) = glcc_harvest(ipts,ivm) + glcc_pftmtc(ipts,ivm,ivma)
455	! ENDIF
456	! ENDDO
457	! ENDDO
458	! ENDDO
459	! CALL histwrite_p (hist_id_stomate, 'glcc_harvest', itime, &
460	! glcc_harvest, npts*nvm, horipft_index)
461
462	glccReal_tmp(:,:) = zero
463	glccReal_tmp(:,1:12) = IncreDeficit
464	CALL histwrite_p (hist_id_stomate, 'IncreDeficit', itime, &
465	glccReal_tmp, npts*nvm, horipft_index)
466	CALL xios_orchidee_send_field ('IncreDeficit', glccReal_tmp) ! kjpindex,nvm
467
468	! glccReal_tmp(:,:) = zero
469	! glccReal_tmp(:,1) = Deficit_pf2yf_final
470	! glccReal_tmp(:,2) = Deficit_sf2yf_final ! is always zero in case of
471	! ! single age class
472	! glccReal_tmp(:,3) = pf2yf_compen_sf2yf ! alawys zero for SinAgeC
473	! glccReal_tmp(:,4) = sf2yf_compen_pf2yf ! always zero for SinAgeC
474
475	! CALL histwrite_p (hist_id_stomate, 'DefiComForHarvest', itime, &
476	! glccReal_tmp, npts*nvm, horipft_index)
477
478	DO ivma = 1, nvmap
479	WRITE(part_str,'(I2)') ivma
480	IF (ivma < 10) part_str(1:1) = '0'
481	CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_'//part_str(1:LEN_TRIM(part_str)), &
482	itime, glcc_pftmtc(:,:,ivma), npts*nvm, horipft_index)
483	ENDDO
484	CALL xios_orchidee_send_field ('glcc_pftmtc', glcc_pftmtc) ! kjpindex,nvm,nvmap
485
486	END SUBROUTINE glcc_SinAgeC
487
488
489	! ================================================================================================================================
490	!! SUBROUTINE : glcc_SinAgeC_firstday
491	!!
492	!>\BRIEF : When necessary, adjust input glcc matrix, and allocate it
493	!! into different contributing age classes and receiving
494	!! youngest age classes.
495	!! \n
496	!_ ================================================================================================================================
497
498	! Note: it has this name because this subroutine will also be called
499	! the first day of each year to precalculate the forest loss for the
500	! deforestation fire module.
501	SUBROUTINE glcc_SinAgeC_firstday(npts,veget_max_org,newvegfrac,harvest_matrix,&
502	glccSecondShift,glccPrimaryShift,glccNetLCC,&
503	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
504	Deficit_pf2yf_final, Deficit_sf2yf_final, &
505	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
506
507	IMPLICIT NONE
508
509	!! 0.1 Input variables
510
511	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
512	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_org !! "maximal" coverage fraction of a PFT on the ground
513	!! May sum to
514	!! less than unity if the pixel has
515	!! nobio area. (unitless, 0-1)
516	REAL(r_std), DIMENSION(npts,12),INTENT(in) :: harvest_matrix !!
517	!!
518	REAL(r_std), DIMENSION (npts,nvmap),INTENT(in) :: newvegfrac !!
519	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccSecondShift !! the land-cover-change (LCC) matrix in case a gross LCC is
520	!! used.
521	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccPrimaryShift !! the land-cover-change (LCC) matrix in case a gross LCC is
522	!! used.
523	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
524	!! used.
525
526	!! 0.2 Output variables
527	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
528	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
529	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
530	!! after considering the consistency between presribed
531	!! glcc matrix and existing vegetation fractions.
532	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
533	!! there are not enough fractions in the source PFTs
534	!! /vegetations to target PFTs/vegetations. I.e., these
535	!! fraction transfers are presribed in LCC matrix but
536	!! not realized.
537	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Deficit_pf2yf_final !!
538	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Deficit_sf2yf_final !!
539	REAL(r_std), DIMENSION(npts), INTENT(inout) :: pf2yf_compen_sf2yf !!
540	REAL(r_std), DIMENSION(npts), INTENT(inout) :: sf2yf_compen_pf2yf !!
541
542
543	!! 0.3 Modified variables
544
545	!! 0.4 Local variables
546	REAL(r_std), DIMENSION (npts,12) :: glcc !! the land-cover-change (LCC) matrix in case a gross LCC is
547	!! used.
548	REAL(r_std), DIMENSION(npts,nvmap) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
549	REAL(r_std), DIMENSION(npts,nvmap) :: veget_mtc_begin !! "maximal" coverage fraction of a PFT on the ground
550	REAL(r_std), DIMENSION(npts,nagec_tree) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
551	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_grass !! fraction of grass age-class groups, in sequence of old->young
552	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_pasture !! fraction of pasture age-class groups, in sequence of old->young
553	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_crop !! fraction of crop age-class groups, in sequence of old->young
554
555
556	REAL(r_std), DIMENSION(npts,4) :: veget_4veg !! "maximal" coverage fraction of a PFT on the ground
557	REAL(r_std), DIMENSION(npts) :: veget_tree !! "maximal" coverage fraction of a PFT on the ground
558	REAL(r_std), DIMENSION(npts) :: veget_grass !! "maximal" coverage fraction of a PFT on the ground
559	REAL(r_std), DIMENSION(npts) :: veget_pasture !! "maximal" coverage fraction of a PFT on the ground
560	REAL(r_std), DIMENSION(npts) :: veget_crop !! "maximal" coverage fraction of a PFT on the ground
561
562	REAL(r_std), DIMENSION(npts,nvm) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
563	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_tmp !! "maximal" coverage fraction of a PFT on the ground
564	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_old !! "maximal" coverage fraction of a PFT on the ground
565	REAL(r_std), DIMENSION(npts,nvm) :: glcc_pft_tmp !! Loss of fraction in each PFT
566
567	! Different indexes for convenient local uses
568	! We define the rules for gross land cover change matrix:
569	! 1 forest->grass
570	! 2 forest->pasture
571	! 3 forest->crop
572	! 4 grass->forest
573	! 5 grass->pasture
574	! 6 grass->crop
575	! 7 pasture->forest
576	! 8 pasture->grass
577	! 9 pasture->crop
578	! 10 crop->forest
579	! 11 crop->grass
580	! 12 crop->pasture
581	INTEGER :: f2g=1, f2p=2, f2c=3
582	INTEGER :: g2f=4, g2p=5, g2c=6, p2f=7, p2g=8, p2c=9, c2f=10, c2g=11, c2p=12
583
584	INTEGER, ALLOCATABLE :: indall_tree(:) !! Indices for all tree PFTs
585	INTEGER, ALLOCATABLE :: indold_tree(:) !! Indices for old tree cohort only
586	INTEGER, ALLOCATABLE :: indagec_tree(:,:) !! Indices for secondary tree cohorts,
587	!! note the sequence is old->young.
588	INTEGER, ALLOCATABLE :: indall_grass(:) !! Indices for all grass PFTs
589	INTEGER, ALLOCATABLE :: indold_grass(:) !! Indices for old grasses only
590	INTEGER, ALLOCATABLE :: indagec_grass(:,:) !! Indices for secondary grass cohorts
591	!! note the sequence is old->young.
592	INTEGER, ALLOCATABLE :: indall_pasture(:) !! Indices for all pasture PFTs
593	INTEGER, ALLOCATABLE :: indold_pasture(:) !! Indices for old pasture only
594	INTEGER, ALLOCATABLE :: indagec_pasture(:,:) !! Indices for secondary pasture cohorts
595	!! note the sequence is old->young.
596	INTEGER, ALLOCATABLE :: indall_crop(:) !! Indices for all crop PFTs
597	INTEGER, ALLOCATABLE :: indold_crop(:) !! Indices for old crops only
598	INTEGER, ALLOCATABLE :: indagec_crop(:,:) !! Indices for secondary crop cohorts
599	!! note the sequence is old->young.
600	INTEGER :: num_tree_sinagec,num_tree_mulagec,num_grass_sinagec,num_grass_mulagec, &
601	num_pasture_sinagec,num_pasture_mulagec,num_crop_sinagec,num_crop_mulagec, &
602	itree,itree2,igrass,igrass2,ipasture,ipasture2,icrop,icrop2,pf2yf,sf2yf
603	INTEGER :: i,j,ivma,staind,endind,ivm
604
605
606	REAL(r_std), DIMENSION(npts,12) :: glccDef !! Gross LCC deficit, negative values mean that there
607	!! are not enough fractions in the source vegetations
608	!! to the target ones as presribed by the LCC matrix.
609	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf !!
610	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf !!
611	REAL(r_std), DIMENSION(npts) :: Surplus_pf2yf !!
612	REAL(r_std), DIMENSION(npts) :: Surplus_sf2yf !!
613	REAL(r_std), DIMENSION(npts,12) :: glccRemain !!
614	REAL(r_std), DIMENSION(npts,12) :: HmatrixReal !!
615	INTEGER :: ipts
616
617
618	!! 1. We first build all different indices that we are going to use
619	!! in handling the PFT exchanges, three types of indices are built:
620	!! - for all age classes
621	!! - include only oldest age classes
622	!! - include all age classes excpet the oldest ones
623	! We have to build these indices because we would like to extract from
624	! donating PFTs in the sequnce of old->young age classes, and add in the
625	! receving PFTs only in the youngest-age-class PFTs. These indicies allow
626	! us to know where the different age classes are.
627
628	num_tree_sinagec=0 ! number of tree PFTs with only one single age class
629	! considered as the oldest age class
630	num_tree_mulagec=0 ! number of tree PFTs having multiple age classes
631	num_grass_sinagec=0
632	num_grass_mulagec=0
633	num_pasture_sinagec=0
634	num_pasture_mulagec=0
635	num_crop_sinagec=0
636	num_crop_mulagec=0
637
638	!! 1.1 Calculate the number of PFTs for different MTCs and allocate
639	!! the old and all indices arrays.
640
641	! [Note here the sequence to identify tree,pasture,grass,crop] is
642	! critical. The similar sequence is used in the subroutine "calc_cover".
643	! Do not forget to change the sequence there if you modify here.
644	DO ivma =2,nvmap
645	staind=start_index(ivma)
646	IF (nagec_pft(ivma)==1) THEN
647	IF (is_tree(staind)) THEN
648	num_tree_sinagec = num_tree_sinagec+1
649	ELSE IF (is_grassland_manag(staind)) THEN
650	num_pasture_sinagec = num_pasture_sinagec+1
651	ELSE IF (natural(staind)) THEN
652	num_grass_sinagec = num_grass_sinagec+1
653	ELSE
654	num_crop_sinagec = num_crop_sinagec+1
655	ENDIF
656
657	ELSE
658	IF (is_tree(staind)) THEN
659	num_tree_mulagec = num_tree_mulagec+1
660	ELSE IF (is_grassland_manag(staind)) THEN
661	num_pasture_mulagec = num_pasture_mulagec+1
662	ELSE IF (natural(staind)) THEN
663	num_grass_mulagec = num_grass_mulagec+1
664	ELSE
665	num_crop_mulagec = num_crop_mulagec+1
666	ENDIF
667	ENDIF
668	ENDDO
669
670	!! Allocate index array
671	! allocate all index
672	ALLOCATE(indall_tree(num_tree_sinagec+num_tree_mulagec*nagec_tree))
673	ALLOCATE(indall_grass(num_grass_sinagec+num_grass_mulagec*nagec_herb))
674	ALLOCATE(indall_pasture(num_pasture_sinagec+num_pasture_mulagec*nagec_herb))
675	ALLOCATE(indall_crop(num_crop_sinagec+num_crop_mulagec*nagec_herb))
676
677	! allocate old-ageclass index
678	ALLOCATE(indold_tree(num_tree_sinagec+num_tree_mulagec))
679	ALLOCATE(indold_grass(num_grass_sinagec+num_grass_mulagec))
680	ALLOCATE(indold_pasture(num_pasture_sinagec+num_pasture_mulagec))
681	ALLOCATE(indold_crop(num_crop_sinagec+num_crop_mulagec))
682
683	!! 1.2 Fill the oldest-age-class and all index arrays
684	itree=0
685	igrass=0
686	ipasture=0
687	icrop=0
688	itree2=1
689	igrass2=1
690	ipasture2=1
691	icrop2=1
692	DO ivma =2,nvmap
693	staind=start_index(ivma)
694	IF (is_tree(staind)) THEN
695	itree=itree+1
696	indold_tree(itree) = staind+nagec_pft(ivma)-1
697	DO j = 0,nagec_pft(ivma)-1
698	indall_tree(itree2+j) = staind+j
699	ENDDO
700	itree2=itree2+nagec_pft(ivma)
701	ELSE IF (natural(staind) .AND. .NOT. is_grassland_manag(staind)) THEN
702	igrass=igrass+1
703	indold_grass(igrass) = staind+nagec_pft(ivma)-1
704	DO j = 0,nagec_pft(ivma)-1
705	indall_grass(igrass2+j) = staind+j
706	ENDDO
707	igrass2=igrass2+nagec_pft(ivma)
708	ELSE IF (is_grassland_manag(staind)) THEN
709	ipasture = ipasture+1
710	indold_pasture(ipasture) = staind+nagec_pft(ivma)-1
711	DO j = 0,nagec_pft(ivma)-1
712	indall_pasture(ipasture2+j) = staind+j
713	ENDDO
714	ipasture2=ipasture2+nagec_pft(ivma)
715	ELSE
716	icrop = icrop+1
717	indold_crop(icrop) = staind+nagec_pft(ivma)-1
718	DO j = 0,nagec_pft(ivma)-1
719	indall_crop(icrop2+j) = staind+j
720	ENDDO
721	icrop2=icrop2+nagec_pft(ivma)
722	ENDIF
723	ENDDO
724
725	!! 1.3 Allocate and fill other age class index
726
727	! [chaoyuejoy@gmail.com 2015-08-05]
728	! note that we treat the case of (num_tree_mulagec==0) differently. In this
729	! case there is no distinction of age groups among tree PFTs. But we still
730	! we want to use the "gross_lcchange" subroutine. In this case we consider
731	! them as having a single age group. In the subroutines
732	! of "type_conversion" and "cross_give_receive", only the youngest-age-group
733	! PFTs of a given MTC or vegetation type could receive the incoming fractions.
734	! To be able to handle this case with least amount of code change, we assign the index
735	! of PFT between youngest and second-oldes (i.e., indagec_tree etc) the same as
736	! those of oldest tree PFTs (or all tree PFTs because in this cases these two indices
737	! are identical) . So that this case could be correctly handled in the subrountines
738	! of "type_conversion" and "cross_give_receive". This treatment allows use
739	! of gross land cover change subroutine with only one single age class. This single
740	! age class is "simultanously the oldest and youngest age class". At the same
741	! time, we also change the num_tree_mulagec as the same of num_crop_sinagec.
742	! The similar case also applies in grass,pasture and crop.
743
744	IF (num_tree_mulagec .EQ. 0) THEN
745	ALLOCATE(indagec_tree(num_tree_sinagec,1))
746	indagec_tree(:,1) = indall_tree(:)
747	num_tree_mulagec = num_tree_sinagec
748	ELSE
749	ALLOCATE(indagec_tree(num_tree_mulagec,nagec_tree-1))
750	END IF
751
752	IF (num_grass_mulagec .EQ. 0) THEN
753	ALLOCATE(indagec_grass(num_grass_sinagec,1))
754	indagec_grass(:,1) = indall_grass(:)
755	num_grass_mulagec = num_grass_sinagec
756	ELSE
757	ALLOCATE(indagec_grass(num_grass_mulagec,nagec_herb-1))
758	END IF
759
760	IF (num_pasture_mulagec .EQ. 0) THEN
761	ALLOCATE(indagec_pasture(num_pasture_sinagec,1))
762	indagec_pasture(:,1) = indall_pasture(:)
763	num_pasture_mulagec = num_pasture_sinagec
764	ELSE
765	ALLOCATE(indagec_pasture(num_pasture_mulagec,nagec_herb-1))
766	END IF
767
768	IF (num_crop_mulagec .EQ. 0) THEN
769	ALLOCATE(indagec_crop(num_crop_sinagec,1))
770	indagec_crop(:,1) = indall_crop(:)
771	num_crop_mulagec = num_crop_sinagec
772	ELSE
773	ALLOCATE(indagec_crop(num_crop_mulagec,nagec_herb-1))
774	END IF
775
776	! fill the non-oldest age class index arrays when number of age classes
777	! is more than 1.
778	! [chaoyuejoy@gmail.com, 2015-08-05]
779	! Note the corresponding part of code will be automatically skipped
780	! when nagec_tree ==1 and/or nagec_herb ==1, i.e., the assginment
781	! in above codes when original num_*_mulagec variables are zero will be retained.
782	itree=0
783	igrass=0
784	ipasture=0
785	icrop=0
786	DO ivma = 2,nvmap
787	staind=start_index(ivma)
788	IF (nagec_pft(ivma) > 1) THEN
789	IF (is_tree(staind)) THEN
790	itree=itree+1
791	DO j = 1,nagec_tree-1
792	indagec_tree(itree,j) = staind+nagec_tree-j-1
793	ENDDO
794	ELSE IF (natural(staind) .AND. .NOT. is_grassland_manag(staind)) THEN
795	igrass=igrass+1
796	DO j = 1,nagec_herb-1
797	indagec_grass(igrass,j) = staind+nagec_herb-j-1
798	ENDDO
799	ELSE IF (is_grassland_manag(staind)) THEN
800	ipasture=ipasture+1
801	DO j = 1,nagec_herb-1
802	indagec_pasture(ipasture,j) = staind+nagec_herb-j-1
803	ENDDO
804	ELSE
805	icrop=icrop+1
806	DO j = 1,nagec_herb-1
807	indagec_crop(icrop,j) = staind+nagec_herb-j-1
808	ENDDO
809	ENDIF
810	ENDIF
811	ENDDO
812
813
814	! we make copies of original input veget_max
815	! veget_max will be modified through different operations in order to
816	! check various purposes, e.g., whether input glcc is compatible with
817	! existing veget_max and how to allocate it etc.
818	! veget_max_old will not be modified
819	veget_max(:,:) = veget_max_org(:,:)
820	veget_max_old(:,:) = veget_max_org(:,:)
821
822	!! 2. Calcuate the fractions covered by tree, grass, pasture and crops
823	!! for each age class
824
825	!************************************************************************!
826	!**block to calculate fractions for basic veg types and age classes *!
827	! Note:
828	! 1. "calc_cover" subroutine does not depend on how many age classes
829	! there are in each MTC.
830	! 2. Fraction of baresoil is excluded here. This means transformation
831	! of baresoil to a vegetated PFT is excluded in gross land cover change.
832	veget_mtc(:,:) = 0.
833	vegagec_tree(:,:) = 0.
834	vegagec_grass(:,:) = 0.
835	vegagec_pasture(:,:) = 0.
836	vegagec_crop(:,:) = 0.
837
838	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
839	vegagec_pasture,vegagec_crop)
840	! In following call of calc_cover, veget_mtc will be updated each time,
841	! but we don't want this, so we put its initial value into veget_mtc_begin
842	! in order to retrieve this initial value later.
843	veget_mtc_begin = veget_mtc
844
845	veget_tree(:) = SUM(vegagec_tree(:,:),DIM=2)
846	veget_grass(:) = SUM(vegagec_grass(:,:),DIM=2)
847	veget_pasture(:) = SUM(vegagec_pasture(:,:),DIM=2)
848	veget_crop(:) = SUM(vegagec_crop(:,:),DIM=2)
849
850	!**end block to calculate fractions for basic veg types and age classes *!
851	!****************************************************************************!
852
853	!******************** block to handle forestry harvest **************
854	!! 2B. Here we handle the forestry wood harvest
855	! Rules:
856	! 1. We take first from second oldest forest, then oldest forest
857
858	pf2yf=1 !primary to young forest conversion because of harvest
859	sf2yf=2 !old secondary to young forest conversion because of harvest
860
861	!! Note that Deficit_pf2yf and Deficit_sf2yf are temporary, intermediate
862	!! variables. The final deficits after mutual compensation are stored in
863	!! Deficit_pf2yf_final and Deficit_sf2yf_final.
864	Deficit_pf2yf(:) = zero
865	Deficit_sf2yf(:) = zero
866	Deficit_pf2yf_final(:) = zero
867	Deficit_sf2yf_final(:) = zero
868
869	!! Note that both Surplus_pf2yf and Surplus_sf2yf and temporary intermediate
870	!! variables, the final surplus after mutual compensation are not outputed.
871	Surplus_pf2yf(:) = zero
872	Surplus_sf2yf(:) = zero
873
874	!! Note in the naming of pf2yf_compen_sf2yf and sf2yf_compen_pf2yf, active
875	!! tense is used.
876	pf2yf_compen_sf2yf(:) = zero !primary->young conversion that compensates
877	!the secondary->young conversion because of deficit
878	!in the latter
879	sf2yf_compen_pf2yf(:) = zero !seondary->young conversion that compensates
880	!the primary->young conversion because of the deficit
881	!in the latter
882
883
884	!! Define the "real" harvest matrix after considering the mutual compenstation
885	!! between primary->young and secondary->young transitions.
886	HmatrixReal(:,:) = zero !Harvest matrix real, used to hold the
887	!harvest matrix after considering the mutual
888	!compensation between primary and old secondary
889	!forest
890
891	! we sum together harvest from primary and secondary forest and consider
892	! as all happening on parimary forest.
893	HmatrixReal(:,1) = harvest_matrix(:,pf2yf) + harvest_matrix(:,sf2yf)
894
895	! Check the availability of forest fractions for harvest
896	WHERE (veget_tree(:) .LE. HmatrixReal(:,1))
897	Deficit_pf2yf_final(:) = veget_tree(:)-HmatrixReal(:,1)
898	HmatrixReal(:,1) = veget_tree(:)
899	ENDWHERE
900
901	glccRemain(:,:) = HmatrixReal(:,:)
902	glcc_pft(:,:) = 0.
903	glcc_pft_tmp(:,:) = 0.
904	glcc_pftmtc(:,:,:) = 0.
905
906	!! Allocate harvest-caused out-going primary and secondary forest fraction
907	!! into different primary and secondary forest PFTs.
908	! [Note: here we need only glcc_pft, but not glcc_pft_tmp and glcc_pftmtc.
909	! The latter two variables will be set to zero again when handling LCC in
910	! later sections.]
911	DO ipts=1,npts
912	!pf2yf
913	CALL type_conversion(ipts,pf2yf,HmatrixReal,veget_mtc,newvegfrac, &
914	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
915	1,nagec_herb, &
916	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
917	glccRemain)
918	ENDDO
919
920	! Because we use the container of type_conversion, now the glcc_pft_tmp
921	! and glcc_pftmtc have wrong information (because harvest loss is assigned
922	! on the newly created youngest-age-class pasture/crop MTCs). So they have
923	! to be re-initialized to zero. Only the information in glcc_pft is what
924	! we need.
925	glcc_pft_tmp(:,:) = 0.
926	glcc_pftmtc(:,:,:) = 0.
927	!Here we need to put glcc_pft into glcc_pftmtc for forestry harvest.
928	!The same MTC will be maintained when forest is harvested.
929	DO ivm =1,nvm
930	IF (is_tree(ivm)) THEN
931	glcc_pftmtc(:,ivm,pft_to_mtc(ivm)) = glcc_pft(:,ivm)
932	ENDIF
933	ENDDO
934	!**************** end block to handle forestry harvest **************
935	veget_max_tmp(:,:) = veget_max(:,:)
936
937
938	!************************************************************************!
939	!**block to calculate fractions for basic veg types and age classes *!
940	! Note:
941	! 1. "calc_cover" subroutine does not depend on how many age classes
942	! there are in each MTC.
943	! 2. Fraction of baresoil is excluded here. This means transformation
944	! of baresoil to a vegetated PFT is excluded in gross land cover change.
945	veget_mtc(:,:) = 0.
946	vegagec_tree(:,:) = 0.
947	vegagec_grass(:,:) = 0.
948	vegagec_pasture(:,:) = 0.
949	vegagec_crop(:,:) = 0.
950
951
952	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
953	vegagec_pasture,vegagec_crop)
954	veget_mtc = veget_mtc_begin
955
956	veget_tree(:) = SUM(vegagec_tree(:,:),DIM=2)
957	veget_grass(:) = SUM(vegagec_grass(:,:),DIM=2)
958	veget_pasture(:) = SUM(vegagec_pasture(:,:),DIM=2)
959	veget_crop(:) = SUM(vegagec_crop(:,:),DIM=2)
960	itree=1
961	igrass=2
962	ipasture=3
963	icrop=4
964	veget_4veg(:,itree) = veget_tree(:)
965	veget_4veg(:,igrass) = veget_grass(:)
966	veget_4veg(:,ipasture) = veget_pasture(:)
967	veget_4veg(:,icrop) = veget_crop(:)
968	!**end block to calculate fractions for basic veg types and age classes *!
969	!****************************************************************************!
970
971	!! 3. Decompose the LCC matrix to different PFTs
972	!! We do this through several steps:
973	! 3.1 Check whether input LCC matrix is feasible with current PFT fractions
974	! (i.e., the fractions of forest,grass,pasture and crops)
975	! and if not, adjust the transfer matrix by compensating the deficits
976	! using the surpluses.
977	! 3.2 Allocate the decreasing fractions of tree/grass/pasture/crop to their
978	! respective age classes, in the sequences of old->young.
979	! 3.3 Allocate the incoming fractions of tree/grass/pasture/crop to their
980	! respective youngest age classes. The incoming fractions are distributed
981	! according to the existing fractions of youngest-age-class PFTs of the
982	! same receiving vegetation type. If none of them exists, the incoming
983	! fraction is distributed equally.
984
985	!! 3.1 Adjust LCC matrix if it's not feasible with current PFT fractions
986
987	!++code freezing++
988	!codes below handle the mutual compenstation of transition matrices
989	!among different land cover types. This is desgined for consistency
990	!with activated DGVM.
991
992	! glcc(:,:) = glccSecondShift+glccPrimaryShift+glccNetLCC
993	! glccReal(:,:) = 0.
994	! glccDef(:,:) = 0.
995
996	! !to crop - sequence: p2c,g2c,f2c
997	! CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
998	! p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
999	! IncreDeficit)
1000
1001	! !to pasture - sequence: g2p,c2p,f2p
1002	! CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
1003	! g2p,igrass,c2p,icrop,f2p,itree,ipasture, &
1004	! IncreDeficit)
1005
1006	! !to grass - sequence: p2g,c2g,f2g
1007	! CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
1008	! p2g,ipasture,c2g,icrop,f2g,itree,igrass, &
1009	! IncreDeficit)
1010
1011	! !to forest - sequence: c2f,p2f,g2f
1012	! CALL glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
1013	! c2f,icrop,p2f,ipasture,g2f,igrass,itree, &
1014	! IncreDeficit)
1015
1016	! !! 3.2 & 3.3 Allocate LCC matrix to different PFTs/age-classes
1017
1018	! ! because we use veget_max as a proxy variable and it has been changed
1019	! ! when we derive the glccReal, so here we have to recover its original
1020	! ! values, which is veget_max_tmp after the forestry harvest.
1021	! veget_max(:,:) = veget_max_tmp(:,:)
1022
1023	! ! Calculate again fractions for different age-classes.
1024	! veget_mtc(:,:) = 0.
1025	! vegagec_tree(:,:) = 0.
1026	! vegagec_grass(:,:) = 0.
1027	! vegagec_pasture(:,:) = 0.
1028	! vegagec_crop(:,:) = 0.
1029
1030	! CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
1031	! vegagec_pasture,vegagec_crop)
1032
1033
1034	!++end codes freezing ++
1035
1036	IncreDeficit(:,:) = 0.
1037	glcc(:,:) = glccSecondShift+glccPrimaryShift+glccNetLCC
1038	glccReal(:,:) = glcc(:,:)
1039	glccRemain(:,:) = glcc(:,:)
1040
1041	! We allocate in the sequences of old->young. Within the same age-class
1042	! group, we allocate in proportion with existing PFT fractions.
1043	DO ipts=1,npts
1044	!f2c
1045	CALL type_conversion(ipts,f2c,glccReal,veget_mtc,newvegfrac, &
1046	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1047	nagec_tree,nagec_herb, &
1048	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1049	glccRemain)
1050	!f2p
1051	CALL type_conversion(ipts,f2p,glccReal,veget_mtc,newvegfrac, &
1052	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec, &
1053	nagec_tree,nagec_herb, &
1054	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1055	glccRemain)
1056	!f2g
1057	CALL type_conversion(ipts,f2g,glccReal,veget_mtc,newvegfrac, &
1058	indold_tree,indagec_tree,indagec_grass,num_grass_mulagec, &
1059	nagec_tree,nagec_herb, &
1060	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1061	glccRemain)
1062	!g2c
1063	CALL type_conversion(ipts,g2c,glccReal,veget_mtc,newvegfrac, &
1064	indold_grass,indagec_grass,indagec_crop,num_crop_mulagec, &
1065	nagec_herb,nagec_herb, &
1066	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1067	glccRemain)
1068	!g2p
1069	CALL type_conversion(ipts,g2p,glccReal,veget_mtc,newvegfrac, &
1070	indold_grass,indagec_grass,indagec_pasture,num_pasture_mulagec, &
1071	nagec_herb,nagec_herb, &
1072	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1073	glccRemain)
1074	!g2f
1075	CALL type_conversion(ipts,g2f,glccReal,veget_mtc,newvegfrac, &
1076	indold_grass,indagec_grass,indagec_tree,num_tree_mulagec, &
1077	nagec_herb,nagec_tree, &
1078	vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1079	glccRemain)
1080	!p2c
1081	CALL type_conversion(ipts,p2c,glccReal,veget_mtc,newvegfrac, &
1082	indold_pasture,indagec_pasture,indagec_crop,num_crop_mulagec, &
1083	nagec_herb,nagec_herb, &
1084	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1085	glccRemain)
1086	!p2g
1087	CALL type_conversion(ipts,p2g,glccReal,veget_mtc,newvegfrac, &
1088	indold_pasture,indagec_pasture,indagec_grass,num_grass_mulagec, &
1089	nagec_herb,nagec_herb, &
1090	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1091	glccRemain)
1092	!p2f
1093	CALL type_conversion(ipts,p2f,glccReal,veget_mtc,newvegfrac, &
1094	indold_pasture,indagec_pasture,indagec_tree,num_tree_mulagec, &
1095	nagec_herb,nagec_tree, &
1096	vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1097	glccRemain)
1098	!c2p
1099	CALL type_conversion(ipts,c2p,glccReal,veget_mtc,newvegfrac, &
1100	indold_crop,indagec_crop,indagec_pasture,num_pasture_mulagec, &
1101	nagec_herb,nagec_herb, &
1102	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1103	glccRemain)
1104	!c2g
1105	CALL type_conversion(ipts,c2g,glccReal,veget_mtc,newvegfrac, &
1106	indold_crop,indagec_crop,indagec_grass,num_grass_mulagec, &
1107	nagec_herb,nagec_herb, &
1108	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1109	glccRemain)
1110	!c2f
1111	CALL type_conversion(ipts,c2f,glccReal,veget_mtc,newvegfrac, &
1112	indold_crop,indagec_crop,indagec_tree,num_tree_mulagec, &
1113	nagec_herb,nagec_tree, &
1114	vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
1115	glccRemain)
1116	ENDDO
1117
1118	WHERE (glccRemain .GT. zero)
1119	glccReal = glcc - glccRemain
1120	IncreDeficit = -1 * glccRemain
1121	ENDWHERE
1122
1123	END SUBROUTINE glcc_SinAgeC_firstday
1124
1125
1126
1127	! ================================================================================================================================
1128	!! SUBROUTINE : type_conversion
1129	!>\BRIEF : Allocate outgoing into different age classes and incoming into
1130	!! yongest age-class of receiving MTCs.
1131	!!
1132	!! REMARK : The current dummy variables give an example of converting forests
1133	!! to crops.
1134	!! \n
1135	!_ ================================================================================================================================
1136	SUBROUTINE type_conversion(ipts,f2c,glccReal,veget_mtc,newvegfrac, &
1137	indold_tree,indagec_tree,indagec_crop,num_crop_mulagec, &
1138	nagec_giving,nagec_receive, &
1139	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
1140	glccRemain, &
1141	iagec_start)
1142
1143	IMPLICIT NONE
1144
1145	!! Input variables
1146	INTEGER, INTENT(in) :: ipts,f2c
1147	REAL(r_std), DIMENSION(:,:), INTENT(in) :: glccReal !! The "real" glcc matrix that we apply in the model
1148	!! after considering the consistency between presribed
1149	!! glcc matrix and existing vegetation fractions.
1150	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
1151	REAL(r_std), DIMENSION(:,:),INTENT(in) :: newvegfrac !!
1152	INTEGER, DIMENSION(:), INTENT(in) :: indold_tree !! Indices for PFTs giving out fractions;
1153	!! here use old tree cohort as an example
1154	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_tree !! Indices for PFTs giving out fractions;
1155	!! here use old tree cohort as an example
1156	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_crop !! Indices for secondary basic-vegetation cohorts; The youngest age classes
1157	!! of these vegetations are going to receive fractions.
1158	!! here we use crop cohorts as an example
1159	INTEGER, INTENT(in) :: num_crop_mulagec !! number of crop MTCs with more than one age classes
1160	INTEGER, INTENT(in) :: nagec_giving !! number of age classes in the giving basic types
1161	!! (i.e., tree, grass, pasture, crop), here we can use tree
1162	!! as an example, nagec=nagec_tree
1163	INTEGER, INTENT(in) :: nagec_receive !! number of age classes in the receiving basic types
1164	!! (i.e., tree, grass, pasture, crop), here we can use crop
1165	!! as an example, nagec=nagec_herb
1166	INTEGER, OPTIONAL, INTENT(in) :: iagec_start !! starting index for iagec, this is added in order to handle
1167	!! the case of secondary forest harvest.
1168
1169	!! 1. Modified variables
1170	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
1171	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1172	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fractions each PFT is going to lose
1173	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
1174	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft_tmp !! Loss of fraction in each PFT
1175	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glccRemain !! The remaining glcc matrix after applying the conversion. I.e., it will
1176	!! record the remaining unrealized transition fraction in case the donation
1177	!! vegetation is not enough compared with prescribed transition fraction.
1178	!! This variable should be initialized the same as glccReal before it's fed
1179	!! into this function.
1180
1181	!! Local vriables
1182	INTEGER :: j,iagec,iagec_start_proxy
1183	REAL(r_std) :: frac_begin,frac_used
1184	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
1185	IF (.NOT. PRESENT(iagec_start)) THEN
1186	iagec_start_proxy=1
1187	ELSE
1188	iagec_start_proxy=iagec_start
1189	ENDIF
1190
1191	! This subroutine handles the conversion from one basic-vegetation type
1192	! to another, by calling the subroutine cross_give_receive, which handles
1193	! allocation of giving-receiving fraction among the giving age classes
1194	! and receiving basic-vegetation young age classes.
1195	! We allocate in the sequences of old->young. Within the same age-class
1196	! group, we allocate in proportion with existing PFT fractions. The same
1197	! also applies in the receiving youngest-age-class PFTs, i.e., the receiving
1198	! total fraction is allocated according to existing fractions of
1199	! MTCs of the same basic vegetation type, otherwise it will be equally
1200	! distributed.
1201
1202	frac_begin = glccReal(ipts,f2c)
1203	!DO WHILE (frac_begin>min_stomate)
1204	DO iagec=iagec_start_proxy,nagec_giving
1205	IF (vegagec_tree(ipts,iagec)>frac_begin) THEN
1206	frac_used = frac_begin
1207	ELSE IF (vegagec_tree(ipts,iagec)>min_stomate) THEN
1208	frac_used = vegagec_tree(ipts,iagec)
1209	ELSE
1210	frac_used = 0.
1211	ENDIF
1212
1213	IF (frac_used>min_stomate) THEN
1214	IF (iagec==1) THEN
1215	! Note that vegagec_tree is fractions of tree age-class groups in the
1216	! the sequence of old->young, so iagec==1 means that we're handling
1217	! first the oldest-age-group tree PFTs.
1218	CALL cross_give_receive(ipts,frac_used,veget_mtc,newvegfrac, &
1219	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
1220	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1221	ELSE
1222	! Note also the sequence of indagec_tree is from old->young, so by
1223	! increasing iagec, we're handling progressively the old to young
1224	! tree age-class PFTs.
1225	CALL cross_give_receive(ipts,frac_used,veget_mtc,newvegfrac, &
1226	indagec_tree(:,iagec-1),indagec_crop,nagec_receive,num_crop_mulagec, &
1227	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
1228	ENDIF
1229	frac_begin = frac_begin-frac_used
1230	vegagec_tree(ipts,iagec)=vegagec_tree(ipts,iagec)-frac_used
1231	glccRemain(ipts,f2c) = glccRemain(ipts,f2c) - frac_used
1232	ENDIF
1233	ENDDO
1234	!ENDDO
1235
1236	END SUBROUTINE type_conversion
1237
1238	END MODULE stomate_glcchange_SinAgeC

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

Context Navigation

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

Download in other formats: