Context Navigation

stomate_fharvest_MulAgeC.f90 @ 6940

Last change on this file since 6940 was 6940, checked in by jinfeng.chang, 4 years ago
add missing files for ORCHIDEE-GMv3.2
File size: 51.3 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_MulAgeC
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_MulAgeC
37	USE stomate_gluc_constants
38	USE xios_orchidee
39	#ifdef CPP_PARA
40	USE mpi
41	#endif
42
43	IMPLICIT NONE
44
45	PRIVATE
46	PUBLIC fharvest_MulAgeC, Get_harvest_matrix
47
48	CONTAINS
49
50	! ================================================================================================================================
51	!! SUBROUTINE gross_lcchange
52	!!
53	!>\BRIEF : Apply gross land cover change.
54	!!
55	!>\DESCRIPTION
56	!_ ================================================================================================================================
57	SUBROUTINE fharvest_MulAgeC (npts, dt_days, harvest_matrix,newvegfrac, &
58	fuelfrac, &
59	def_fuel_1hr_remain, def_fuel_10hr_remain, &
60	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
61	deforest_litter_remain, deforest_biomass_remain, &
62	convflux, &
63	glcc_pft, glcc_pftmtc, &
64	veget_max, prod10, prod100, &
65	PFTpresent, senescence, moiavail_month, moiavail_week, &
66	gpp_week, ngd_minus5, resp_maint, resp_growth, &
67	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
68	ind, lm_lastyearmax, everywhere, age, &
69	co2_to_bm, gpp_daily, co2_fire, &
70	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
71	gdd_m5_dormance, ncd_dormance, &
72	lignin_struc, carbon, leaf_frac, &
73	deepC_a, deepC_s, deepC_p, &
74	leaf_age, bm_to_litter, biomass, litter, &
75	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
76
77	IMPLICIT NONE
78
79	!! 0.1 Input variables
80
81	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
82	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
83	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: harvest_matrix !!
84	REAL(r_std), DIMENSION (npts),INTENT(in) :: fuelfrac !!
85	REAL(r_std), DIMENSION (npts,nvmap), INTENT(in) :: newvegfrac
86	!!
87
88	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1hr_remain
89	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_10hr_remain
90	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_100hr_remain
91	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in) :: def_fuel_1000hr_remain
92	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
93	!! deforestation region.
94	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
95	!! deforestation region.
96
97
98	!! 0.2 Output variables
99	REAL(r_std), DIMENSION(npts,nwp), INTENT(inout) :: convflux !! release during first year following land cover
100	!! change
101	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
102	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout):: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
103	!! i.e., the contribution of each PFT to the youngest age-class of MTC
104
105	!! 0.3 Modified variables
106	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT (LAI ->
107	!! infinity) on ground (unitless)
108	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
109	!! pool after the annual release for each
110	!! compartment (10 + 1 : input from year of land
111	!! cover change)
112	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
113	!! pool after the annual release for each
114	!! compartment (100 + 1 : input from year of land
115	!! cover change)
116	!! pool compartments
117	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
118	!! each pixel
119	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
120	!! for deciduous trees)
121	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
122	!! unitless)
123	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
124	!! (0 to 1, unitless)
125	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
126	!! @tex $(gC m^{-2} day^{-1})$ @endtex
127	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
128	!! -5 deg C (for phenology)
129	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
130	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
131	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
132	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
133	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
134	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
135	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
136	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
137	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
138	!! the growing season (days)
139	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
140	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
141	!! @tex $(m^{-2})$ @endtex
142	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
143	!! @tex ($gC m^{-2}$) @endtex
144	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
145	!! very localized (after its introduction) (?)
146	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
147	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
148	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
149	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
150	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
151	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_fire !! Fire carbon emissions
152	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
153	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
154	!! availability (days)
155	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
156	!! (for phenology) - this is written to the
157	!! history files
158	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
159	!! for crops
160	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
161	!! C (for phenology)
162	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
163	!! leaves were lost (for phenology)
164	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
165	!! above and below ground
166	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
167	!! @tex ($gC m^{-2}$) @endtex
168	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
169	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
170	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
171	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
172	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
173	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
174	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
175	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
176	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
177	!! below ground @tex ($gC m^{-2}$) @endtex
178	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
179	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
180	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
181	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
182
183	!! 0.4 Local variables
184	REAL(r_std), DIMENSION(nvmap,nparts,nelements) :: bm_to_litter_pro !! conversion of biomass to litter
185	!! @tex ($gC m^{-2} day^{-1}$) @endtex
186	REAL(r_std), DIMENSION(nvmap,nparts,nelements) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
187	REAL(r_std), DIMENSION(nvmap) :: veget_max_pro !! "maximal" coverage fraction of a PFT (LAI ->
188	!! infinity) on ground (unitless)
189	REAL(r_std), DIMENSION(nvmap,ncarb) :: carbon_pro !! carbon pool: active, slow, or passive
190	!! @tex ($gC m^{-2}$) @endtex
191	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_a_pro !! Permafrost carbon pool: active, slow, or passive
192	!! @tex ($gC m^{-3}$) @endtex
193	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_s_pro !! Permafrost carbon pool: active, slow, or passive
194	!! @tex ($gC m^{-3}$) @endtex
195	REAL(r_std), DIMENSION(nvmap,ndeep) :: deepC_p_pro !! Permafrost carbon pool: active, slow, or passive
196	!! @tex ($gC m^{-3}$) @endtex
197	REAL(r_std), DIMENSION(nvmap,nlitt,nlevs,nelements) :: litter_pro !! metabolic and structural litter, above and
198	!! below ground @tex ($gC m^{-2}$) @endtex
199	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_1hr_pro
200	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_10hr_pro
201	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_100hr_pro
202	REAL(r_std), DIMENSION(nvmap,nlitt,nelements) :: fuel_1000hr_pro
203	REAL(r_std), DIMENSION(nvmap,nlevs) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
204	!! above and below ground
205	REAL(r_std), DIMENSION(nvmap,nleafages) :: leaf_frac_pro !! fraction of leaves in leaf age class
206	REAL(r_std), DIMENSION(nvmap,nleafages) :: leaf_age_pro !! fraction of leaves in leaf age class
207	LOGICAL, DIMENSION(nvmap) :: PFTpresent_pro, senescence_pro !! Is pft there (unitless)
208	REAL(r_std), DIMENSION(nvmap) :: ind_pro, age_pro, lm_lastyearmax_pro, npp_longterm_pro
209	REAL(r_std), DIMENSION(nvmap) :: everywhere_pro
210	REAL(r_std), DIMENSION(nvmap) :: gpp_daily_pro, npp_daily_pro, co2_to_bm_pro
211	REAL(r_std), DIMENSION(nvmap) :: resp_maint_pro, resp_growth_pro
212	REAL(r_std), DIMENSION(nvmap) :: resp_hetero_pro, co2_fire_pro
213
214	INTEGER :: ipts,ivm,ivma,l,m,ipft_young_agec
215	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
216
217	REAL(r_std), DIMENSION(npts,nvmap) :: glcc_mtc !! Increase in fraction of each MTC in its youngest age-class
218	REAL(r_std), DIMENSION(npts,nvm) :: glccReal_tmp !! A temporary variable to hold glccReal
219	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf_final !!
220	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf_final !!
221	REAL(r_std), DIMENSION(npts) :: pf2yf_compen_sf2yf !!
222	REAL(r_std), DIMENSION(npts) :: sf2yf_compen_pf2yf !!
223	REAL(r_std), DIMENSION(npts,nvm) :: glcc_harvest !! Loss of fraction due to forestry harvest
224	REAL(r_std), DIMENSION(npts) :: harvest_wood !! Loss of fraction due to forestry harvest
225
226	WRITE(numout,*) 'Entering fharvest_MulAgeC'
227	glcc_harvest(:,:) = zero
228	harvest_wood(:) = zero
229	glccReal_tmp(:,:) = zero
230
231	CALL MulAgeC_fh_firstday(npts,veget_max,newvegfrac,harvest_matrix, &
232	glcc_pft,glcc_pftmtc, &
233	Deficit_pf2yf_final, Deficit_sf2yf_final, &
234	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
235
236	glcc_mtc(:,:) = SUM(glcc_pftmtc,DIM=2)
237
238	! Write out forestry harvest variables
239	DO ipts = 1,npts
240	DO ivm = 1,nvm
241	DO ivma = 1,nvmap
242	IF (is_tree(ivm) .AND. is_tree(start_index(ivma))) THEN
243	glcc_harvest(ipts,ivm) = glcc_harvest(ipts,ivm) + glcc_pftmtc(ipts,ivm,ivma)
244	ENDIF
245	ENDDO
246	harvest_wood(ipts) = harvest_wood(ipts)+ (biomass(ipts,ivm,isapabove,icarbon)+ &
247	biomass(ipts,ivm,iheartabove,icarbon))*glcc_harvest(ipts,ivm)
248	ENDDO
249	ENDDO
250
251	DO ipts=1,npts
252
253	!! Initialize the _pro variables
254	bm_to_litter_pro(:,:,:)=zero
255	biomass_pro(:,:,:)=zero
256	veget_max_pro(:)=zero
257	carbon_pro(:,:)=zero
258	deepC_a_pro(:,:)=zero
259	deepC_s_pro(:,:)=zero
260	deepC_p_pro(:,:)=zero
261	litter_pro(:,:,:,:)=zero
262	fuel_1hr_pro(:,:,:)=zero
263	fuel_10hr_pro(:,:,:)=zero
264	fuel_100hr_pro(:,:,:)=zero
265	fuel_1000hr_pro(:,:,:)=zero
266	lignin_struc_pro(:,:)=zero
267
268	leaf_frac_pro = zero
269	leaf_age_pro = zero
270	PFTpresent_pro(:) = .FALSE.
271	senescence_pro(:) = .TRUE.
272	ind_pro = zero
273	age_pro = zero
274	lm_lastyearmax_pro = zero
275	npp_longterm_pro = zero
276	everywhere_pro = zero
277
278	gpp_daily_pro=zero
279	npp_daily_pro=zero
280	co2_to_bm_pro=zero
281	resp_maint_pro=zero
282	resp_growth_pro=zero
283	resp_hetero_pro=zero
284	co2_fire_pro=zero
285
286	! Note that we assume people don't intentionally change baresoil to
287	! vegetated land.
288	DO ivma = 2,nvmap
289	! we assume only the youngest age class receives the incoming PFT
290	! [chaoyuejoy@gmail.com 2015-08-04] This line is commented to allow
291	! the case of only single age class being handled.
292
293	! here we set glcc_mtc(ipts,ivma) > min_stomate as a condition,
294	! this is necessary because later on in the subroutine of
295	! add_incoming_proxy_pft we have to merge the newly established
296	! youngest proxy with potentially exisiting youngest receiving MTC,
297	! thus have to devide a new fraction of (frac_proxy + frac_exist),
298	! but in case frac_exist = zero, we risk deviding by a very small value
299	! of frac_proxy and thus we want it to be bigger than min_stomate.
300	IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
301
302	! 1. we accumulate the scalar variables that will be inherited
303	! note we don't handle the case of harvesting forest because
304	! we assume glcc_pftmtc(forest->forest) would be zero and this
305	! case won't occur as it's filtered by the condition of
306	! (frac>min_stomate)
307	CALL collect_legacy_pft_forestry(npts, ipts, ivma, glcc_pftmtc, &
308	fuelfrac, &
309	biomass, bm_to_litter, carbon, litter, &
310	deepC_a, deepC_s, deepC_p, &
311	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
312	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
313	resp_maint, resp_growth, resp_hetero, co2_fire, &
314	def_fuel_1hr_remain, def_fuel_10hr_remain, &
315	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
316	deforest_litter_remain, deforest_biomass_remain, &
317	veget_max_pro(ivma), carbon_pro(ivma,:), &
318	lignin_struc_pro(ivma,:), litter_pro(ivma,:,:,:), &
319	deepC_a_pro(ivma,:), deepC_s_pro(ivma,:), deepC_p_pro(ivma,:), &
320	fuel_1hr_pro(ivma,:,:), fuel_10hr_pro(ivma,:,:), &
321	fuel_100hr_pro(ivma,:,:), fuel_1000hr_pro(ivma,:,:), &
322	bm_to_litter_pro(ivma,:,:), co2_to_bm_pro(ivma), &
323	gpp_daily_pro(ivma), npp_daily_pro(ivma), &
324	resp_maint_pro(ivma), resp_growth_pro(ivma), &
325	resp_hetero_pro(ivma), co2_fire_pro(ivma), &
326	convflux,prod10,prod100)
327
328	!++TEMP++
329	! Here we substract the outgoing fraction from the source PFT.
330	! If a too small fraction remains in this source PFT, then it is
331	! exhausted, we empty it. The subroutine 'empty_pft' might be
332	! combined with 'collect_legacy_pft', but now we just put it here.
333	DO ivm = 1,nvm
334	IF( glcc_pftmtc(ipts,ivm,ivma)>zero ) THEN
335	veget_max(ipts,ivm) = veget_max(ipts,ivm)-glcc_pftmtc(ipts,ivm,ivma)
336	IF ( veget_max(ipts,ivm)<min_stomate ) THEN
337	CALL empty_pft(ipts, ivm, veget_max, biomass, ind, &
338	carbon, litter, lignin_struc, bm_to_litter, &
339	deepC_a, deepC_s, deepC_p, &
340	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
341	gpp_daily, npp_daily, gpp_week, npp_longterm, &
342	co2_to_bm, resp_maint, resp_growth, resp_hetero, &
343	lm_lastyearmax, leaf_frac, leaf_age, age, &
344	everywhere, PFTpresent, when_growthinit, &
345	senescence, gdd_from_growthinit, gdd_midwinter, &
346	time_hum_min, gdd_m5_dormance, ncd_dormance, &
347	moiavail_month, moiavail_week, ngd_minus5)
348	ENDIF
349	ENDIF
350	ENDDO
351
352	ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
353	ENDDO !(DO ivma = 2,nvmap)
354
355	! We can only establish new youngest proxy and add it to the
356	! existing youngest-age PFT after all the harvest is done, to
357	! avoid the dilution of harvestable biomass by the young proxy
358	! and ensure consistency. Therefore now we have to loop again
359	! over nvmap.
360	DO ivma = 2,nvmap
361	IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
362
363	ipft_young_agec = start_index(ivma)
364
365	! 2. we establish a proxy PFT with the fraction of veget_max_pro,
366	! which is going to be either merged with existing target
367	! `ipft_young_agec` PFT, or fill the place if no existing target PFT
368	! exits.
369	CALL initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro(ivma), &
370	biomass_pro(ivma,:,:), co2_to_bm_pro(ivma), ind_pro(ivma), &
371	age_pro(ivma), &
372	senescence_pro(ivma), PFTpresent_pro(ivma), &
373	lm_lastyearmax_pro(ivma), everywhere_pro(ivma), &
374	npp_longterm_pro(ivma), &
375	leaf_frac_pro(ivma,:),leaf_age_pro(ivma,:))
376
377	CALL sap_take (ipts,ivma,veget_max,biomass_pro(ivma,:,:), &
378	biomass,co2_to_bm_pro(ivma))
379
380	! 3. we merge the newly initiazlized proxy PFT into existing one
381	! or use it to fill an empty PFT slot.
382	CALL add_incoming_proxy_pft(npts, ipts, ipft_young_agec, veget_max_pro(ivma),&
383	carbon_pro(ivma,:), litter_pro(ivma,:,:,:), lignin_struc_pro(ivma,:), &
384	bm_to_litter_pro(ivma,:,:), &
385	deepC_a_pro(ivma,:), deepC_s_pro(ivma,:), deepC_p_pro(ivma,:), &
386	fuel_1hr_pro(ivma,:,:), fuel_10hr_pro(ivma,:,:), &
387	fuel_100hr_pro(ivma,:,:), fuel_1000hr_pro(ivma,:,:), &
388	biomass_pro(ivma,:,:), co2_to_bm_pro(ivma), &
389	npp_longterm_pro(ivma), ind_pro(ivma), &
390	lm_lastyearmax_pro(ivma), age_pro(ivma), everywhere_pro(ivma), &
391	leaf_frac_pro(ivma,:), leaf_age_pro(ivma,:), &
392	PFTpresent_pro(ivma), senescence_pro(ivma), &
393	gpp_daily_pro(ivma), npp_daily_pro(ivma), &
394	resp_maint_pro(ivma), resp_growth_pro(ivma), &
395	resp_hetero_pro(ivma), co2_fire_pro(ivma), &
396	veget_max, carbon, litter, lignin_struc, bm_to_litter, &
397	deepC_a, deepC_s, deepC_p, &
398	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
399	biomass, co2_to_bm, npp_longterm, ind, &
400	lm_lastyearmax, age, everywhere, &
401	leaf_frac, leaf_age, PFTpresent, senescence, &
402	gpp_daily, npp_daily, resp_maint, resp_growth, &
403	resp_hetero, co2_fire)
404
405	ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
406	ENDDO !(DO ivma=1,nvmap)
407
408	ENDDO !(DO ipts=1,npts)
409
410
411	CALL histwrite_p (hist_id_stomate, 'glcc_harvest', itime, &
412	glcc_harvest, npts*nvm, horipft_index)
413	CALL xios_orchidee_send_field ('glcc_harvest', glcc_harvest) ! kjpindex,nvm
414
415	CALL histwrite_p (hist_id_stomate, 'harvest_wood', itime, &
416	harvest_wood, npts, hori_index)
417	CALL xios_orchidee_send_field ('harvest_wood', harvest_wood) ! kjpindex
418
419	glccReal_tmp(:,:) = zero
420	glccReal_tmp(:,1) = Deficit_pf2yf_final
421	glccReal_tmp(:,2) = Deficit_sf2yf_final
422	glccReal_tmp(:,3) = pf2yf_compen_sf2yf
423	glccReal_tmp(:,4) = sf2yf_compen_pf2yf ! this is zero as currently the deficit
424	! in primary forest harvest is not
425	! compensated for by the surplus in
426	! secondary forest harvest.
427
428	CALL histwrite_p (hist_id_stomate, 'DefiComForHarvest', itime, &
429	glccReal_tmp, npts*nvm, horipft_index)
430	CALL xios_orchidee_send_field ('DefiComForHarvest', glccReal_tmp) ! kjpindex,nvm
431
432	END SUBROUTINE fharvest_MulAgeC
433
434	! ================================================================================================================================
435	!! SUBROUTINE : gross_lcc_firstday
436	!!
437	!>\BRIEF : When necessary, adjust input glcc matrix, and allocate it
438	!! into different contributing age classes and receiving
439	!! youngest age classes.
440	!! \n
441	!_ ================================================================================================================================
442
443	! Note: it has this name because this subroutine will also be called
444	! the first day of each year to precalculate the forest loss for the
445	! deforestation fire module.
446	SUBROUTINE MulAgeC_fh_firstday(npts,veget_max_org,newvegfrac,harvest_matrix, &
447	glcc_pft,glcc_pftmtc, &
448	Deficit_pf2yf_final, Deficit_sf2yf_final, &
449	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
450
451	IMPLICIT NONE
452
453	!! 0.1 Input variables
454
455	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
456	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_org !! "maximal" coverage fraction of a PFT on the ground
457	!! May sum to
458	!! less than unity if the pixel has
459	!! nobio area. (unitless, 0-1)
460	REAL(r_std), DIMENSION(npts,nvmap), INTENT(in) :: newvegfrac !! used to guid the allocation of new PFTs.
461	REAL(r_std), DIMENSION(npts,12),INTENT(in) :: harvest_matrix !!
462	!!
463	!! 0.2 Output variables
464	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout) :: glcc_pftmtc !! the fractions each PFT is going to lose
465	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
466
467	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Deficit_pf2yf_final !!
468	REAL(r_std), DIMENSION(npts), INTENT(inout) :: Deficit_sf2yf_final !!
469	REAL(r_std), DIMENSION(npts), INTENT(inout) :: pf2yf_compen_sf2yf !!
470	REAL(r_std), DIMENSION(npts), INTENT(inout) :: sf2yf_compen_pf2yf !!
471
472	!! 0.3 Modified variables
473
474	!! 0.4 Local variables
475	REAL(r_std), DIMENSION(npts,nvmap) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
476	REAL(r_std), DIMENSION(npts,nvmap) :: veget_mtc_begin !! "maximal" coverage fraction of a PFT on the ground
477	REAL(r_std), DIMENSION(npts,nagec_tree) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
478	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_grass !! fraction of grass age-class groups, in sequence of old->young
479	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_pasture !! fraction of pasture age-class groups, in sequence of old->young
480	REAL(r_std), DIMENSION(npts,nagec_herb) :: vegagec_crop !! fraction of crop age-class groups, in sequence of old->young
481
482
483	REAL(r_std), DIMENSION(npts,4) :: veget_4veg !! "maximal" coverage fraction of a PFT on the ground
484	REAL(r_std), DIMENSION(npts) :: veget_tree !! "maximal" coverage fraction of a PFT on the ground
485	REAL(r_std), DIMENSION(npts) :: veget_grass !! "maximal" coverage fraction of a PFT on the ground
486	REAL(r_std), DIMENSION(npts) :: veget_pasture !! "maximal" coverage fraction of a PFT on the ground
487	REAL(r_std), DIMENSION(npts) :: veget_crop !! "maximal" coverage fraction of a PFT on the ground
488
489	REAL(r_std), DIMENSION(npts,nvm) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
490	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_tmp !! "maximal" coverage fraction of a PFT on the ground
491	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_old !! "maximal" coverage fraction of a PFT on the ground
492	REAL(r_std), DIMENSION(npts,nvm) :: glcc_pft_tmp !! Loss of fraction in each PFT
493
494	REAL(r_std), DIMENSION(npts,nvm,nvmap) :: glcc_pftmtc_harvest !! a temporary variable to hold the fractions each PFT is going to lose
495	REAL(r_std), DIMENSION(npts,12) :: glccRealHarvest !! real matrix applied for forestry harvest.
496
497	REAL(r_std), DIMENSION(npts,nvm) :: glccReal_tmp !! A temporary variable to hold glccReal
498
499	LOGICAL, SAVE :: MulAgeC_fh_firstday_done = .FALSE.
500
501	INTEGER :: ivma, pf2yf, sf2yf, ivm
502
503	REAL(r_std), DIMENSION(npts,12) :: FHmatrix_remainA !!
504	REAL(r_std), DIMENSION(npts,12) :: FHmatrix_remainB !!
505	REAL(r_std), DIMENSION(npts,12) :: glccRemain !!
506
507	INTEGER :: ipts,IndStart_f,IndEnd_f
508	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
509
510	!Some more local configurations
511	LOGICAL :: allow_youngest_forest_WoodHarvest = .FALSE.
512
513
514	! we make copies of original input veget_max (which is veget_max_org
515	! in the subroutine parameter list).
516	! veget_max will be modified through different operations in order to
517	! check for various purposes, e.g., whether input harvest matrix
518	! is compatible with existing veget_max and how to allocate it etc.
519	! veget_max_old will not be modified
520	veget_max(:,:) = veget_max_org(:,:)
521	veget_max_old(:,:) = veget_max_org(:,:)
522
523	!******************** block to handle forestry harvest **************
524	!! 1. Handle the forestry harvest process
525
526	!! 1.0 Some preparation
527
528	pf2yf=1 !primary to young forest conversion because of harvest
529	sf2yf=2 !old secondary to young forest conversion because of harvest
530
531	Deficit_pf2yf_final(:) = zero
532	Deficit_sf2yf_final(:) = zero
533
534	! Note in the naming of pf2yf_compen_sf2yf and sf2yf_compen_pf2yf, active
535	! tense is used. I.e., pf2yf_compen_sf2yf means the fraction which pf2yf
536	! compenstates for sf2yf
537	pf2yf_compen_sf2yf(:) = zero !primary->young conversion that compensates
538	!the secondary->young conversion because of deficit
539	!in the latter
540	sf2yf_compen_pf2yf(:) = zero !seondary->young conversion that compensates
541	!the primary->young conversion because of the deficit
542	!in the latter
543
544	! we now have to fill the transtion of forest->forest in the process of harvest
545	! into our target matrix glcc_pftmtc. Thus we will initiliaze them first.
546	glcc_pft(:,:) = 0.
547	glcc_pft_tmp(:,:) = 0.
548	glcc_pftmtc(:,:,:) = 0.
549	glccRemain(:,:) = harvest_matrix(:,:)
550
551	!! 2.1 Handle secondary forest harvest
552
553	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
554	vegagec_pasture,vegagec_crop)
555
556	! In following call of calc_cover, veget_mtc will be updated each time,
557	! but we don't want this, so we put its initial value into veget_mtc_begin
558	! in order to retrieve this initial value later. veget_mtc is used in the
559	! subroutine type_conversion to guild the allocation of newly created
560	! lands of a certain vegetation type (e.g., pasture) into its compoenent
561	! meta-classes (e.g. C3 and C4 pasture).
562	veget_mtc_begin = veget_mtc
563
564	! Allocate harvest-caused out-going primary and secondary forest fraction
565	! into different primary and secondary (all other younger age classes) forest PFTs.
566
567	! [Note]
568	! Below we used the tempelate of type_conversion but in fact we need
569	! only glcc_pft, which means the fraction loss in each PFT. We then need to
570	! use glcc_pft to fill glcc_pftmtc (our final target matrix), assuming that
571	! the loss of forest PFT will go to the youngest age class of its forest MTC.
572	! Therefore we assume no changes of forset species (meta-class) in forestry
573	! harvest and following tree-replanting.
574	! Although glcc_pftmtc and glcc_pft_tmp will be automatically filled when
575	! we use the tempelate `type_conversion` by calling it as below, but it makes
576	! no sense because they will be reset later.
577
578	!! 1.1 Secondary forest harvest.
579
580	! We first handle within the secondary forest age classes, in the sequence
581	! of old->young
582
583	IndStart_f = 2 ! note the indecies for tree age class are in the
584	! sequence of from old to young, thus index=2 means the
585	! 2nd oldest age class.
586	IndEnd_f = nagec_tree-1 ! nagec_tree-2: The 3rd youngest age class
587	! nagec_tree-1: The 2nd youngest age class
588	! nagec_tree: The youngest age class
589
590	IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
591	write(numout,*) 'Forest harvest: Age class index cannot be negative or zero!'
592	STOP
593	ENDIF
594
595	DO ipts=1,npts
596	!sf2yf
597	CALL type_conversion(ipts,sf2yf,harvest_matrix,veget_mtc,newvegfrac, &
598	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec,&
599	IndEnd_f,nagec_herb, &
600	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
601	glccRemain, &
602	.TRUE., iagec_start=IndStart_f)
603	ENDDO
604	FHmatrix_remainA(:,:) = glccRemain
605
606	!! 1.2 Use primary forest harvest to compensate the deficit in secondary
607	!! forest harvest. Note such compensation happens automatically here.
608
609	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
610	vegagec_pasture,vegagec_crop)
611
612	! retrieve the initial veget_mtc value, as explained above.
613	veget_mtc = veget_mtc_begin
614
615	! we check whether the required harvest of secondary forest
616	! is met by the existing secondary forest fractions. Otherwise
617	! we use the oldest-age-class forest to compenstate it.
618	DO ipts=1,npts
619	IF (FHmatrix_remainA(ipts,sf2yf) .GT. zero) THEN
620	! in this case, the existing secondary forest fraction
621	! is not enough for secondary forest harvest, we have to
622	! use primary (oldest age class) foret to compensate it.
623
624	IndStart_f = 1 ! Oldest age class
625	IndEnd_f = 1 ! Oldest age class
626
627	!sf2yf
628	CALL type_conversion(ipts,sf2yf,FHmatrix_remainA,veget_mtc,newvegfrac, &
629	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec,&
630	IndEnd_f,nagec_herb, &
631	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
632	glccRemain, &
633	.TRUE., iagec_start=IndStart_f)
634	pf2yf_compen_sf2yf(ipts) = FHmatrix_remainA(ipts,sf2yf) - glccRemain(ipts,sf2yf)
635
636	!!++Temp++
637	! Normally we don't allow youngest foret cohort to be harvested,
638	! but when external input driving data are not consistent, we harvest
639	! even the youngest forest cohort to make sure it's consistent with
640	! the input data and comparable with the case of SinAgeC. However, except
641	! for the thse two reasons, we should not harvest the youngest age class
642	! as the wood there is not feasible for usage in most cases.
643	IF (allow_youngest_forest_WoodHarvest) THEN
644	FHmatrix_remainA(ipts,sf2yf) = glccRemain(ipts,sf2yf)
645	IF (FHmatrix_remainA(ipts,sf2yf) .GT. zero) THEN
646
647	IndStart_f = nagec_tree ! youngest age class
648	IndEnd_f = nagec_tree ! youngest age class
649
650	IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
651	write(numout,*) 'Age class index cannot be negative or zero!'
652	STOP
653	ENDIF
654
655	!sf2yf
656	CALL type_conversion(ipts,sf2yf,FHmatrix_remainA,veget_mtc,newvegfrac, &
657	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec,&
658	IndEnd_f,nagec_herb, &
659	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
660	glccRemain, &
661	.TRUE., iagec_start=IndStart_f)
662	ENDIF
663	ENDIF !(IF allow_youngest_forest_WoodHarvest)
664
665	ENDIF
666	ENDDO
667	FHmatrix_remainB(:,:) = glccRemain
668
669	!! 1.3 Handle primary forest harvest
670
671	CALL calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
672	vegagec_pasture,vegagec_crop)
673	veget_mtc = veget_mtc_begin
674
675	! We check first whether there is still deficit in the required secondary
676	! harvest even after being compensated for by the primary forests.
677	! If yes, that means all existing harvestable forests
678	! are depleted, thus required primary harvest will be suppressed.
679	! Otherwise we will handle primary forest harvest starting from the
680	! oldest-age-class forest
681
682	DO ipts=1,npts
683	IF (FHmatrix_remainB(ipts,sf2yf) .GT. min_stomate) THEN
684	! in this case, all forest fraction is depleted in handling
685	! required secondary forest harvest. We thus suppress the
686	! the required primary forest harvest.
687	Deficit_sf2yf_final(ipts) = -1 * FHmatrix_remainB(ipts,sf2yf)
688	Deficit_pf2yf_final(ipts) = -1 * FHmatrix_remainB(ipts,pf2yf)
689
690
691	ELSE
692	! there are still forest can be used for required primary forest harvest.
693	! we assume primary harvest occurs in the oldest age class. Here,
694	! we will start from the oldest forest and then move to the younger ones,
695	! until the second youngest when the youngest forest cohort is not allowed
696	! for harvesting, and to youngest one when it is allowed.
697
698	IndStart_f = 1 ! Oldest age class
699
700	!!++Temp++
701	IF (allow_youngest_forest_WoodHarvest) THEN
702	IndEnd_f = nagec_tree ! youngest age class
703	ELSE
704	IndEnd_f = nagec_tree-1 ! 2nd youngest age class
705	ENDIF
706
707	IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
708	write(numout,*) 'Age class index cannot be negative or zero!'
709	STOP
710	ENDIF
711
712	!pf2yf
713	CALL type_conversion(ipts,pf2yf,FHmatrix_remainB,veget_mtc,newvegfrac, &
714	indold_tree,indagec_tree,indagec_pasture,num_pasture_mulagec,&
715	IndEnd_f,nagec_herb, &
716	vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
717	glccRemain, &
718	.TRUE., iagec_start=IndStart_f)
719	ENDIF
720
721	IF (glccRemain(ipts,pf2yf) .GT. min_stomate) THEN
722	Deficit_pf2yf_final(ipts) = -1 * glccRemain(ipts,pf2yf)
723	ENDIF
724	ENDDO
725	glccRealHarvest = harvest_matrix - glccRemain
726
727	! Because we use the template of `type_conversion, now the glcc_pft_tmp
728	! and glcc_pftmtc have wrong information (because harvest loss is assigned
729	! on the newly created youngest-age-class pasture/crop MTCs). So they have
730	! to be re-initialized to zero. Only the information in glcc_pft is what
731	! we need, as explained above.
732	glcc_pft_tmp(:,:) = 0.
733	glcc_pftmtc(:,:,:) = 0.
734	! Here we need to put glcc_pft into glcc_pftmtc for forestry harvest.
735	! The same MTC will be maintained when forest is harvested.
736	DO ivm =1,nvm
737	IF (is_tree(ivm)) THEN
738	glcc_pftmtc(:,ivm,pft_to_mtc(ivm)) = glcc_pft(:,ivm)
739	ENDIF
740	ENDDO
741	glcc_pftmtc_harvest = glcc_pftmtc
742	!**************** end block to handle forestry harvest **************
743
744	IF (.NOT. MulAgeC_fh_firstday_done) THEN
745
746	glccReal_tmp = zero
747	glccReal_tmp(:,1:12) = glccRealHarvest
748	CALL histwrite_p (hist_id_stomate, 'glccRealHarvest', itime, &
749	glccReal_tmp, npts*nvm, horipft_index)
750	CALL xios_orchidee_send_field ('glccRealHarvest', glccReal_tmp) ! kjpindex,nvm
751
752	DO ivma = 1, nvmap
753	WRITE(part_str,'(I2)') ivma
754	IF (ivma < 10) part_str(1:1) = '0'
755	CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_H_'//part_str(1:LEN_TRIM(part_str)), &
756	itime, glcc_pftmtc_harvest(:,:,ivma), npts*nvm, horipft_index)
757	ENDDO
758	CALL xios_orchidee_send_field ('glcc_pftmtc_H', glcc_pftmtc_harvest) ! kjpindex,nvm,nvmap
759
760	MulAgeC_fh_firstday_done = .TRUE.
761	ENDIF
762
763	END SUBROUTINE MulAgeC_fh_firstday
764
765
766	SUBROUTINE Get_harvest_matrix(npts,veget_max,newvegfrac,harvest_matrix, &
767	biomass, harvest_biomass, area_land_m2, &
768	glcc_pft,glcc_pftmtc, &
769	Deficit_pf2yf_final, Deficit_sf2yf_final, &
770	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
771
772	IMPLICIT NONE
773
774	!! 0.1 Input variables
775
776	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
777	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !!
778	REAL(r_std), DIMENSION(npts), INTENT(in) :: area_land_m2
779
780	REAL(r_std), DIMENSION(npts,nvmap), INTENT(in) :: newvegfrac
781	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
782
783	REAL(r_std), DIMENSION(npts,12),INTENT(in) :: harvest_biomass !! the 1st dimension is for industrial wood, 2nd dimension for fuelwood.
784
785	!! 0.2 Output variables
786	REAL(r_std), DIMENSION(npts,12),INTENT(inout) :: harvest_matrix !!
787
788
789	!! 0.3 Modified variables
790	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! Loss of fraction in each PFT
791	REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
792
793	!! 0.4 Local variables
794	REAL(r_std), DIMENSION(npts) :: mean_abwood_dens !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
795	REAL(r_std), DIMENSION(npts) :: total_abwood !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
796	REAL(r_std), DIMENSION(npts) :: total_abwood_harvest !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
797	REAL(r_std), DIMENSION(npts) :: treefrac !! "maximal" coverage fraction of a PFT on the ground
798
799	INTEGER :: ipts,ivm,ivma,num,ierr
800
801	! Note: these parameters are useless here but to conform with the form they
802	! are needed to call the subroutine MulAgeC_fh_firstday.
803	REAL(r_std), DIMENSION(npts,nvmap) :: glcc_mtc !! Increase in fraction of each MTC in its youngest age-class
804	REAL(r_std), DIMENSION(npts,nvm) :: glccReal_tmp !! A temporary variable to hold glccReal
805	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf_final !!
806	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf_final !!
807	REAL(r_std), DIMENSION(npts) :: pf2yf_compen_sf2yf !!
808	REAL(r_std), DIMENSION(npts) :: sf2yf_compen_pf2yf !!
809	REAL(r_std) :: ratio,input_total,model_total
810	REAL(r_std) :: input_total_allproc,model_total_allproc
811
812
813	input_total = SUM(harvest_biomass(:,1:2))
814
815	treefrac(:) = zero
816	total_abwood(:) = zero
817	DO ivm = 1,nvm
818	IF (is_tree(ivm) ) THEN
819	total_abwood(:) = total_abwood(:) + veget_max(:,ivm) * &
820	(biomass(:,ivm,isapabove,icarbon)+biomass(:,ivm,iheartabove,icarbon))
821	treefrac(:) = treefrac(:) + veget_max(:,ivm)
822	ENDIF
823	ENDDO
824
825	mean_abwood_dens(:) = zero
826	WHERE( treefrac(:) .GT. min_stomate )
827	mean_abwood_dens(:) = total_abwood(:)/treefrac(:)
828	ENDWHERE
829
830	! harvest_matrix(:,:) = zero
831	! WHERE( mean_abwood_dens(:) .GT. min_stomate .AND. area_land_m2(:) .GT. min_stomate )
832	! harvest_matrix(:,1) = SUM(harvest_biomass(:,1:2),DIM=2)/(mean_abwood_dens(:)*area_land_m2(:))
833	! ENDWHERE
834
835	CALL MulAgeC_fh_firstday(npts,veget_max,newvegfrac,harvest_matrix, &
836	glcc_pft,glcc_pftmtc, &
837	Deficit_pf2yf_final, Deficit_sf2yf_final, &
838	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
839
840
841	! calculate model_total
842	total_abwood_harvest(:) = zero
843	DO ivm = 1,nvm
844	IF (is_tree(ivm) ) THEN
845	total_abwood_harvest(:) = total_abwood_harvest(:) + glcc_pft(:,ivm) area_land_m2(:) &
846	(biomass(:,ivm,isapabove,icarbon)+biomass(:,ivm,iheartabove,icarbon))
847	ENDIF
848	ENDDO
849
850	model_total = zero
851	DO ipts = 1,npts
852	model_total = model_total + total_abwood_harvest(ipts)
853	ENDDO
854	WRITE(numout,*) 'model_total: ',model_total
855
856	CALL allreduce_sum(input_total, input_total_allproc)
857	WRITE(numout,*) 'Regional total input wood harvest,', input_total_allproc
858	IF ( input_total_allproc .GT. min_stomate) THEN
859	CALL allreduce_sum(model_total, model_total_allproc)
860	ratio = model_total_allproc/input_total_allproc
861	WRITE(numout,*) 'Initial ratio is: ',ratio
862	WRITE(numout,*) 'Initial regional total output wood harvest,', model_total_allproc
863	ENDIF
864
865
866	! We adjust the harvest_matrix if the output total harvet wood
867	! geos beyond 10% difference with the input harvest biomass.
868	num=1
869	DO WHILE ( (ABS(ratio-1) .GT. 0.01) .AND. (ratio.GT.min_stomate) .AND. num .LT. 20)
870
871	! Note here we treat all harvest as primary harvest
872	harvest_matrix(:,1) = harvest_matrix(:,1)/ratio
873	CALL MulAgeC_fh_firstday(npts,veget_max,newvegfrac,harvest_matrix, &
874	glcc_pft,glcc_pftmtc, &
875	Deficit_pf2yf_final, Deficit_sf2yf_final, &
876	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
877
878	! calculate model_total
879	total_abwood_harvest(:) = zero
880	DO ivm = 1,nvm
881	IF (is_tree(ivm) ) THEN
882	total_abwood_harvest(:) = total_abwood_harvest(:) + glcc_pft(:,ivm) area_land_m2(:) &
883	(biomass(:,ivm,isapabove,icarbon)+biomass(:,ivm,iheartabove,icarbon))
884	ENDIF
885	ENDDO
886
887	model_total = zero
888	DO ipts = 1,npts
889	model_total = model_total + total_abwood_harvest(ipts)
890	ENDDO
891
892	!! recalculating the ratio.
893	CALL allreduce_sum(input_total, input_total_allproc)
894	IF ( input_total_allproc .GT. min_stomate) THEN
895	CALL allreduce_sum(model_total, model_total_allproc)
896	ratio = model_total_allproc/input_total_allproc
897	WRITE(numout,*) 'Get_harvest_matrix:: ratio = ',ratio,'model_total_allproc,',model_total_allproc
898	ENDIF
899
900	num=num+1
901	WRITE(numout,*) 'Get_harvest_matrix:: Wood harvest matrix iteration: ',num
902
903	ENDDO
904
905	END SUBROUTINE Get_harvest_matrix
906	END MODULE stomate_fharvest_MulAgeC

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

Context Navigation

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

Download in other formats: