Context Navigation

stomate_fharvest_SinAgeC.f90 @ 6940

Last change on this file since 6940 was 6940, checked in by jinfeng.chang, 4 years ago
add missing files for ORCHIDEE-GMv3.2
File size: 62.1 KB

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

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

Context Navigation

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

Download in other formats: