Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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

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

Context Navigation

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

Download in other formats: