Context Navigation

stomate_gluc_common.f90 @ 6737

Last change on this file since 6737 was 5258, checked in by chao.yue, 6 years ago
adjustment mass balance threshold
File size: 126.6 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_gluc_common
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 contains common fuctions and subroutines used by
10	! gross land use change and forestry harvest modules.
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S):
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_gluc_common
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
36	IMPLICIT NONE
37
38	PRIVATE
39	PUBLIC calc_cover, cross_give_receive, &
40	initialize_proxy_pft, sap_take, collect_legacy_pft, &
41	collect_legacy_pft_forestry, &
42	add_incoming_proxy_pft, empty_pft, build_age_index, &
43	prepare_balance_check, luc_balance_check
44
45
46	CONTAINS
47
48	!! ================================================================================================================================
49	!! SUBROUTINE : build_age_index
50	!!
51	!>\BRIEF
52	!!
53	!! DESCRIPTION : None
54	!!
55	!! RECENT CHANGE(S) : None
56	!!
57	!! MAIN OUTPUT VARIABLE(S): None
58	!!
59	!! REFERENCES : None
60	!!
61	!! FLOWCHART : None
62	!! \n
63	!_ ================================================================================================================================
64
65	SUBROUTINE build_age_index
66
67	INTEGER, ALLOCATABLE, SAVE :: indall_tree(:) !! Indices for all tree PFTs
68	INTEGER, ALLOCATABLE, SAVE :: indold_tree(:) !! Indices for old tree cohort only
69	INTEGER, ALLOCATABLE, SAVE :: indagec_tree(:,:) !! Indices for secondary tree cohorts,
70	!! note the sequence is old->young.
71	INTEGER, ALLOCATABLE, SAVE :: indall_grass(:) !! Indices for all grass PFTs
72	INTEGER, ALLOCATABLE, SAVE :: indold_grass(:) !! Indices for old grasses only
73	INTEGER, ALLOCATABLE, SAVE :: indagec_grass(:,:) !! Indices for secondary grass cohorts
74	!! note the sequence is old->young.
75	INTEGER, ALLOCATABLE, SAVE :: indall_pasture(:) !! Indices for all pasture PFTs
76	INTEGER, ALLOCATABLE, SAVE :: indold_pasture(:) !! Indices for old pasture only
77	INTEGER, ALLOCATABLE, SAVE :: indagec_pasture(:,:) !! Indices for secondary pasture cohorts
78	!! note the sequence is old->young.
79	INTEGER, ALLOCATABLE, SAVE :: indall_crop(:) !! Indices for all crop PFTs
80	INTEGER, ALLOCATABLE, SAVE :: indold_crop(:) !! Indices for old crops only
81	INTEGER, ALLOCATABLE, SAVE :: indagec_crop(:,:) !! Indices for secondary crop cohorts
82	!! note the sequence is old->young.
83	INTEGER :: num_tree_mulagec,num_grass_mulagec, &
84	num_pasture_mulagec,num_crop_mulagec, &
85	itree,itree2,igrass,igrass2,ipasture,ipasture2,icrop,icrop2
86	INTEGER :: i,j,ivma,staind,endind,ivm
87
88
89	!! 1. We first build all different indices that we are going to use
90	!! in handling the PFT exchanges, three types of indices are built:
91	!! - for all age classes
92	!! - include only oldest age classes
93	!! - include all age classes excpet the oldest ones
94	! We have to build these indices because we would like to extract from
95	! donating PFTs in the sequnce of old->young age classes or the revserse,
96	! and add in the receving PFTs only in the youngest-age-class PFTs. These
97	! indicies allow us to know where the different age classes are.
98
99	! calculate the total number of MTCs for each vegetation type.
100	num_tree_mulagec=0
101	num_grass_mulagec=0
102	num_pasture_mulagec=0
103	num_crop_mulagec=0
104
105	!! 1.1 Calculate the number of PFTs for different MTCs and allocate
106	!! the old and all indices arrays.
107
108	! [Note here the sequence to identify tree,pasture,grass,crop] is
109	! critical. The similar sequence is used in the subroutine "calc_cover".
110	! Do not forget to change the sequence there if you modify here.
111	DO ivma =2,nvmap
112	staind=start_index(ivma)
113	IF (nagec_pft(ivma)==1) THEN
114	WRITE(numout,*) "Error: metaclass has only a single age group: ",ivma
115	STOP
116	ELSE
117	IF (is_tree(staind)) THEN
118	num_tree_mulagec = num_tree_mulagec+1
119	ELSE IF (is_grassland_manag(staind)) THEN
120	num_pasture_mulagec = num_pasture_mulagec+1
121	ELSE IF (natural(staind)) THEN
122	num_grass_mulagec = num_grass_mulagec+1
123	ELSE
124	num_crop_mulagec = num_crop_mulagec+1
125	ENDIF
126	ENDIF
127	ENDDO
128
129	!! Allocate index array
130	! allocate all index
131	ALLOCATE(indall_tree(num_tree_mulagec*nagec_tree))
132	ALLOCATE(indall_grass(num_grass_mulagec*nagec_herb))
133	ALLOCATE(indall_pasture(num_pasture_mulagec*nagec_herb))
134	ALLOCATE(indall_crop(num_crop_mulagec*nagec_herb))
135
136	! allocate old-ageclass index
137	ALLOCATE(indold_tree(num_tree_mulagec))
138	ALLOCATE(indold_grass(num_grass_mulagec))
139	ALLOCATE(indold_pasture(num_pasture_mulagec))
140	ALLOCATE(indold_crop(num_crop_mulagec))
141
142	!! 1.2 Fill the oldest-age-class and all index arrays
143	itree=0
144	igrass=0
145	ipasture=0
146	icrop=0
147	itree2=1
148	igrass2=1
149	ipasture2=1
150	icrop2=1
151	DO ivma =2,nvmap
152	staind=start_index(ivma)
153	IF (is_tree(staind)) THEN
154	itree=itree+1
155	indold_tree(itree) = staind+nagec_pft(ivma)-1
156	DO j = 0,nagec_pft(ivma)-1
157	indall_tree(itree2+j) = staind+j
158	ENDDO
159	itree2=itree2+nagec_pft(ivma)
160	ELSE IF (natural(staind) .AND. .NOT. is_grassland_manag(staind)) THEN
161	igrass=igrass+1
162	indold_grass(igrass) = staind+nagec_pft(ivma)-1
163	DO j = 0,nagec_pft(ivma)-1
164	indall_grass(igrass2+j) = staind+j
165	ENDDO
166	igrass2=igrass2+nagec_pft(ivma)
167	ELSE IF (is_grassland_manag(staind)) THEN
168	ipasture = ipasture+1
169	indold_pasture(ipasture) = staind+nagec_pft(ivma)-1
170	DO j = 0,nagec_pft(ivma)-1
171	indall_pasture(ipasture2+j) = staind+j
172	ENDDO
173	ipasture2=ipasture2+nagec_pft(ivma)
174	ELSE
175	icrop = icrop+1
176	indold_crop(icrop) = staind+nagec_pft(ivma)-1
177	DO j = 0,nagec_pft(ivma)-1
178	indall_crop(icrop2+j) = staind+j
179	ENDDO
180	icrop2=icrop2+nagec_pft(ivma)
181	ENDIF
182	ENDDO
183
184	!! 1.3 Allocate and fill other age class index
185
186	ALLOCATE(indagec_tree(num_tree_mulagec,nagec_tree-1))
187	ALLOCATE(indagec_grass(num_grass_mulagec,nagec_herb-1))
188	ALLOCATE(indagec_pasture(num_pasture_mulagec,nagec_herb-1))
189	ALLOCATE(indagec_crop(num_crop_mulagec,nagec_herb-1))
190
191	! fill the non-oldest age class index arrays when number of age classes
192	! is more than 1.
193	itree=0
194	igrass=0
195	ipasture=0
196	icrop=0
197	DO ivma = 2,nvmap
198	staind=start_index(ivma)
199	IF (nagec_pft(ivma) > 1) THEN
200	IF (is_tree(staind)) THEN
201	itree=itree+1
202	DO j = 1,nagec_tree-1
203	indagec_tree(itree,j) = staind+nagec_tree-j-1
204	ENDDO
205	ELSE IF (natural(staind) .AND. .NOT. is_grassland_manag(staind)) THEN
206	igrass=igrass+1
207	DO j = 1,nagec_herb-1
208	indagec_grass(igrass,j) = staind+nagec_herb-j-1
209	ENDDO
210	ELSE IF (is_grassland_manag(staind)) THEN
211	ipasture=ipasture+1
212	DO j = 1,nagec_herb-1
213	indagec_pasture(ipasture,j) = staind+nagec_herb-j-1
214	ENDDO
215	ELSE
216	icrop=icrop+1
217	DO j = 1,nagec_herb-1
218	indagec_crop(icrop,j) = staind+nagec_herb-j-1
219	ENDDO
220	ENDIF
221	ENDIF
222	ENDDO
223
224	END SUBROUTINE build_age_index
225
226	! ================================================================================================================================
227	!! SUBROUTINE : calc_cover
228	!!
229	!>\BRIEF Calculate coverage fraction for different age classes of forest,
230	!! grass, pasture and crops and also for each metaclass. Note baresoil is excluded.
231	!!
232	!! DESCRIPTION :
233	!! Note:
234	!! 1. "calc_cover" subroutine does not depend on how many age classes
235	!! there are in each MTC.
236	!! 2. Fraction of baresoil is excluded here. This means transformation
237	!! of baresoil to a vegetated PFT is excluded in gross land cover change.
238	!!
239	!!
240	!! MAIN OUTPUT VARIABLE(S) :
241	!!
242	!! \n
243	!_ ================================================================================================================================
244	SUBROUTINE calc_cover(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
245	vegagec_pasture,vegagec_crop)
246
247
248	IMPLICIT NONE
249
250	!! Input variables
251	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
252	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
253
254	!! Output variables
255	REAL(r_std), DIMENSION(npts,nvmap), INTENT(inout) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
256	REAL(r_std), DIMENSION(npts,nagec_tree), INTENT(inout) :: vegagec_tree !! fraction of tree age-class groups, in sequence of old->young
257	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_grass !! fraction of grass age-class groups, in sequence of old->young
258	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_pasture !! fraction of pasture age-class groups, in sequence of old->young
259	REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout) :: vegagec_crop !! fraction of crop age-class groups, in sequence of old->young
260
261	!! Local variables
262	INTEGER(i_std) :: ivma,staind,endind,j !! indices (unitless)
263
264	veget_mtc(:,:) = 0.
265	vegagec_tree(:,:) = 0.
266	vegagec_grass(:,:) = 0.
267	vegagec_pasture(:,:) = 0.
268	vegagec_crop(:,:) = 0.
269
270	! Calculate veget_max for MTCs
271	DO ivma = 1,nvmap
272	staind = start_index(ivma)
273	IF (nagec_pft(ivma) == 1) THEN
274	veget_mtc(:,ivma) = veget_max(:,staind)
275	ELSE
276	veget_mtc(:,ivma) = \
277	SUM(veget_max(:,staind:staind+nagec_pft(ivma)-1),DIM=2)
278	ENDIF
279	ENDDO
280
281	! Calculate veget_max for each age class
282	DO ivma = 2,nvmap !here we start with 2 to exclude baresoil (always PFT1)
283	staind = start_index(ivma)
284	endind = staind+nagec_pft(ivma)-1
285
286	! Single-age-class MTC goest to oldest age class.
287	IF (nagec_pft(ivma) == 1) THEN
288	IF (is_tree(staind)) THEN
289	vegagec_tree(:,1) = vegagec_tree(:,1)+veget_max(:,staind)
290	ELSE IF (is_grassland_manag(staind)) THEN
291	vegagec_pasture(:,1) = vegagec_pasture(:,1)+veget_max(:,staind)
292	ELSE IF (natural(staind)) THEN
293	vegagec_grass(:,1) = vegagec_grass(:,1)+veget_max(:,staind)
294	ELSE
295	vegagec_crop(:,1) = vegagec_crop(:,1)+veget_max(:,staind)
296	ENDIF
297
298	ELSE
299	IF (is_tree(staind)) THEN
300	DO j=1,nagec_tree
301	vegagec_tree(:,j) = vegagec_tree(:,j)+veget_max(:,endind-j+1)
302	ENDDO
303	ELSE IF (is_grassland_manag(staind)) THEN
304	DO j=1,nagec_herb
305	vegagec_pasture(:,j) = vegagec_pasture(:,j)+veget_max(:,endind-j+1)
306	ENDDO
307	ELSE IF (natural(staind)) THEN
308	DO j=1,nagec_herb
309	vegagec_grass(:,j) = vegagec_grass(:,j)+veget_max(:,endind-j+1)
310	ENDDO
311	ELSE
312	DO j=1,nagec_herb
313	vegagec_crop(:,j) = vegagec_crop(:,j)+veget_max(:,endind-j+1)
314	ENDDO
315	ENDIF
316	ENDIF
317	ENDDO
318
319	END SUBROUTINE calc_cover
320
321
322	! ================================================================================================================================
323	!! SUBROUTINE : cross_give_receive
324	!!
325	!>\BRIEF : Allocate the outgoing and receving fractions in respective
326	!! PFTs.
327	!! \n
328	!! Notes:
329	!! 1. veget_max is subtracted when fractions are taken out, but newly added
330	!! fractions in the youngest age class is not added, to avoid this newly
331	!! created fractions being used again the following transitions. This is
332	!! is reasonable because the newly created youngest-age-class PFT fractions
333	!! have nothing but small sapling biomass and it's unreasonable to use it
334	!! for any further land use conversion activities.
335	!_ ================================================================================================================================
336	SUBROUTINE cross_give_receive(ipts,frac_used,veget_mtc,newvegfrac, &
337	indold_tree,indagec_crop,nagec_receive,num_crop_mulagec, &
338	veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp)
339
340
341	IMPLICIT NONE
342
343	!! 0. Input variables
344	INTEGER, INTENT(in) :: ipts
345	REAL(r_std), INTENT(in) :: frac_used !! fraction that the giving PFTs are going to collectively give
346	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_mtc !! "maximal" coverage fraction of a PFT on the ground
347	INTEGER, DIMENSION(:), INTENT(in) :: indold_tree !! Indices for PFTs giving out fractions;
348	!! here use old tree cohort as an example
349	INTEGER, DIMENSION(:,:), INTENT(in) :: indagec_crop !! Indices for secondary basic-vegetation cohorts; The youngest age classes
350	!! of these vegetations are going to receive fractions.
351	!! here we use crop cohorts as an example
352	INTEGER, INTENT(in) :: num_crop_mulagec !! number of crop MTCs with more than one age classes
353	INTEGER, INTENT(in) :: nagec_receive !! number of age classes in the receiving basic types
354	!! (i.e., tree, grass, pasture, crop), here we can use crop
355	!! as an example, nagec_receive=nagec_herb
356	REAL(r_std), DIMENSION(:,:),INTENT(in) :: newvegfrac !!
357
358	!! 1. Modified variables
359	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
360	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fractions each PFT is going to lose
361	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: glcc_pftmtc !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
362	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
363	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: glcc_pft_tmp !! a temporary variable to hold the fractions each PFT is going to lose
364
365	!! Local vriables
366	INTEGER :: j,ipft, iyoung
367	REAL(r_std) :: totalveg, tveg_now
368	LOGICAL :: guide_newfrac
369
370	guide_newfrac = gluc_newfrac_guide
371
372	! Out final objective is to know glcc_pftmtc, i.e., the fraction from each PFT
373	! to the youngest age group of each MTC. We separate this task into two steps:
374	! 1. we allocate the total outgoing fraction into the same age-class PFTs of
375	! the a basic-vegetation (for example, the same age-calss PFTs of forest);
376	! 2. we further allocate the outgoing fraction of each age-class PFT to
377	! the different receiving youngest age-class PFTs of the same basic-vegetation
378	! type, for example, the youngest age-calss PFTs of cropland.
379
380	! glcc_pft_tmp used only as a temporary variable to store the value
381	glcc_pft_tmp(ipts,indold_tree) = veget_max(ipts,indold_tree)/SUM(veget_max(ipts,indold_tree))*frac_used
382	glcc_pft(ipts,indold_tree) = glcc_pft(ipts,indold_tree) + glcc_pft_tmp(ipts,indold_tree)
383	!we have to remove the outgoing fraction from veget_max in order to use this information for next loop
384	veget_max(ipts,indold_tree) = veget_max(ipts,indold_tree) - glcc_pft_tmp(ipts,indold_tree)
385
386	! when receiving basic-vegetation type has a single age group, it will be considered as
387	! both old and young age group (thus recevie the fraction donation), otherwise the youngest
388	! age group is always the final element of indagec_crop.
389	IF (nagec_receive == 1) THEN
390	iyoung = 1
391	ELSE
392	iyoung = nagec_receive - 1
393	ENDIF
394
395	! [20160728] Here we have two options:
396	! 1. allocate the newly created young age class according to existing fractions
397	! the MTCs.
398	! 2. Use the fractions of MTCs from the current-day PFT map to guid the
399	! allocation of newly created young-age-class MTCs.
400
401	totalveg = 0. ! [20160130 note here totalveg is the total fraction
402	! of all existing MTCs that are going to recieve newly
403	! convervted fractions.]
404	tveg_now = 0. ! total vegetation fraction in the current-day MTC
405	! input map.
406	DO j=1,num_crop_mulagec
407	totalveg = totalveg + veget_mtc(ipts,agec_group(indagec_crop(j,iyoung)))
408	tveg_now = tveg_now + newvegfrac(ipts,agec_group(indagec_crop(j,iyoung)))
409	ENDDO
410
411	IF (guide_newfrac) THEN
412	IF (tveg_now>min_stomate) THEN
413	DO j=1,num_crop_mulagec
414	ipft = indagec_crop(j,iyoung)
415	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) = &
416	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) + &
417	glcc_pft_tmp(ipts,indold_tree) * newvegfrac(ipts,agec_group(ipft))/tveg_now
418	ENDDO
419	ELSE
420	DO j=1,num_crop_mulagec
421	ipft = indagec_crop(j,iyoung)
422	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) = &
423	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) + glcc_pft_tmp(ipts,indold_tree)/num_crop_mulagec
424	ENDDO
425	ENDIF
426
427	ELSE
428	IF (totalveg>min_stomate) THEN
429	DO j=1,num_crop_mulagec
430	ipft = indagec_crop(j,iyoung)
431	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) = &
432	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) + &
433	glcc_pft_tmp(ipts,indold_tree) * veget_mtc(ipts,agec_group(ipft))/totalveg
434	ENDDO
435	ELSE
436	DO j=1,num_crop_mulagec
437	ipft = indagec_crop(j,iyoung)
438	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) = &
439	glcc_pftmtc(ipts,indold_tree,agec_group(ipft)) + glcc_pft_tmp(ipts,indold_tree)/num_crop_mulagec
440	ENDDO
441	ENDIF
442
443	ENDIF
444
445	END SUBROUTINE cross_give_receive
446
447	! ================================================================================================================================
448	!! SUBROUTINE : clear_forest
449	!!
450	!>\BRIEF : Handle forest harvest before its legacy is transferred to
451	! newly initialized youngest-age-class PFT.
452	!!
453	!>\DESCRIPTION
454	!_ ================================================================================================================================
455	!!++TEMP++ biomass,veget_frac are not used because the remaining biomass to be
456	!! harvested is calculated within the deforestation fire module.
457	SUBROUTINE clear_forest (npts,ipts,ivm,biomass,frac, &
458	instant_loss, &
459	litter, deforest_biomass_remain,&
460	fuel_1hr,fuel_10hr,&
461	fuel_100hr,fuel_1000hr,&
462	lignin_struc,&
463	bm_to_litter_pro,convflux,prod10,prod100,&
464	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
465	fuel_1000hr_pro, lignin_content_pro)
466
467
468	IMPLICIT NONE
469
470	!! 0.1 Input variables
471	INTEGER, INTENT(in) :: npts
472	INTEGER, INTENT(in) :: ipts
473	INTEGER, INTENT(in) :: ivm
474	REAL(r_std), INTENT(in) :: instant_loss
475	REAL(r_std), INTENT(in) :: frac !! the fraction of land covered by forest to be deforested
476	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
477	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
478	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
479	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
480	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
481	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! Vegetmax-weighted remaining litter on the ground for
482	!! deforestation region.
483	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
484	!! deforestation region.
485	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
486	!! above and below ground
487
488	!! 0.2 Modified variables
489	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro !! conversion of biomass to litter
490	!! @tex ($gC m^{-2} day^{-1}$) @endtex
491	REAL(r_std), DIMENSION(:), INTENT(inout) :: convflux !! release during first year following land cover
492	!! change
493
494	REAL(r_std), DIMENSION(npts,0:10), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
495	!! pool after the annual release for each
496	!! compartment (10 + 1 : input from year of land
497	!! cover change)
498	REAL(r_std), DIMENSION(npts,0:100), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
499	!! pool after the annual release for each
500	!! compartment (100 + 1 : input from year of land
501	!! cover change)
502
503	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: litter_pro
504	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1hr_pro
505	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_10hr_pro
506	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_100hr_pro
507	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1000hr_pro
508	REAL(r_std), DIMENSION(:),INTENT(inout) :: lignin_content_pro
509
510
511
512	!! 0.4 Local variables
513	REAL(r_std) :: above
514
515	! harvest of aboveground sap- and heartwood biomass after taking into
516	! account of deforestation fire
517	IF (allow_deforest_fire) THEN
518	above = deforest_biomass_remain(ipts,ivm,isapabove,icarbon)+ &
519	deforest_biomass_remain(ipts,ivm,iheartabove,icarbon)
520	convflux(ipts) = convflux(ipts) + 0
521	prod10(ipts,0) = prod10(ipts,0) + 0.4*above
522	prod100(ipts,0) = prod100(ipts,0) + 0.6*above
523	ELSE
524	above = (biomass(ipts,ivm,isapabove,icarbon)+ &
525	biomass(ipts,ivm,ileaf,icarbon) + &
526	biomass(ipts,ivm,ifruit,icarbon) + &
527	biomass(ipts,ivm,iheartabove,icarbon))*frac
528	convflux(ipts) = convflux(ipts) + instant_loss * above
529	!prod10(ipts,0) = prod10(ipts,0) + coeff_lcchange_10(ivm) * above
530	!prod100(ipts,0) = prod100(ipts,0) + coeff_lcchange_100(ivm) * above
531	ENDIF
532
533	! the transfer of dead biomass to litter
534
535	bm_to_litter_pro(isapabove,:) = bm_to_litter_pro(isapabove,:) + &
536	biomass(ipts,ivm,isapabove,:)frac(1-instant_loss)
537	bm_to_litter_pro(iheartabove,:) = bm_to_litter_pro(iheartabove,:) + &
538	biomass(ipts,ivm,iheartabove,:)frac(1-instant_loss)
539
540	bm_to_litter_pro(isapbelow,:) = bm_to_litter_pro(isapbelow,:) + &
541	biomass(ipts,ivm,isapbelow,:)*frac
542	bm_to_litter_pro(iheartbelow,:) = bm_to_litter_pro(iheartbelow,:) + &
543	biomass(ipts,ivm,iheartbelow,:)*frac
544	bm_to_litter_pro(iroot,:) = bm_to_litter_pro(iroot,:) + &
545	biomass(ipts,ivm,iroot,:)*frac
546	bm_to_litter_pro(ifruit,:) = bm_to_litter_pro(ifruit,:) + &
547	biomass(ipts,ivm,ifruit,:)frac(1-instant_loss)
548	bm_to_litter_pro(icarbres,:) = bm_to_litter_pro(icarbres,:) + &
549	biomass(ipts,ivm,icarbres,:)*frac
550	bm_to_litter_pro(ileaf,:) = bm_to_litter_pro(ileaf,:) + &
551	biomass(ipts,ivm,ileaf,:)frac(1-instant_loss)
552
553	!update litter_pro
554	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
555	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
556	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
557	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
558	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
559	!don't forget to hanle litter lignin content
560	lignin_content_pro(:)= lignin_content_pro(:) + &
561	litter(ipts,istructural,ivm,:,icarbon)fraclignin_struc(ipts,ivm,:)
562
563	END SUBROUTINE clear_forest
564
565	! ================================================================================================================================
566	!! SUBROUTINE : harvest_industrial
567	!!
568	!>\BRIEF : Handle forest harvest before its legacy is transferred to
569	! newly initialized youngest-age-class PFT.
570	!!
571	!>\DESCRIPTION
572	!_ ================================================================================================================================
573	!!++TEMP++ biomass,veget_frac are not used because the remaining biomass to be
574	!! harvested is calculated within the deforestation fire module.
575	SUBROUTINE harvest_industrial (npts,ipts,ivm,biomass,frac, &
576	litter, deforest_biomass_remain,&
577	fuel_1hr,fuel_10hr,&
578	fuel_100hr,fuel_1000hr,&
579	lignin_struc,&
580	bm_to_litter_pro,convflux,prod10,prod100,&
581	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
582	fuel_1000hr_pro, lignin_content_pro)
583
584
585	IMPLICIT NONE
586
587	!! 0.1 Input variables
588	INTEGER, INTENT(in) :: npts
589	INTEGER, INTENT(in) :: ipts
590	INTEGER, INTENT(in) :: ivm
591	REAL(r_std), INTENT(in) :: frac !! the fraction of land covered by forest to be deforested
592	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
593	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
594	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
595	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
596	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
597	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! Vegetmax-weighted remaining litter on the ground for
598	!! deforestation region.
599	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
600	!! deforestation region.
601	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
602	!! above and below ground
603
604	!! 0.2 Modified variables
605	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro !! conversion of biomass to litter
606	!! @tex ($gC m^{-2} day^{-1}$) @endtex
607	REAL(r_std), DIMENSION(:), INTENT(inout) :: convflux !! release during first year following land cover
608	!! change
609
610	REAL(r_std), DIMENSION(npts,0:10), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
611	!! pool after the annual release for each
612	!! compartment (10 + 1 : input from year of land
613	!! cover change)
614	REAL(r_std), DIMENSION(npts,0:100), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
615	!! pool after the annual release for each
616	!! compartment (100 + 1 : input from year of land
617	!! cover change)
618
619	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: litter_pro
620	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1hr_pro
621	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_10hr_pro
622	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_100hr_pro
623	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1000hr_pro
624	REAL(r_std), DIMENSION(:),INTENT(inout) :: lignin_content_pro
625
626
627
628	!! 0.4 Local variables
629	REAL(r_std) :: above
630	REAL(r_std) :: slash_frac
631	REAL(r_std) :: harvest_efficiency
632
633	! harvest_efficiency indicates the percentage of the gross area that
634	! are effetively harvested to reach the target.
635	harvest_efficiency = 0.667
636	! harvest of aboveground sap- and heartwood biomass after taking into
637	! account of deforestation fire
638	IF (allow_deforest_fire) THEN
639	above = deforest_biomass_remain(ipts,ivm,isapabove,icarbon)+ &
640	deforest_biomass_remain(ipts,ivm,iheartabove,icarbon)
641	convflux(ipts) = convflux(ipts) + 0
642	prod10(ipts,0) = prod10(ipts,0) + 0.4*above
643	prod100(ipts,0) = prod100(ipts,0) + 0.6*above
644	ELSE
645	above = (biomass(ipts,ivm,isapabove,icarbon)+ &
646	biomass(ipts,ivm,iheartabove,icarbon))*frac
647	convflux(ipts) = convflux(ipts) + harvest_efficiency * coeff_indwood_1(ivm) * above
648	prod10(ipts,0) = prod10(ipts,0) + harvest_efficiency * coeff_indwood_10(ivm) * above
649	prod100(ipts,0) = prod100(ipts,0) + harvest_efficiency * coeff_indwood_100(ivm) * above
650	ENDIF
651
652	slash_frac = 1 - (coeff_indwood_10(ivm) + coeff_indwood_100(ivm) &
653	+ coeff_indwood_1(ivm)) * harvest_efficiency
654	! the transfer of dead biomass to litter
655	bm_to_litter_pro(isapabove,:) = bm_to_litter_pro(isapabove,:) + &
656	biomass(ipts,ivm,isapabove,:)fracslash_frac
657	bm_to_litter_pro(iheartabove,:) = bm_to_litter_pro(iheartabove,:) + &
658	biomass(ipts,ivm,iheartabove,:)fracslash_frac
659	bm_to_litter_pro(isapbelow,:) = bm_to_litter_pro(isapbelow,:) + &
660	biomass(ipts,ivm,isapbelow,:)*frac
661	bm_to_litter_pro(iheartbelow,:) = bm_to_litter_pro(iheartbelow,:) + &
662	biomass(ipts,ivm,iheartbelow,:)*frac
663	bm_to_litter_pro(iroot,:) = bm_to_litter_pro(iroot,:) + &
664	biomass(ipts,ivm,iroot,:)*frac
665	bm_to_litter_pro(ifruit,:) = bm_to_litter_pro(ifruit,:) + &
666	biomass(ipts,ivm,ifruit,:)*frac
667	bm_to_litter_pro(icarbres,:) = bm_to_litter_pro(icarbres,:) + &
668	biomass(ipts,ivm,icarbres,:)*frac
669	bm_to_litter_pro(ileaf,:) = bm_to_litter_pro(ileaf,:) + &
670	biomass(ipts,ivm,ileaf,:)*frac
671
672	!update litter_pro
673	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
674	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
675	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
676	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
677	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
678	!don't forget to hanle litter lignin content
679	lignin_content_pro(:)= lignin_content_pro(:) + &
680	litter(ipts,istructural,ivm,:,icarbon)fraclignin_struc(ipts,ivm,:)
681
682	END SUBROUTINE harvest_industrial
683
684	! ================================================================================================================================
685	!! SUBROUTINE : harvest_fuelwood
686	!!
687	!>\BRIEF : Handle forest harvest before its legacy is transferred to
688	! newly initialized youngest-age-class PFT.
689	!!
690	!>\DESCRIPTION
691	!_ ================================================================================================================================
692	!!++TEMP++ biomass,veget_frac are not used because the remaining biomass to be
693	!! harvested is calculated within the deforestation fire module.
694	SUBROUTINE harvest_fuelwood (npts,ipts,ivm,biomass,frac, &
695	litter, deforest_biomass_remain,&
696	fuel_1hr,fuel_10hr,&
697	fuel_100hr,fuel_1000hr,&
698	lignin_struc,&
699	bm_to_litter_pro,convflux,prod10,prod100,&
700	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
701	fuel_1000hr_pro, lignin_content_pro)
702
703
704	IMPLICIT NONE
705
706	!! 0.1 Input variables
707	INTEGER, INTENT(in) :: npts
708	INTEGER, INTENT(in) :: ipts
709	INTEGER, INTENT(in) :: ivm
710	REAL(r_std), INTENT(in) :: frac !! the fraction of land covered by forest to be deforested
711	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
712	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
713	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
714	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
715	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
716	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! Vegetmax-weighted remaining litter on the ground for
717	!! deforestation region.
718	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
719	!! deforestation region.
720	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
721	!! above and below ground
722
723	!! 0.2 Modified variables
724	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro !! conversion of biomass to litter
725	!! @tex ($gC m^{-2} day^{-1}$) @endtex
726	REAL(r_std), DIMENSION(:), INTENT(inout) :: convflux !! release during first year following land cover
727	!! change
728
729	REAL(r_std), DIMENSION(npts,0:10), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
730	!! pool after the annual release for each
731	!! compartment (10 + 1 : input from year of land
732	!! cover change)
733	REAL(r_std), DIMENSION(npts,0:100), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
734	!! pool after the annual release for each
735	!! compartment (100 + 1 : input from year of land
736	!! cover change)
737
738	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: litter_pro
739	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1hr_pro
740	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_10hr_pro
741	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_100hr_pro
742	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1000hr_pro
743	REAL(r_std), DIMENSION(:),INTENT(inout) :: lignin_content_pro
744
745
746
747	!! 0.4 Local variables
748	REAL(r_std) :: above
749
750	! harvest of aboveground sap- and heartwood biomass after taking into
751	! account of deforestation fire
752	IF (allow_deforest_fire) THEN
753	above = deforest_biomass_remain(ipts,ivm,isapabove,icarbon)+ &
754	deforest_biomass_remain(ipts,ivm,iheartabove,icarbon)
755	convflux(ipts) = convflux(ipts) + 0
756	prod10(ipts,0) = prod10(ipts,0) + 0.4*above
757	prod100(ipts,0) = prod100(ipts,0) + 0.6*above
758	ELSE
759	above = (biomass(ipts,ivm,isapabove,icarbon)+ &
760	biomass(ipts,ivm,iheartabove,icarbon))*frac
761	! we assume no wood goes to 10-year and 100-year pool in fuelwood collection.
762	convflux(ipts) = convflux(ipts) + above
763	ENDIF
764
765	! the transfer of dead biomass to litter
766	bm_to_litter_pro(isapbelow,:) = bm_to_litter_pro(isapbelow,:) + &
767	biomass(ipts,ivm,isapbelow,:)*frac
768	bm_to_litter_pro(iheartbelow,:) = bm_to_litter_pro(iheartbelow,:) + &
769	biomass(ipts,ivm,iheartbelow,:)*frac
770	bm_to_litter_pro(iroot,:) = bm_to_litter_pro(iroot,:) + &
771	biomass(ipts,ivm,iroot,:)*frac
772	bm_to_litter_pro(ifruit,:) = bm_to_litter_pro(ifruit,:) + &
773	biomass(ipts,ivm,ifruit,:)*frac
774	bm_to_litter_pro(icarbres,:) = bm_to_litter_pro(icarbres,:) + &
775	biomass(ipts,ivm,icarbres,:)*frac
776	bm_to_litter_pro(ileaf,:) = bm_to_litter_pro(ileaf,:) + &
777	biomass(ipts,ivm,ileaf,:)*frac
778
779	!update litter_pro
780	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
781	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
782	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
783	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
784	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
785	!don't forget to hanle litter lignin content
786	lignin_content_pro(:)= lignin_content_pro(:) + &
787	litter(ipts,istructural,ivm,:,icarbon)fraclignin_struc(ipts,ivm,:)
788
789	END SUBROUTINE harvest_fuelwood
790
791	! ================================================================================================================================
792	!! SUBROUTINE : harvest_herb
793	!!
794	!>\BRIEF : Handle herbaceous PFT clearing before its legacy is transferred to
795	! newly initialized youngest-age-class PFT.
796	!!
797	!>\DESCRIPTION
798	!_ ================================================================================================================================
799	SUBROUTINE harvest_herb (ipts,ivm,biomass,veget_frac,bm_to_litter_pro)
800
801	IMPLICIT NONE
802
803	!! 0.1 Input variables
804	INTEGER, INTENT(in) :: ipts
805	INTEGER, INTENT(in) :: ivm
806	REAL(r_std), INTENT(in) :: veget_frac !! the fraction of land covered by herbaceous PFT to be cleared
807	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
808
809	!! 0.2 Modified variables
810	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro
811
812
813	INTEGER :: ipart, num
814
815	! the transfer of dead biomass to litter
816	DO ipart = 1,nparts
817	bm_to_litter_pro(ipart,:) = bm_to_litter_pro(ipart,:) + biomass(ipts,ivm,ipart,:)*veget_frac
818	ENDDO
819
820	END SUBROUTINE harvest_herb
821
822
823	! ================================================================================================================================
824	!! SUBROUTINE : initialize_proxy_pft
825	!!
826	!>\BRIEF Initialize a proxy new youngest age class PFT.
827	!!
828	!>\DESCRIPTION Initialize a proxy new youngest age class PFT that will be
829	!! merged with existing yongest age class, or fill the empty
830	!! niche of the youngest age class PFT.
831	!_ ================================================================================================================================
832	SUBROUTINE initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro, &
833	biomass_pro, co2_to_bm_pro, ind_pro, age_pro, &
834	senescence_pro, PFTpresent_pro, &
835	lm_lastyearmax_pro, everywhere_pro, npp_longterm_pro, &
836	leaf_frac_pro,leaf_age_pro)
837
838	IMPLICIT NONE
839
840	!! 0.1 Input variables
841	INTEGER, INTENT(in) :: ipts !!
842	INTEGER, INTENT(in) :: ipft_young_agec !! index of the concerned youngest-age-class PFT
843	REAL(r_std), INTENT(in) :: veget_max_pro !! fraction of grid cell land area that's to be occupied
844
845	!! 0.2 Modified variables
846	REAL(r_std), INTENT(inout) :: co2_to_bm_pro
847
848	!! 0.3 Output variables
849	REAL(r_std), DIMENSION(:,:), INTENT(out) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
850	REAL(r_std), DIMENSION(:), INTENT(out) :: leaf_frac_pro !! fraction of leaves in leaf age class
851	REAL(r_std), DIMENSION(:), INTENT(out) :: leaf_age_pro !! fraction of leaves in leaf age class
852	REAL(r_std), INTENT(out) :: age_pro, ind_pro, lm_lastyearmax_pro
853	REAL(r_std), INTENT(out) :: npp_longterm_pro
854	REAL(r_std), INTENT(out) :: everywhere_pro !! is the PFT everywhere in the grid box or very
855	LOGICAL, INTENT(out) :: senescence_pro !! plant senescent (only for deciduous trees) Set
856	!! to .FALSE. if PFT is introduced or killed
857	LOGICAL, INTENT(out) :: PFTpresent_pro !! Is pft there (unitless)
858
859	!! 0.4 Local variables
860	!REAL(r_std), DIMENSION(npts,nvm) :: when_growthinit !! how many days ago was the beginning of the
861	! !! growing season (days)
862
863	REAL(r_std), DIMENSION(nparts,nelements) :: bm_new !! biomass increase @tex ($gC m^{-2}$) @endtex
864	REAL(r_std) :: cn_ind,ind
865	INTEGER :: i,j,k,l
866
867	! -Note-
868	! This part of codes are copied from the original lcchange_main subroutine
869	! that initialize a new PFT.
870
871	i=ipts
872	j=ipft_young_agec
873
874	!! Initialization of some variables
875	leaf_frac_pro(:) = zero
876	leaf_age_pro(:) = zero
877
878	!! Initial setting of new establishment
879	IF (is_tree(j)) THEN
880	! cn_sapl(j)=0.5; stomate_data.f90
881	cn_ind = cn_sapl(j)
882	ELSE
883	cn_ind = un
884	ENDIF
885	ind = veget_max_pro / cn_ind
886	ind_pro = ind*veget_max_pro
887	PFTpresent_pro = .TRUE.
888	senescence_pro = .FALSE.
889	everywhere_pro = 1.*veget_max_pro
890	age_pro = zero
891
892	! large_value = 1.E33_r_std
893	! when_growthinit(i,j) = large_value
894	leaf_frac_pro(1) = 1.0 * veget_max_pro
895	leaf_age_pro(1) = 1.0 * veget_max_pro !This was not included in original lcchange_main subroutine
896	npp_longterm_pro = npp_longterm_init * veget_max_pro
897	lm_lastyearmax_pro = bm_sapl(j,ileaf,icarbon) * ind * veget_max_pro
898
899	!! Update of biomass in each each carbon stock component (leaf, sapabove, sapbelow,
900	!> heartabove, heartbelow, root, fruit, and carbres)\n
901	DO k = 1, nparts ! loop over # carbon stock components, nparts = 8; stomate_constant.f90
902	DO l = 1,nelements ! loop over # elements
903	biomass_pro(k,l) = ind * bm_sapl(j,k,l)
904	END DO ! loop over # elements
905	co2_to_bm_pro = co2_to_bm_pro + ind * bm_sapl(j,k,icarbon)
906	ENDDO ! loop over # carbon stock components
907
908	END SUBROUTINE initialize_proxy_pft
909
910	! ================================================================================================================================
911	!! SUBROUTINE sap_take
912	!!
913	!>\BRIEF : Take the sapling biomass of the new PFTs from the existing biomass, otherwise
914	! take from co2_to_bm
915	!!
916	!>\DESCRIPTION
917	!_ ================================================================================================================================
918	SUBROUTINE sap_take (ipts,ivma,veget_max,biomass_pro,biomass,co2_to_bm_pro)
919
920	INTEGER, INTENT(in) :: ipts !!
921	INTEGER, INTENT(in) :: ivma
922	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT (LAI ->
923	REAL(r_std), DIMENSION(:,:), INTENT(in) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
924
925	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
926	REAL(r_std), INTENT(inout) :: co2_to_bm_pro
927
928
929	REAL(r_std), DIMENSION(nparts,nelements) :: biomass_total !! biomass @tex ($gC m^{-2}$) @endtex
930	REAL(r_std) :: bm_org,bmpro_share
931	INTEGER :: i,ivm,ipart
932
933	biomass_total(:,:) = zero
934	bm_org = zero
935	bmpro_share = zero
936
937	DO i = 1,nagec_pft(ivma)
938	ivm = start_index(ivma)+i-1
939	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
940	biomass_total = biomass_total + biomass(ipts,ivm,:,:)*veget_max(ipts,ivm)
941	ENDIF
942	ENDDO
943
944	DO ipart = 1, nparts
945	IF (biomass_total(ipart,icarbon) .GT. biomass_pro(ipart,icarbon)) THEN
946	co2_to_bm_pro = co2_to_bm_pro - biomass_pro(ipart,icarbon)
947	!treat each PFT of the MTC
948	DO i = 1,nagec_pft(ivma)
949	ivm = start_index(ivma)+i-1
950	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
951	bm_org = biomass(ipts,ivm,ipart,icarbon) * veget_max(ipts,ivm)
952	bmpro_share = bm_org/biomass_total(ipart,icarbon) * biomass_pro(ipart,icarbon)
953	biomass(ipts,ivm,ipart,icarbon) = (bm_org - bmpro_share)/veget_max(ipts,ivm)
954	ENDIF
955	ENDDO
956	ENDIF
957	ENDDO
958
959	END SUBROUTINE sap_take
960
961	! ================================================================================================================================
962	!! SUBROUTINE collect_legacy_pft
963	!!
964	!>\BRIEF : Collect the legacy variables that are going to be included
965	! in the newly initialized PFT.
966	!!
967	!>\DESCRIPTION
968	!_ ================================================================================================================================
969	SUBROUTINE collect_legacy_pft(npts, ipts, ivma, glcc_pftmtc, &
970	biomass, bm_to_litter, carbon, litter, &
971	deepC_a, deepC_s, deepC_p, &
972	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
973	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
974	resp_maint, resp_growth, resp_hetero, co2_fire, &
975	def_fuel_1hr_remain, def_fuel_10hr_remain, &
976	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
977	deforest_litter_remain, deforest_biomass_remain, &
978	veget_max_pro, carbon_pro, lignin_struc_pro, litter_pro, &
979	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
980	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
981	bm_to_litter_pro, co2_to_bm_pro, gpp_daily_pro, &
982	npp_daily_pro, resp_maint_pro, resp_growth_pro, &
983	resp_hetero_pro, co2_fire_pro, &
984	convflux,prod10,prod100)
985
986	IMPLICIT NONE
987
988	!! 0.1 Input variables
989	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
990	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
991	INTEGER, INTENT(in) :: ivma !! Index for metaclass
992	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
993	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
994	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: bm_to_litter !! Transfer of biomass to litter
995	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
996	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: carbon !! carbon pool: active, slow, or passive
997	!! @tex ($gC m^{-2}$) @endtex
998	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_a !! Permafrost soil carbon (g/m**3) active
999	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1000	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1001	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! metabolic and structural litter, above and
1002	!! below ground @tex ($gC m^{-2}$) @endtex
1003	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
1004	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
1005	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
1006	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
1007	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1008	!! above and below ground
1009	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_to_bm !! biomass uptaken
1010	!! @tex ($gC m^{-2} day^{-1}$) @endtex
1011	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gpp_daily !! Daily gross primary productivity
1012	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1013	REAL(r_std), DIMENSION(:,:), INTENT(in) :: npp_daily !! Net primary productivity
1014	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1015	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_maint !! Maintenance respiration
1016	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1017	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_growth !! Growth respiration
1018	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1019	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_hetero !! Heterotrophic respiration
1020	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1021	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_fire !! Heterotrophic respiration
1022	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1023	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_1hr_remain
1024	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_10hr_remain
1025	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_100hr_remain
1026	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_1000hr_remain
1027	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
1028	!! deforestation region.
1029	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
1030	!! deforestation region.
1031
1032	!! 0.2 Output variables
1033	REAL(r_std), DIMENSION(:), INTENT(inout) :: carbon_pro
1034	REAL(r_std), DIMENSION(:), INTENT(inout) :: deepC_a_pro
1035	REAL(r_std), DIMENSION(:), INTENT(inout) :: deepC_s_pro
1036	REAL(r_std), DIMENSION(:), INTENT(inout) :: deepC_p_pro
1037	REAL(r_std), DIMENSION(:), INTENT(inout) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
1038	!! above and below ground
1039	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: litter_pro
1040	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1hr_pro
1041	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_10hr_pro
1042	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_100hr_pro
1043	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1000hr_pro
1044	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro
1045	REAL(r_std), INTENT(inout) :: veget_max_pro, co2_to_bm_pro
1046	REAL(r_std), INTENT(inout) :: gpp_daily_pro, npp_daily_pro
1047	REAL(r_std), INTENT(inout) :: resp_maint_pro, resp_growth_pro
1048	REAL(r_std), INTENT(inout) :: resp_hetero_pro, co2_fire_pro
1049
1050	!! 0.3 Modified variables
1051	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: convflux !! release during first year following land cover
1052	!! change
1053
1054	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
1055	!! pool after the annual release for each
1056	!! compartment (10 + 1 : input from year of land
1057	!! cover change)
1058	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
1059	!! pool after the annual release for each
1060	!! compartment (100 + 1 : input from year of land
1061	!! cover change)
1062
1063	!! 0.4 Local variables
1064	REAL(r_std), DIMENSION(nlevs) :: lignin_content_pro
1065	REAL(r_std) :: frac, instant_loss
1066	INTEGER :: ivm
1067
1068
1069	! All *_pro variables collect the legacy pools/fluxes of the ancestor
1070	! PFTs for the receiving youngest age class. All *_pro variables
1071	! represent the quantity weighted by the fraction of ancestor contributing
1072	! PFTs.
1073	! Exceptions:
1074	! lignin_struc_pro:: the ratio of lignin content in structural litter.
1075
1076	lignin_content_pro(:)=zero
1077
1078	DO ivm = 1,nvm
1079	frac = glcc_pftmtc(ipts,ivm,ivma)
1080	IF (frac>zero) THEN
1081	veget_max_pro = veget_max_pro+frac
1082
1083	IF (is_tree(ivm)) THEN
1084	IF (.NOT. is_tree(start_index(ivma))) THEN
1085
1086	IF (is_grassland_manag(start_index(ivma))) THEN
1087	instant_loss = 0.67 !forest to pasture
1088	ELSE IF (.NOT. natural(start_index(ivma))) THEN
1089	instant_loss = 0.5 !forest to crop
1090	ELSE
1091	instant_loss = 0. !forest to natural grassland
1092	ENDIF
1093
1094	CALL clear_forest (npts,ipts,ivm,biomass,frac, &
1095	instant_loss, &
1096	litter, deforest_biomass_remain,&
1097	fuel_1hr,fuel_10hr,&
1098	fuel_100hr,fuel_1000hr,&
1099	lignin_struc,&
1100	bm_to_litter_pro,convflux(:,iwplcc),prod10(:,:,iwplcc),prod100(:,:,iwplcc),&
1101	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
1102	fuel_1000hr_pro, lignin_content_pro)
1103	ENDIF
1104	ELSE
1105	!CALL harvest_herb(ipts,ivm,biomass,frac, &
1106	! bm_to_litter_pro)
1107	![2016-04-19] We put the transfer of biomass to litter for herbaceous
1108	! PFT directly here, because seprating them in a module harvest_herb
1109	! gives some error.
1110	bm_to_litter_pro(:,:) = bm_to_litter_pro(:,:) + biomass(ipts,ivm,:,:)*frac
1111
1112	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
1113	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
1114	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
1115	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
1116	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
1117	!don't forget to hanle litter lignin content
1118	lignin_content_pro(:)= lignin_content_pro(:) + &
1119	litter(ipts,istructural,ivm,:,icarbon)lignin_struc(ipts,ivm,:)frac
1120	ENDIF
1121
1122	!! scalar variables to be accumulated and inherited
1123	!! by the destination PFT
1124	bm_to_litter_pro(:,:) = bm_to_litter_pro(:,:) + &
1125	bm_to_litter(ipts,ivm,:,:)*frac
1126	carbon_pro(:) = carbon_pro(:)+carbon(ipts,:,ivm)*frac
1127	deepC_a_pro(:) = deepC_a_pro(:)+deepC_a(ipts,:,ivm)*frac
1128	deepC_s_pro(:) = deepC_s_pro(:)+deepC_s(ipts,:,ivm)*frac
1129	deepC_p_pro(:) = deepC_p_pro(:)+deepC_p(ipts,:,ivm)*frac
1130	co2_to_bm_pro = co2_to_bm_pro + co2_to_bm(ipts,ivm)*frac
1131
1132	gpp_daily_pro = gpp_daily_pro + gpp_daily(ipts,ivm)*frac
1133	npp_daily_pro = npp_daily_pro + npp_daily(ipts,ivm)*frac
1134	resp_maint_pro = resp_maint_pro + resp_maint(ipts,ivm)*frac
1135	resp_growth_pro = resp_growth_pro + resp_growth(ipts,ivm)*frac
1136	resp_hetero_pro = resp_hetero_pro + resp_hetero(ipts,ivm)*frac
1137	co2_fire_pro = co2_fire_pro + co2_fire(ipts,ivm)*frac
1138	ENDIF
1139	ENDDO
1140
1141	WHERE (litter_pro(istructural,:,icarbon) .GT. min_stomate)
1142	lignin_struc_pro(:) = lignin_content_pro(:)/litter_pro(istructural,:,icarbon)
1143	ELSEWHERE
1144	lignin_struc_pro(:) = zero
1145	ENDWHERE
1146
1147	END SUBROUTINE collect_legacy_pft
1148
1149	! ================================================================================================================================
1150	!! SUBROUTINE collect_legacy_pft_forestry
1151	!!
1152	!>\BRIEF : Collect the legacy variables that are going to be included
1153	! in the newly initialized PFT.
1154	!!
1155	!>\DESCRIPTION
1156	!_ ================================================================================================================================
1157	SUBROUTINE collect_legacy_pft_forestry(npts, ipts, ivma, glcc_pftmtc, &
1158	fuelfrac, &
1159	biomass, bm_to_litter, carbon, litter, &
1160	deepC_a, deepC_s, deepC_p, &
1161	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1162	lignin_struc, co2_to_bm, gpp_daily, npp_daily, &
1163	resp_maint, resp_growth, resp_hetero, co2_fire, &
1164	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1165	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1166	deforest_litter_remain, deforest_biomass_remain, &
1167	veget_max_pro, carbon_pro, lignin_struc_pro, litter_pro, &
1168	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
1169	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
1170	bm_to_litter_pro, co2_to_bm_pro, gpp_daily_pro, &
1171	npp_daily_pro, resp_maint_pro, resp_growth_pro, &
1172	resp_hetero_pro, co2_fire_pro, &
1173	convflux,prod10,prod100)
1174
1175	IMPLICIT NONE
1176
1177	!! 0.1 Input variables
1178	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
1179	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
1180	INTEGER, INTENT(in) :: ivma !! Index for metaclass
1181	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
1182	REAL(r_std), DIMENSION(:), INTENT(in) :: fuelfrac !!
1183	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
1184	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: bm_to_litter !! Transfer of biomass to litter
1185	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1186	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: carbon !! carbon pool: active, slow, or passive
1187	!! @tex ($gC m^{-2}$) @endtex
1188	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_a !! Permafrost soil carbon (g/m**3) active
1189	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1190	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1191	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: litter !! metabolic and structural litter, above and
1192	!! below ground @tex ($gC m^{-2}$) @endtex
1193	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1hr
1194	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_10hr
1195	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_100hr
1196	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: fuel_1000hr
1197	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1198	!! above and below ground
1199	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_to_bm !! biomass uptaken
1200	!! @tex ($gC m^{-2} day^{-1}$) @endtex
1201	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gpp_daily !! Daily gross primary productivity
1202	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1203	REAL(r_std), DIMENSION(:,:), INTENT(in) :: npp_daily !! Net primary productivity
1204	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1205	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_maint !! Maintenance respiration
1206	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1207	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_growth !! Growth respiration
1208	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1209	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_hetero !! Heterotrophic respiration
1210	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1211	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_fire !! Heterotrophic respiration
1212	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1213	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_1hr_remain
1214	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_10hr_remain
1215	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_100hr_remain
1216	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: def_fuel_1000hr_remain
1217	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
1218	!! deforestation region.
1219	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
1220	!! deforestation region.
1221
1222	!! 0.2 Output variables
1223	REAL(r_std), DIMENSION(:), INTENT(inout) :: carbon_pro
1224	REAL(r_std), DIMENSION(:), INTENT(inout) :: deepC_a_pro
1225	REAL(r_std), DIMENSION(:), INTENT(inout) :: deepC_s_pro
1226	REAL(r_std), DIMENSION(:), INTENT(inout) :: deepC_p_pro
1227	REAL(r_std), DIMENSION(:), INTENT(inout) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
1228	!! above and below ground
1229	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: litter_pro
1230	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1hr_pro
1231	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_10hr_pro
1232	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_100hr_pro
1233	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: fuel_1000hr_pro
1234	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bm_to_litter_pro
1235	REAL(r_std), INTENT(inout) :: veget_max_pro, co2_to_bm_pro
1236	REAL(r_std), INTENT(inout) :: gpp_daily_pro, npp_daily_pro
1237	REAL(r_std), INTENT(inout) :: resp_maint_pro, resp_growth_pro
1238	REAL(r_std), INTENT(inout) :: resp_hetero_pro, co2_fire_pro
1239
1240	!! 0.3 Modified variables
1241	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: convflux !! release during first year following land cover
1242	!! change
1243
1244	REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! products remaining in the 10 year-turnover
1245	!! pool after the annual release for each
1246	!! compartment (10 + 1 : input from year of land
1247	!! cover change)
1248	REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! products remaining in the 100 year-turnover
1249	!! pool after the annual release for each
1250	!! compartment (100 + 1 : input from year of land
1251	!! cover change)
1252
1253	!! 0.4 Local variables
1254	REAL(r_std), DIMENSION(nlevs) :: lignin_content_pro
1255	REAL(r_std) :: frac,pfuelfrac
1256	INTEGER :: ivm
1257
1258
1259	! All *_pro variables collect the legacy pools/fluxes of the ancestor
1260	! PFTs for the receiving youngest age class. All *_pro variables
1261	! represent the quantity weighted by the fraction of ancestor contributing
1262	! PFTs.
1263	! Exceptions:
1264	! lignin_struc_pro:: the ratio of lignin content in structural litter.
1265
1266	lignin_content_pro(:)=zero
1267	pfuelfrac = fuelfrac(ipts)
1268
1269	DO ivm = 1,nvm
1270	frac = glcc_pftmtc(ipts,ivm,ivma)
1271	IF (frac>zero) THEN
1272	veget_max_pro = veget_max_pro+frac
1273
1274	IF (is_tree(ivm)) THEN
1275	IF (is_tree(start_index(ivma))) THEN
1276	CALL harvest_industrial (npts,ipts,ivm,biomass,frac*(1-pfuelfrac), &
1277	litter, deforest_biomass_remain,&
1278	fuel_1hr,fuel_10hr,&
1279	fuel_100hr,fuel_1000hr,&
1280	lignin_struc,&
1281	bm_to_litter_pro,convflux(:,iwphar),prod10(:,:,iwphar),prod100(:,:,iwphar),&
1282	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
1283	fuel_1000hr_pro, lignin_content_pro)
1284
1285	CALL harvest_fuelwood (npts,ipts,ivm,biomass,frac*pfuelfrac, &
1286	litter, deforest_biomass_remain,&
1287	fuel_1hr,fuel_10hr,&
1288	fuel_100hr,fuel_1000hr,&
1289	lignin_struc,&
1290	bm_to_litter_pro,convflux(:,iwphar),prod10(:,:,iwphar),prod100(:,:,iwphar),&
1291	litter_pro, fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, &
1292	fuel_1000hr_pro, lignin_content_pro)
1293	ENDIF
1294	ELSE
1295	!CALL harvest_herb(ipts,ivm,biomass,frac, &
1296	! bm_to_litter_pro)
1297	![2016-04-19] We put the transfer of biomass to litter for herbaceous
1298	! PFT directly here, because seprating them in a module harvest_herb
1299	! gives some error.
1300	bm_to_litter_pro(:,:) = bm_to_litter_pro(:,:) + biomass(ipts,ivm,:,:)*frac
1301
1302	litter_pro(:,:,:) = litter_pro(:,:,:) + litter(ipts,:,ivm,:,:)*frac
1303	fuel_1hr_pro(:,:) = fuel_1hr_pro(:,:) + fuel_1hr(ipts,ivm,:,:)*frac
1304	fuel_10hr_pro(:,:) = fuel_10hr_pro(:,:) + fuel_10hr(ipts,ivm,:,:)*frac
1305	fuel_100hr_pro(:,:) = fuel_100hr_pro(:,:) + fuel_100hr(ipts,ivm,:,:)*frac
1306	fuel_1000hr_pro(:,:) = fuel_1000hr_pro(:,:) + fuel_1000hr(ipts,ivm,:,:)*frac
1307	!don't forget to hanle litter lignin content
1308	lignin_content_pro(:)= lignin_content_pro(:) + &
1309	litter(ipts,istructural,ivm,:,icarbon)lignin_struc(ipts,ivm,:)frac
1310	ENDIF
1311
1312	!! scalar variables to be accumulated and inherited
1313	!! by the destination PFT
1314	bm_to_litter_pro(:,:) = bm_to_litter_pro(:,:) + &
1315	bm_to_litter(ipts,ivm,:,:)*frac
1316	carbon_pro(:) = carbon_pro(:)+carbon(ipts,:,ivm)*frac
1317	deepC_a_pro(:) = deepC_a_pro(:)+deepC_a(ipts,:,ivm)*frac
1318	deepC_s_pro(:) = deepC_s_pro(:)+deepC_s(ipts,:,ivm)*frac
1319	deepC_p_pro(:) = deepC_p_pro(:)+deepC_p(ipts,:,ivm)*frac
1320	co2_to_bm_pro = co2_to_bm_pro + co2_to_bm(ipts,ivm)*frac
1321
1322	gpp_daily_pro = gpp_daily_pro + gpp_daily(ipts,ivm)*frac
1323	npp_daily_pro = npp_daily_pro + npp_daily(ipts,ivm)*frac
1324	resp_maint_pro = resp_maint_pro + resp_maint(ipts,ivm)*frac
1325	resp_growth_pro = resp_growth_pro + resp_growth(ipts,ivm)*frac
1326	resp_hetero_pro = resp_hetero_pro + resp_hetero(ipts,ivm)*frac
1327	co2_fire_pro = co2_fire_pro + co2_fire(ipts,ivm)*frac
1328	ENDIF
1329	ENDDO
1330
1331	WHERE (litter_pro(istructural,:,icarbon) .GT. min_stomate)
1332	lignin_struc_pro(:) = lignin_content_pro(:)/litter_pro(istructural,:,icarbon)
1333	ELSEWHERE
1334	lignin_struc_pro(:) = zero
1335	ENDWHERE
1336
1337	END SUBROUTINE collect_legacy_pft_forestry
1338
1339
1340	! ================================================================================================================================
1341	!! SUBROUTINE : add_incoming_proxy_pft
1342	!!
1343	!>\BRIEF : Merge the newly incoming proxy PFT cohort with the exisiting
1344	!! cohort.
1345	!! \n
1346	!
1347	!_ ================================================================================================================================
1348	SUBROUTINE add_incoming_proxy_pft(npts, ipts, ipft, veget_max_pro, &
1349	carbon_pro, litter_pro, lignin_struc_pro, bm_to_litter_pro, &
1350	deepC_a_pro, deepC_s_pro, deepC_p_pro, &
1351	fuel_1hr_pro, fuel_10hr_pro, fuel_100hr_pro, fuel_1000hr_pro, &
1352	biomass_pro, co2_to_bm_pro, npp_longterm_pro, ind_pro, &
1353	lm_lastyearmax_pro, age_pro, everywhere_pro, &
1354	leaf_frac_pro, leaf_age_pro, PFTpresent_pro, senescence_pro, &
1355	gpp_daily_pro, npp_daily_pro, resp_maint_pro, resp_growth_pro, &
1356	resp_hetero_pro, co2_fire_pro, &
1357	veget_max, carbon, litter, lignin_struc, bm_to_litter, &
1358	deepC_a, deepC_s, deepC_p, &
1359	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1360	biomass, co2_to_bm, npp_longterm, ind, &
1361	lm_lastyearmax, age, everywhere, &
1362	leaf_frac, leaf_age, PFTpresent, senescence, &
1363	gpp_daily, npp_daily, resp_maint, resp_growth, &
1364	resp_hetero, co2_fire)
1365
1366	IMPLICIT NONE
1367
1368	!! 0.1 Input variables
1369	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
1370	INTEGER, INTENT(in) :: ipts !! Domain size - number of pixels (unitless)
1371	INTEGER, INTENT(in) :: ipft
1372	REAL(r_std), INTENT(in) :: veget_max_pro !! The land fraction of incoming new PFTs that are
1373	!! the sum of all its ancestor PFTs
1374	REAL(r_std), DIMENSION(:), INTENT(in) :: carbon_pro
1375	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_a_pro
1376	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_s_pro
1377	REAL(r_std), DIMENSION(:), INTENT(in) :: deepC_p_pro
1378	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: litter_pro
1379	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_1hr_pro
1380	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_10hr_pro
1381	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_100hr_pro
1382	REAL(r_std), DIMENSION(:,:), INTENT(in) :: fuel_1000hr_pro
1383	REAL(r_std), DIMENSION(:,:), INTENT(in) :: bm_to_litter_pro
1384	REAL(r_std), DIMENSION(:), INTENT(in) :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
1385	!! above and below ground
1386	REAL(r_std), DIMENSION(:,:), INTENT(in) :: biomass_pro !! biomass @tex ($gC m^{-2}$) @endtex
1387	REAL(r_std), DIMENSION(:), INTENT(in) :: leaf_frac_pro !! fraction of leaves in leaf age class
1388	REAL(r_std), DIMENSION(:), INTENT(in) :: leaf_age_pro !! fraction of leaves in leaf age class
1389	REAL(r_std), INTENT(in) :: ind_pro, age_pro, lm_lastyearmax_pro
1390	REAL(r_std), INTENT(in) :: npp_longterm_pro, co2_to_bm_pro
1391	REAL(r_std), INTENT(in) :: everywhere_pro !! is the PFT everywhere in the grid box or very
1392	LOGICAL, INTENT(in) :: PFTpresent_pro, senescence_pro !! Is pft there (unitless)
1393
1394	REAL(r_std), INTENT(in) :: gpp_daily_pro, npp_daily_pro
1395	REAL(r_std), INTENT(in) :: resp_maint_pro, resp_growth_pro
1396	REAL(r_std), INTENT(in) :: resp_hetero_pro, co2_fire_pro
1397
1398	!! 0.2 Output variables
1399
1400	!! 0.3 Modified variables
1401
1402	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1403	!! May sum to
1404	!! less than unity if the pixel has
1405	!! nobio area. (unitless, 0-1)
1406
1407	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
1408	!! @tex ($gC m^{-2}$) @endtex
1409	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
1410	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1411	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1412	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
1413	!! below ground @tex ($gC m^{-2}$) @endtex
1414	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
1415	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
1416	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
1417	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
1418	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1419	!! above and below ground
1420	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
1421	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1422	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
1423	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
1424	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
1425
1426	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
1427	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
1428	!! @tex $(m^{-2})$ @endtex
1429	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
1430	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
1431	!! each pixel
1432	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
1433	!! for deciduous trees)
1434	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
1435	!! @tex ($gC m^{-2}$) @endtex
1436	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
1437	!! very localized (after its introduction) (?)
1438	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
1439	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
1440
1441	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
1442	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1443	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
1444	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1445	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
1446	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1447	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
1448	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1449	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
1450	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1451	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_fire !! Heterotrophic respiration
1452	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1453
1454	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
1455	! !! unitless)
1456	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
1457	! !! (0 to 1, unitless)
1458	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
1459	! !! @tex $(gC m^{-2} day^{-1})$ @endtex
1460	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
1461	! !! -5 deg C (for phenology)
1462	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
1463	! !! the growing season (days)
1464	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
1465	! !! availability (days)
1466	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
1467	! !! (for phenology) - this is written to the
1468	! !! history files
1469	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
1470	! !! for crops
1471	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
1472	! !! C (for phenology)
1473	! REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
1474	! !! leaves were lost (for phenology)
1475
1476	!! 0.4 Local variables
1477
1478	INTEGER(i_std) :: iele !! Indeces(unitless)
1479	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
1480	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements) :: litter_old !! metabolic and structural litter, above and
1481	!! below ground @tex ($gC m^{-2}$) @endtex
1482	REAL(r_std) :: veget_old,veget_total
1483
1484
1485	! Back up some variables in case they're needed later
1486	litter_old(:,:,:,:,:) = litter(:,:,:,:,:)
1487
1488	!! General idea
1489	! The established proxy vegetation has a fraction of 'veget_max_pro'; the
1490	! existing iPFT has a fraction of veget_max(ipts,ipft).
1491	! Suppose we want to merge a scalar variable B, the value of B after merging
1492	! is (BiVi+BjVj)/(Vi+Vj), where Vi is the original veget_max, Vj is the
1493	! incoming veget_max. Note that in case Vi=0, this equation remains solid,
1494	! i.e. the veget_max after merging is Vj and B after merging is Bj. In other
1495	! words, the proxy vegetation "fills" up the empty niche of iPFT.
1496	! Also note that for many scalar variables our input value is Bj*Vj, which
1497	! is accumulated from multiple ancestor PFTs.
1498	veget_old = veget_max(ipts,ipft)
1499	veget_total = veget_old+veget_max_pro
1500
1501	!! Different ways of handling merging depending on nature of variables:
1502
1503	!! 1. Area-based scalar variables, use the equation above
1504	! biomass,carbon, litter, bm_to_litter, co2_to_bm, ind,
1505	! lm_lastyearmax, npp_longterm, lm_lastyearmax,
1506	! lignin_struc (ratio variable depending on area-based variable)
1507
1508	!! 2. Variables are tentatively handled like area-based variables:
1509	! leaf_frac, leaf_age,
1510
1511	!! 3. Variables that are overwritten by the newly initialized PFT:
1512	! PFTpresent, senescence
1513
1514	!! 4. Variables whose operation is uncertain and are not handled currently:
1515	! when_growthinit :: how many days ago was the beginning of the growing season (days)
1516	! gdd_from_growthinit :: growing degree days, since growthinit
1517	! gdd_midwinter, time_hum_min, gdd_m5_dormance, ncd_dormance,
1518	! moiavail_month, moiavail_week, ngd_minus5
1519
1520	!! 5. Variables that concern with short-term fluxes that do not apply in
1521	! this case:
1522	! gpp_daily, npp_daily etc.
1523
1524	! Add the coming veget_max_pro into existing veget_max
1525	veget_max(ipts,ipft) = veget_total
1526
1527	IF (veget_total .GT. min_stomate) THEN
1528	! Merge scalar variables which are defined on area basis
1529	carbon(ipts,:,ipft) = (veget_old * carbon(ipts,:,ipft) + &
1530	carbon_pro(:))/veget_total
1531	deepC_a(ipts,:,ipft) = (veget_old * deepC_a(ipts,:,ipft) + &
1532	deepC_a_pro(:))/veget_total
1533	deepC_s(ipts,:,ipft) = (veget_old * deepC_s(ipts,:,ipft) + &
1534	deepC_s_pro(:))/veget_total
1535	deepC_p(ipts,:,ipft) = (veget_old * deepC_p(ipts,:,ipft) + &
1536	deepC_p_pro(:))/veget_total
1537	litter(ipts,:,ipft,:,:) = (veget_old * litter(ipts,:,ipft,:,:) + &
1538	litter_pro(:,:,:))/veget_total
1539	fuel_1hr(ipts,ipft,:,:) = (veget_old * fuel_1hr(ipts,ipft,:,:) + &
1540	fuel_1hr_pro(:,:))/veget_total
1541	fuel_10hr(ipts,ipft,:,:) = (veget_old * fuel_10hr(ipts,ipft,:,:) + &
1542	fuel_10hr_pro(:,:))/veget_total
1543	fuel_100hr(ipts,ipft,:,:) = (veget_old * fuel_100hr(ipts,ipft,:,:) + &
1544	fuel_100hr_pro(:,:))/veget_total
1545	fuel_1000hr(ipts,ipft,:,:) = (veget_old * fuel_1000hr(ipts,ipft,:,:) + &
1546	fuel_1000hr_pro(:,:))/veget_total
1547
1548	WHERE (litter(ipts,istructural,ipft,:,icarbon) .GT. min_stomate)
1549	lignin_struc(ipts,ipft,:) = (veget_oldlitter_old(ipts,istructural,ipft,:,icarbon) &
1550	lignin_struc(ipts,ipft,:) + litter_pro(istructural,:,icarbon)* &
1551	lignin_struc_pro(:))/(veget_total*litter(ipts,istructural,ipft,:,icarbon))
1552	ENDWHERE
1553	bm_to_litter(ipts,ipft,:,:) = (veget_old * bm_to_litter(ipts,ipft,:,:) + &
1554	bm_to_litter_pro(:,:))/veget_total
1555
1556	biomass(ipts,ipft,:,:) = (biomass(ipts,ipft,:,:)*veget_old + &
1557	biomass_pro(:,:))/veget_total
1558	co2_to_bm(ipts,ipft) = (veget_old*co2_to_bm(ipts,ipft) + &
1559	co2_to_bm_pro)/veget_total
1560	ind(ipts,ipft) = (ind(ipts,ipft)*veget_old + ind_pro)/veget_total
1561	lm_lastyearmax(ipts,ipft) = (lm_lastyearmax(ipts,ipft)*veget_old + &
1562	lm_lastyearmax_pro)/veget_total
1563	npp_longterm(ipts,ipft) = (veget_old * npp_longterm(ipts,ipft) + &
1564	npp_longterm_pro)/veget_total
1565
1566	!CHECK: Here follows the original idea in DOFOCO, more strictly,
1567	! leas mass should be considered together. The same also applies on
1568	! leaf age.
1569	leaf_frac(ipts,ipft,:) = (leaf_frac(ipts,ipft,:)*veget_old + &
1570	leaf_frac_pro(:))/veget_total
1571	leaf_age(ipts,ipft,:) = (leaf_age(ipts,ipft,:)*veget_old + &
1572	leaf_age_pro(:))/veget_total
1573	age(ipts,ipft) = (veget_old * age(ipts,ipft) + &
1574	age_pro)/veget_total
1575
1576	! Everywhere deals with the migration of vegetation. Copy the
1577	! status of the most migrated vegetation for the whole PFT
1578	everywhere(ipts,ipft) = MAX(everywhere(ipts,ipft), everywhere_pro)
1579
1580	! Overwrite the original variables with that from newly initialized
1581	! proxy PFT
1582	PFTpresent(ipts,ipft) = PFTpresent_pro
1583	senescence(ipts,ipft) = senescence_pro
1584
1585	! This is to close carbon loop when writing history variables.
1586	gpp_daily(ipts,ipft) = (veget_old * gpp_daily(ipts,ipft) + &
1587	gpp_daily_pro)/veget_total
1588	npp_daily(ipts,ipft) = (veget_old * npp_daily(ipts,ipft) + &
1589	npp_daily_pro)/veget_total
1590	resp_maint(ipts,ipft) = (veget_old * resp_maint(ipts,ipft) + &
1591	resp_maint_pro)/veget_total
1592	resp_growth(ipts,ipft) = (veget_old * resp_growth(ipts,ipft) + &
1593	resp_growth_pro)/veget_total
1594	resp_hetero(ipts,ipft) = (veget_old * resp_hetero(ipts,ipft) + &
1595	resp_hetero_pro)/veget_total
1596	co2_fire(ipts,ipft) = (veget_old * co2_fire(ipts,ipft) + &
1597	co2_fire_pro)/veget_total
1598
1599	! Phenology- or time-related variables will be copied from original values if
1600	! there is already youngest-age-class PFT there, otherwise they're left
1601	! untouched, because 1. to initiliaze all new PFTs here is wrong and
1602	! phenology is not explicitly considered, so we cannot assign a value
1603	! to these variables. 2. We assume they will be correctly filled if
1604	! other variables are in place (e.g., non-zero leaf mass will lead to
1605	! onset of growing season). In this case, merging a newly initialized PFT
1606	! to an existing one is not the same as merging PFTs when they grow
1607	! old enough to exceed thresholds.
1608
1609	! gpp_week(ipts,ipft) = (veget_old * gpp_week(ipts,ipft) + &
1610	! gpp_week_pro)/veget_total
1611	! when_growthinit(ipts,ipft) = (veget_old * when_growthinit(ipts,ipft) + &
1612	! when_growthinit_pro)/veget_total
1613	! gdd_from_growthinit(ipts,ipft) = (veget_old * gdd_from_growthinit(ipts,ipft) + &
1614	! gdd_from_growthinit_pro)/veget_total
1615	! gdd_midwinter(ipts,ipft) = (veget_old * gdd_midwinter(ipts,ipft) + &
1616	! gdd_midwinter_pro)/veget_total
1617	! time_hum_min(ipts,ipft) = (veget_old * time_hum_min(ipts,ipft) + &
1618	! time_hum_min_pro)/veget_total
1619	! gdd_m5_dormance(ipts,ipft) = (veget_old * gdd_m5_dormance(ipts,ipft) + &
1620	! gdd_m5_dormance_pro)/veget_total
1621	! ncd_dormance(ipts,ipft) = (veget_old * ncd_dormance(ipts,ipft) + &
1622	! ncd_dormance_pro)/veget_total
1623	! moiavail_month(ipts,ipft) = (veget_old * moiavail_month(ipts,ipft) + &
1624	! moiavail_month_pro)/veget_total
1625	! moiavail_week(ipts,ipft) = (veget_old * moiavail_week(ipts,ipft) + &
1626	! moiavail_week_pro)/veget_total
1627	! ngd_minus5(ipts,ipft) = (veget_old * ngd_minus5(ipts,ipft) + &
1628	! ngd_minus5_pro)/veget_total
1629	ENDIF
1630
1631
1632	END SUBROUTINE add_incoming_proxy_pft
1633
1634	! ================================================================================================================================
1635	!! SUBROUTINE : empty_pft
1636	!!
1637	!>\BRIEF : Empty a PFT when,
1638	!! - it is exhausted because of land cover change.
1639	!! - it moves to the next age class
1640	!! \n
1641	!_ ================================================================================================================================
1642	SUBROUTINE empty_pft(ipts, ivm, veget_max, biomass, ind, &
1643	carbon, litter, lignin_struc, bm_to_litter, &
1644	deepC_a, deepC_s, deepC_p, &
1645	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1646	gpp_daily, npp_daily, gpp_week, npp_longterm, &
1647	co2_to_bm, resp_maint, resp_growth, resp_hetero, &
1648	lm_lastyearmax, leaf_frac, leaf_age, age, &
1649	everywhere, PFTpresent, when_growthinit, &
1650	senescence, gdd_from_growthinit, gdd_midwinter, &
1651	time_hum_min, gdd_m5_dormance, ncd_dormance, &
1652	moiavail_month, moiavail_week, ngd_minus5)
1653
1654	IMPLICIT NONE
1655
1656	!! 0.1 Input variables
1657	INTEGER, INTENT(in) :: ipts !! index for grid cell
1658	INTEGER, INTENT(in) :: ivm !! index for pft
1659
1660	!! 0.2 Output variables
1661
1662	!! 0.3 Modified variables
1663
1664	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
1665	!! May sum to
1666	!! less than unity if the pixel has
1667	!! nobio area. (unitless, 0-1)
1668	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
1669	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Number of individuals at the stand level
1670	!! @tex $(m^{-2})$ @endtex
1671	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
1672	!! @tex ($gC m^{-2}$) @endtex
1673	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
1674	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
1675	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
1676	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
1677	!! below ground @tex ($gC m^{-2}$) @endtex
1678	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1hr
1679	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_10hr
1680	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_100hr
1681	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: fuel_1000hr
1682	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
1683	!! above and below ground
1684	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
1685	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1686	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
1687	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1688	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
1689	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1690	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
1691	!! @tex $(gC m^{-2} day^{-1})$ @endtex
1692	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary productivity
1693	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
1694	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
1695	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
1696	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1697	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
1698	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1699	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
1700	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1701	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
1702	!! @tex ($gC m^{-2}$) @endtex
1703	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless;0-1)
1704	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
1705	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
1706	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
1707	!! very localized (after its introduction) (?)
1708	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
1709	!! each pixel
1710	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
1711	!! the growing season (days)
1712	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
1713	!! for deciduous trees)
1714	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
1715	!! for crops
1716	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
1717	!! (for phenology) - this is written to the
1718	!! history files
1719	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
1720	!! availability (days)
1721	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
1722	!! C (for phenology)
1723	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
1724	!! leaves were lost (for phenology)
1725	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
1726	!! unitless)
1727	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
1728	!! (0 to 1, unitless)
1729	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
1730	!! -5 deg C (for phenology)
1731
1732	!! 0.4 Local variables
1733	INTEGER(i_std) :: iele !! Indeces(unitless)
1734	INTEGER(i_std) :: ilit,ilev,icarb !! Indeces(unitless)
1735
1736	veget_max(ipts,ivm) = zero
1737	ind(ipts,ivm) = zero
1738	biomass(ipts,ivm,:,:) = zero
1739	litter(ipts,:,ivm,:,:) = zero
1740	fuel_1hr(ipts,ivm,:,:) = zero
1741	fuel_10hr(ipts,ivm,:,:) = zero
1742	fuel_100hr(ipts,ivm,:,:) = zero
1743	fuel_1000hr(ipts,ivm,:,:) = zero
1744	carbon(ipts,:,ivm) = zero
1745	deepC_a(ipts,:,ivm) = zero
1746	deepC_s(ipts,:,ivm) = zero
1747	deepC_p(ipts,:,ivm) = zero
1748	bm_to_litter(ipts,ivm,:,:) = zero
1749	DO ilev=1,nlevs
1750	lignin_struc(ipts,ivm,ilev) = zero
1751	ENDDO
1752	npp_longterm(ipts,ivm) = zero
1753	gpp_daily(ipts,ivm) = zero
1754	gpp_week(ipts,ivm) = zero
1755	resp_maint(ipts,ivm) = zero
1756	resp_growth(ipts,ivm) = zero
1757	resp_hetero(ipts,ivm) = zero
1758	npp_daily(ipts,ivm) = zero
1759	co2_to_bm(ipts,ivm) = zero
1760	lm_lastyearmax(ipts,ivm) = zero
1761	age(ipts,ivm) = zero
1762	leaf_frac(ipts,ivm,:) = zero
1763	leaf_age(ipts,ivm,:) = zero
1764	everywhere(ipts,ivm) = zero
1765	when_growthinit(ipts,ivm) = zero
1766	gdd_from_growthinit(ipts,ivm) = zero
1767	gdd_midwinter(ipts,ivm) = zero
1768	time_hum_min(ipts,ivm) = zero
1769	gdd_m5_dormance(ipts,ivm) = zero
1770	ncd_dormance(ipts,ivm) = zero
1771	moiavail_month(ipts,ivm) = zero
1772	moiavail_week(ipts,ivm) = zero
1773	ngd_minus5(ipts,ivm) = zero
1774	PFTpresent(ipts,ivm) = .FALSE.
1775	senescence(ipts,ivm) = .FALSE.
1776
1777	END SUBROUTINE empty_pft
1778
1779	SUBROUTINE prepare_balance_check(outflux_sta,influx_sta,pool_sta, &
1780	veget_cov_max, &
1781	co2_to_bm,gpp_daily,npp_daily, &
1782	biomass,litter,carbon,prod10,prod100, &
1783	bm_to_litter,resp_maint,resp_growth,resp_hetero, &
1784	convflux,cflux_prod10,cflux_prod100,co2_fire)
1785
1786	IMPLICIT NONE
1787
1788	!! 0.1 Input variables
1789	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_cov_max !! "Maximal" coverage fraction of a PFT (LAI
1790	!! -> infinity) on ground
1791	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_to_bm !! CO2 taken up from atmosphere when
1792	!! introducing a new PFT (introduced for
1793	!! carbon balance closure)
1794	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1795	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gpp_daily !! Daily gross primary productivity
1796	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1797	REAL(r_std), DIMENSION(:,:), INTENT(in) :: npp_daily !! Net primary productivity
1798	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1799	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! Biomass @tex $(gC m^{-2})$ @endtex
1800	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in):: litter !! Metabolic and structural litter, above
1801	!! and below ground
1802	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: carbon !! Carbon pool: active, slow, or passive,
1803	!! @tex $(gC m^{-2})$ @endtex
1804	REAL(r_std),DIMENSION(:,:,:), INTENT(in) :: prod10 !! Products remaining in the 10
1805	!! year-turnover pool after the annual
1806	!! release for each compartment (10
1807	!! + 1 : input from year of land cover
1808	!! change) @tex $(gC m^{-2})$ @endtex
1809	REAL(r_std),DIMENSION(:,:,:), INTENT(in) :: prod100 !! Products remaining in the 100
1810	!! year-turnover pool after the annual
1811	!! release for each compartment (100
1812	!! + 1 : input from year of land cover
1813	!! change) @tex $(gC m^{-2})$ @endtex
1814	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in):: bm_to_litter !! Conversion of biomass to litter
1815	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1816	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_maint !! Maintenance respiration
1817	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1818
1819	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_growth !! Growth respiration
1820	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1821	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_hetero !! Heterotrophic respiration
1822	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1823	REAL(r_std),DIMENSION(:,:), INTENT(in) :: convflux !! Release during first year following land
1824	!! cover change @tex $(gC m^{-2})$ @endtex
1825	REAL(r_std),DIMENSION(:,:), INTENT(in) :: cflux_prod10 !! Total annual release from the 10
1826	!! year-turnover pool
1827	!! @tex $(gC m^{-2})$ @endtex
1828	REAL(r_std),DIMENSION(:,:), INTENT(in) :: cflux_prod100 !! Total annual release from the 100
1829	!! year-turnover pool
1830	!! @tex $(gC m^{-2})$ @endtex
1831	REAL(r_std), DIMENSION(:,:), INTENT(in) :: co2_fire !! Carbon emitted into the atmosphere by
1832	!! fire (living and dead biomass)
1833
1834
1835	!! 0.2 Output variables
1836	REAL(r_std), DIMENSION(:,:) :: outflux_sta
1837	REAL(r_std), DIMENSION(:,:) :: influx_sta
1838	REAL(r_std), DIMENSION(:,:) :: pool_sta
1839
1840
1841	!! 0.3 Modified variables
1842
1843	!! 0.4 Local variables
1844	INTEGER(i_std) :: ind_biomass,ind_litter,ind_soil,ind_prod,ind_co2tobm,ind_gpp,ind_npp,&
1845	ind_bm2lit,ind_resph,ind_respm,ind_respg,ind_convf,ind_cflux,ind_fire
1846
1847
1848
1849	ind_biomass = 1
1850	ind_litter = 2
1851	ind_soil = 3
1852	ind_prod = 4
1853	pool_sta(:,ind_biomass) = SUM(SUM(biomass(:,:,:,icarbon),DIM=3) * veget_cov_max,DIM=2)
1854	pool_sta(:,ind_litter) = SUM(SUM(SUM(litter(:,:,:,:,icarbon),DIM=4),DIM=2) * veget_cov_max,DIM=2)
1855	pool_sta(:,ind_soil) = SUM(SUM(carbon(:,:,:),DIM=2) * veget_cov_max,DIM=2)
1856	pool_sta(:,ind_prod) = SUM(SUM(prod10,DIM=3),DIM=2) + SUM(SUM(prod100,DIM=3),DIM=2)
1857
1858	ind_co2tobm = 1
1859	ind_gpp = 2
1860	ind_npp = 3
1861	influx_sta(:,ind_co2tobm) = SUM(co2_to_bm*veget_cov_max,DIM=2)
1862	influx_sta(:,ind_gpp) = SUM(gpp_daily*veget_cov_max,DIM=2)
1863	influx_sta(:,ind_npp) = SUM(npp_daily*veget_cov_max,DIM=2)
1864
1865	ind_bm2lit = 1
1866	ind_respm = 2
1867	ind_respg = 3
1868	ind_resph = 4
1869	ind_convf = 5
1870	ind_cflux = 6
1871	ind_fire = 7
1872	outflux_sta(:,ind_bm2lit) = SUM(SUM(bm_to_litter(:,:,:,icarbon),DIM=3) * veget_cov_max,DIM=2)
1873	outflux_sta(:,ind_respm) = SUM(resp_maint*veget_cov_max,DIM=2)
1874	outflux_sta(:,ind_respg) = SUM(resp_growth*veget_cov_max,DIM=2)
1875	outflux_sta(:,ind_resph) = SUM(resp_hetero*veget_cov_max,DIM=2)
1876	outflux_sta(:,ind_convf) = SUM(convflux,DIM=2)
1877	outflux_sta(:,ind_cflux) = SUM(cflux_prod10,DIM=2) + SUM(cflux_prod100,DIM=2)
1878	outflux_sta(:,ind_fire) = SUM(co2_fire*veget_cov_max,DIM=2)
1879
1880	END SUBROUTINE prepare_balance_check
1881
1882
1883	SUBROUTINE luc_balance_check(outflux_sta,influx_sta,pool_sta, &
1884	outflux_end,influx_end,pool_end, &
1885	npts,lalo,identifier)
1886
1887	IMPLICIT NONE
1888
1889	!! 0.1 Input variables
1890	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
1891	REAL(r_std), DIMENSION(:,:), INTENT(in) :: outflux_sta
1892	REAL(r_std), DIMENSION(:,:), INTENT(in) :: influx_sta
1893	REAL(r_std), DIMENSION(:,:), INTENT(in) :: pool_sta
1894	REAL(r_std), DIMENSION(:,:), INTENT(in) :: outflux_end
1895	REAL(r_std), DIMENSION(:,:), INTENT(in) :: influx_end
1896	REAL(r_std), DIMENSION(:,:), INTENT(in) :: pool_end
1897	REAL(r_std), DIMENSION(:,:), INTENT(in) :: lalo !! Geographical coordinates (latitude,longitude)
1898	!! for pixels (degrees)
1899	CHARACTER(LEN=*), INTENT(in), OPTIONAL :: identifier !! string with to identify where this routine was called form
1900
1901	!! 0.2 Output variables
1902
1903	!! 0.3 Modified variables
1904
1905	!! 0.4 Local variables
1906	REAL(r_std), DIMENSION(npts,nelements) :: mass_before_2rd !! Temporary variable
1907	REAL(r_std), DIMENSION(npts,nelements) :: mass_after_2rd !! Temporary variable
1908	REAL(r_std), DIMENSION(npts,nelements) :: mass_change_2rd !! positive
1909	REAL(r_std), DIMENSION(npts,nelements) :: mass_balance_2rd !! Temporary variable
1910	INTEGER(i_std) :: ipts
1911
1912
1913	mass_before_2rd = zero
1914	mass_after_2rd = zero
1915	mass_change_2rd = zero
1916	mass_balance_2rd = zero
1917
1918	mass_before_2rd(:,icarbon) = SUM(pool_sta(:,:),DIM=2)
1919	mass_after_2rd(:,icarbon) = SUM(pool_end(:,:),dim=2)
1920
1921	! mass_change_2rd is the mass increase
1922	mass_change_2rd(:,icarbon) = SUM(influx_end(:,:),DIM=2) - SUM(influx_sta(:,:),DIM=2) + &
1923	+ SUM(outflux_sta(:,:),DIM=2) - SUM(outflux_end(:,:),DIM=2)
1924
1925	mass_balance_2rd(:,icarbon) = mass_before_2rd(:,icarbon) - mass_after_2rd(:,icarbon) &
1926	+ mass_change_2rd(:,icarbon)
1927
1928	DO ipts = 1,npts
1929	IF (ABS(mass_balance_2rd(ipts,icarbon)) .GE. 1e-2) THEN
1930	WRITE (numout,*) ' FATAL Error'
1931	WRITE (numout,*) ' Mass conservation failed after ',identifier
1932	WRITE (numout,*) ' limit: ' , 1e-2
1933	WRITE (numout,*) ' Mismatch :' , mass_balance_2rd(ipts,icarbon)
1934	WRITE (numout,*) ' Coordinates :' , lalo(ipts,:)
1935	WRITE (numout,*) ' gridpoint: ' , ipts , ' of ngrids: ',npts
1936	STOP
1937	ENDIF
1938	ENDDO
1939	WRITE (numout,*) 'mass balance check successful after ',identifier
1940
1941	END SUBROUTINE luc_balance_check
1942
1943
1944
1945
1946	! Note this subroutine does not depend on how many age classes there are
1947	! in different MTCs.
1948	SUBROUTINE glcc_compensation_full(npts,veget_4veg,glcc,glccReal,glccDef, &
1949	p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
1950	IncreDeficit)
1951
1952	IMPLICIT NONE
1953
1954	!! 0.1 Input variables
1955	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
1956	INTEGER, INTENT(in) :: p2c,ipasture,g2c,igrass,f2c,itree,icrop
1957	REAL(r_std), DIMENSION (npts,12),INTENT(in) :: glcc !! the land-cover-change (LCC) matrix in case a gross LCC is
1958	!! used.
1959
1960	!! 0.2 Output variables
1961
1962
1963	!! 0.3 Modified variables
1964	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: veget_4veg !! "maximal" coverage of tree/grass/pasture/crop
1965	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccDef !! Gross LCC deficit, negative values mean that there
1966	!! are not enough fractions in the source vegetations
1967	!! to the target ones as presribed by the LCC matrix.
1968	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccReal !! The "real" glcc matrix that we apply in the model
1969	!! after considering the consistency between presribed
1970	!! glcc matrix and existing vegetation fractions.
1971	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
1972	!! there are not enough fractions in the source PFTs
1973	!! /vegetations to target PFTs/vegetations. I.e., these
1974	!! fraction transfers are presribed in LCC matrix but
1975	!! not realized.
1976
1977	!! 0.4 Local variables
1978	REAL(r_std), DIMENSION(npts) :: tmpdef !! LCC deficits by summing up all the deficits to the
1979	!! the same target vegetation.
1980
1981	!! 0. We first handle the cases where veget_4veg might be very small
1982	!tree
1983	WHERE(veget_4veg(:,itree) > min_stomate)
1984	glccDef(:,f2c) = veget_4veg(:,itree)-glcc(:,f2c)
1985	WHERE(veget_4veg(:,itree)>glcc(:,f2c))
1986	glccReal(:,f2c) = glcc(:,f2c)
1987	ELSEWHERE
1988	glccReal(:,f2c) = veget_4veg(:,itree)
1989	ENDWHERE
1990	ELSEWHERE
1991	glccReal(:,f2c) = 0.
1992	glccDef(:,f2c) = -1*glcc(:,f2c)
1993	ENDWHERE
1994
1995	!pasture
1996	WHERE(veget_4veg(:,ipasture) > min_stomate)
1997	glccDef(:,p2c) = veget_4veg(:,ipasture)-glcc(:,p2c)
1998	WHERE(veget_4veg(:,ipasture)>glcc(:,p2c))
1999	glccReal(:,p2c) = glcc(:,p2c)
2000	ELSEWHERE
2001	glccReal(:,p2c) = veget_4veg(:,ipasture)
2002	ENDWHERE
2003	ELSEWHERE
2004	glccReal(:,p2c) = 0.
2005	glccDef(:,p2c) = -1*glcc(:,p2c)
2006	ENDWHERE
2007
2008	!grass
2009	WHERE(veget_4veg(:,igrass) > min_stomate)
2010	glccDef(:,g2c) = veget_4veg(:,igrass)-glcc(:,g2c)
2011	WHERE(veget_4veg(:,igrass)>glcc(:,g2c))
2012	glccReal(:,g2c) = glcc(:,g2c)
2013	ELSEWHERE
2014	glccReal(:,g2c) = veget_4veg(:,igrass)
2015	ENDWHERE
2016	ELSEWHERE
2017	glccReal(:,g2c) = 0.
2018	glccDef(:,g2c) = -1*glcc(:,g2c)
2019	ENDWHERE
2020
2021	!! 1. Compensation sequence: pasture,grass,forest
2022	tmpdef(:) = glccDef(:,f2c)+glccDef(:,g2c)+glccDef(:,p2c)
2023	WHERE(glccDef(:,p2c)<0)
2024	WHERE(glccDef(:,g2c)<0)
2025	WHERE(glccDef(:,f2c)<0) ! 1 (-,-,-)
2026	IncreDeficit(:,icrop) = tmpdef(:)
2027	ELSEWHERE ! 2 (-,-,+)
2028	WHERE(tmpdef(:)>=min_stomate)
2029	glccReal(:,f2c) = glccReal(:,f2c)-glccDef(:,g2c)-glccDef(:,p2c)
2030	ELSEWHERE
2031	glccReal(:,f2c) = veget_4veg(:,itree)
2032	IncreDeficit(:,icrop) = tmpdef(:)
2033	ENDWHERE
2034	ENDWHERE
2035	ELSEWHERE
2036	WHERE(glccDef(:,f2c)<0) ! 3 (-,+,-)
2037	WHERE(tmpdef(:)>=min_stomate)
2038	glccReal(:,g2c) = glccReal(:,g2c)-glccDef(:,p2c)-glccDef(:,f2c)
2039	ELSEWHERE
2040	glccReal(:,g2c) = veget_4veg(:,igrass)
2041	IncreDeficit(:,icrop) = tmpdef(:)
2042	ENDWHERE
2043	ELSEWHERE ! 4 (-,+,+)
2044	WHERE(tmpdef(:)>=min_stomate)
2045	WHERE((glccDef(:,g2c)+glccDef(:,p2c))>=min_stomate)
2046	glccReal(:,g2c) = glccReal(:,g2c)-glccDef(:,p2c)
2047	ELSEWHERE
2048	glccReal(:,g2c) = veget_4veg(:,igrass)
2049	glccReal(:,f2c) = glccReal(:,f2c)-(glccDef(:,p2c)+glccDef(:,g2c))
2050	ENDWHERE
2051	ELSEWHERE
2052	glccReal(:,g2c) = veget_4veg(:,igrass)
2053	glccReal(:,f2c) = veget_4veg(:,itree)
2054	IncreDeficit(:,icrop) = tmpdef(:)
2055	ENDWHERE
2056	ENDWHERE
2057	ENDWHERE
2058	ELSEWHERE
2059	WHERE(glccDef(:,g2c)<0)
2060	WHERE(glccDef(:,f2c)<0) ! 5 (+,-,-)
2061	WHERE(tmpdef(:)>=min_stomate)
2062	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,g2c)-glccDef(:,f2c)
2063	ELSEWHERE
2064	IncreDeficit(:,icrop) = tmpdef(:)
2065	glccReal(:,p2c) = veget_4veg(:,ipasture)
2066	ENDWHERE
2067	ELSEWHERE ! 6 (+,-,+)
2068	WHERE(tmpdef(:)>=min_stomate)
2069	WHERE((glccDef(:,p2c)+glccDef(:,g2c))>=min_stomate)
2070	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,g2c)
2071	ELSEWHERE
2072	glccReal(:,p2c) = veget_4veg(:,ipasture)
2073	glccReal(:,f2c) = glccReal(:,f2c)-(glccDef(:,g2c)+glccDef(:,p2c))
2074	ENDWHERE
2075	ELSEWHERE
2076	IncreDeficit(:,icrop) = tmpdef(:)
2077	glccReal(:,p2c) = veget_4veg(:,ipasture)
2078	glccReal(:,f2c) = veget_4veg(:,itree)
2079	ENDWHERE
2080	ENDWHERE
2081	ELSEWHERE
2082	WHERE(glccDef(:,f2c)<0) ! 7 (+,+,-)
2083	WHERE(tmpdef(:)>=min_stomate)
2084	WHERE((glccDef(:,p2c)+glccDef(:,f2c))>=min_stomate)
2085	glccReal(:,p2c) = glccReal(:,p2c)-glccDef(:,f2c)
2086	ELSEWHERE
2087	glccReal(:,p2c) = veget_4veg(:,ipasture)
2088	glccReal(:,g2c) = glccReal(:,g2c)-(glccDef(:,f2c)+glccDef(:,p2c))
2089	ENDWHERE
2090	ELSEWHERE
2091	IncreDeficit(:,icrop) = tmpdef(:)
2092	glccReal(:,g2c) = veget_4veg(:,igrass)
2093	glccReal(:,p2c) = veget_4veg(:,ipasture)
2094	ENDWHERE
2095	ELSEWHERE ! 8 (+,+,+)
2096	!do nothing
2097	ENDWHERE
2098	ENDWHERE
2099	ENDWHERE
2100	veget_4veg(:,itree) = veget_4veg(:,itree) - glccReal(:,f2c)
2101	veget_4veg(:,igrass) = veget_4veg(:,igrass) - glccReal(:,g2c)
2102	veget_4veg(:,ipasture) = veget_4veg(:,ipasture) - glccReal(:,p2c)
2103
2104	END SUBROUTINE glcc_compensation_full
2105
2106
2107
2108	!! This subroutine implements non-full compensation, is currently
2109	!! abandoned.
2110	SUBROUTINE glcc_compensation(npts,veget_4veg,glcc,glccDef, &
2111	p2c,ipasture,g2c,igrass,f2c,itree,icrop, &
2112	IncreDeficit)
2113
2114	IMPLICIT NONE
2115
2116	!! 0.1 Input variables
2117	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
2118	REAL(r_std), DIMENSION(npts,4), INTENT(in) :: veget_4veg !! "maximal" coverage fraction of a PFT on the ground
2119	INTEGER, INTENT(in) :: p2c,ipasture,g2c,igrass,f2c,itree,icrop
2120
2121	!! 0.2 Output variables
2122
2123
2124	!! 0.3 Modified variables
2125	REAL(r_std), DIMENSION (npts,12),INTENT(inout) :: glcc !! the land-cover-change (LCC) matrix in case a gross LCC is
2126	!! used.
2127	REAL(r_std), DIMENSION(npts,12), INTENT(inout) :: glccDef !! Gross LCC deficit, negative values mean that there
2128	!! are not enough fractions in the source vegetations
2129	!! to the target ones as presribed by the LCC matrix.
2130	REAL(r_std), DIMENSION(npts,4), INTENT(inout) :: IncreDeficit !! "Increment" deficits, negative values mean that
2131	!! there are not enough fractions in the source PFTs
2132	!! /vegetations to target PFTs/vegetations. I.e., these
2133	!! fraction transfers are presribed in LCC matrix but
2134	!! not realized.
2135
2136	!! 0.4 Local variables
2137	REAL(r_std), DIMENSION(npts) :: glccDef_all !! LCC deficits by summing up all the deficits to the
2138	!! the same target vegetation.
2139
2140
2141	WHERE(veget_4veg(:,itree) > min_stomate)
2142	glccDef(:,f2c) = veget_4veg(:,itree)-glcc(:,f2c)
2143	ELSEWHERE
2144	glccDef(:,f2c) = -1*glcc(:,f2c)
2145	glcc(:,f2c) = 0.
2146	ENDWHERE
2147
2148	WHERE(veget_4veg(:,ipasture) > min_stomate)
2149	glccDef(:,p2c) = veget_4veg(:,ipasture)-glcc(:,p2c)
2150	ELSEWHERE
2151	glccDef(:,p2c) = -1*glcc(:,p2c)
2152	glcc(:,p2c) = 0.
2153	ENDWHERE
2154
2155	WHERE(veget_4veg(:,igrass) > min_stomate)
2156	glccDef(:,g2c) = veget_4veg(:,igrass)-glcc(:,g2c)
2157	ELSEWHERE
2158	glccDef(:,g2c) = -1*glcc(:,g2c)
2159	glcc(:,g2c) = 0.
2160	ENDWHERE
2161
2162	glccDef_all(:) = glccDef(:,f2c)+glccDef(:,p2c)+glccDef(:,g2c)
2163
2164	! We allow the surpluses/deficits in p2c and g2c mutually compensating
2165	! for each other. If there are still deficits after this compensation,
2166	! they will be further compensated for by the surpluses from f2c (if there are any
2167	! surpluses). The ultimate deficits that cannot be compensated for
2168	! will be recorded and dropped.
2169
2170	! Because we assume the "pasture rule" is used, i.e., the crops
2171	! are supposed to come primarily from pastures and grasses, normally
2172	! we expect the deficits to occur in p2c or g2c rather than in f2c. But
2173	! if it happens that f2c has deficits while p2c or g2c has surpluse,
2174	! the surpluses will not be used to compensate for the f2c-deficits,
2175	! instead, we will just record and drop the f2c-deficits.
2176
2177	! In following codes for convenience we're not going to check
2178	! whether surpluses in f2c are enough to compensate for deficits
2179	! in p2c or g2c or both. Instead, we just add their deficits on top
2180	! of f2c. The issues of not-enough surpluses in f2c will be left for
2181	! the codes after this section to handle.
2182	WHERE (glccDef(:,p2c) < 0.)
2183	glcc(:,p2c) = veget_4veg(:,ipasture)
2184	WHERE (glccDef(:,g2c) < 0.)
2185	glcc(:,g2c) = veget_4veg(:,igrass)
2186	ELSEWHERE
2187	WHERE (glccDef(:,g2c)+glccDef(:,p2c) > min_stomate)
2188	glcc(:,g2c) = glcc(:,g2c)-glccDef(:,p2c)
2189	ELSEWHERE
2190	glcc(:,g2c) = veget_4veg(:,igrass)
2191	! whatever the case, we simply add the dificts to f2c
2192	glcc(:,f2c) = glcc(:,f2c)-glccDef(:,p2c)-glccDef(:,g2c)
2193	ENDWHERE
2194	ENDWHERE
2195
2196	ELSEWHERE
2197	WHERE(glccDef(:,g2c) < 0.)
2198	glcc(:,g2c) = veget_4veg(:,igrass)
2199	WHERE(glccDef(:,p2c)+glccDef(:,g2c) > min_stomate)
2200	glcc(:,p2c) = glcc(:,p2c)-glccDef(:,g2c)
2201	ELSEWHERE
2202	glcc(:,p2c) = veget_4veg(:,ipasture)
2203	! whatever the case, we simply add the dificts to f2c
2204	glcc(:,f2c) = glcc(:,f2c)-glccDef(:,p2c)-glccDef(:,g2c)
2205	ENDWHERE
2206	ELSEWHERE
2207	!Here p2c and g2c both show surplus, we're not going to check whether
2208	!glccDef(:,f2c) has negative values because we assume a "pasture rule"
2209	!is applied when constructing the gross LCC matrix, so deficits in
2210	!f2c will just be dropped but not be compensated for by the surpluses in
2211	!p2c or g2c.
2212	ENDWHERE
2213	ENDWHERE
2214
2215	! 1. We calculate again the f2c-deficit because f2c-glcc is adjusted in the
2216	! codes above as we allocated the deficits of p2c and g2c into f2c.
2217	! In cases where glccDef_all is less than zero, f2c-glcc will be larger
2218	! than available forest veget_max and we therefore limit the f2c-glcc to
2219	! available forest cover.
2220	! 2. There is (probably) a second case where glccDef_all is larger then zero,
2221	! but f2c-glcc is higher than veget_tree, i.e., Originally f2c is given a
2222	! high value that there is deficit in f2c but surpluses exist for p2c and g2c.
2223	! Normally we
2224	! assume this won't happen as explained above, given that a "pasture rule" was
2225	! used in constructing the gross LCC matrix. Nevertheless if this deos
2226	! happen, we will just drop the f2c deficit without being compensated
2227	! for by the surplus in p2c or g2c.
2228
2229	! we handle the 2nd case first
2230	WHERE(veget_4veg(:,itree) > min_stomate )
2231	WHERE(glccDef(:,f2c) < 0.)
2232	glcc(:,f2c) = veget_4veg(:,itree)
2233	WHERE (glccDef(:,p2c)+glccDef(:,g2c) > min_stomate)
2234	IncreDeficit(:,icrop) = glccDef(:,f2c)
2235	ELSEWHERE
2236	IncreDeficit(:,icrop) = glccDef_all(:)
2237	ENDWHERE
2238	ELSEWHERE
2239	WHERE(glccDef_all(:) < 0.) !handle the 1st case
2240	glcc(:,f2c) = veget_4veg(:,itree)
2241	IncreDeficit(:,icrop) = glccDef_all(:)
2242	ENDWHERE
2243	ENDWHERE
2244	ELSEWHERE
2245	WHERE(glccDef(:,p2c)+glccDef(:,g2c)>min_stomate)
2246	IncreDeficit(:,icrop) = glccDef(:,f2c)
2247	ELSEWHERE
2248	IncreDeficit(:,icrop) = glccDef_all(:)
2249	ENDWHERE
2250	ENDWHERE
2251
2252	END SUBROUTINE glcc_compensation
2253
2254	END MODULE stomate_gluc_common

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

Context Navigation

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

Download in other formats: