Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

stomate_glcchange_SinAgeC.f90 @ 6737

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

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

Context Navigation

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

Download in other formats: