Context Navigation

stomate_fharvest_MulAgeC.f90 @ 6737

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

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

Context Navigation

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

Download in other formats: