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

Last change on this file since 6737 was 5208, checked in by chao.yue, 6 years ago

Initialize several variables in each of glcc subroutines

File size: 48.6 KB
Line 
1! Remove definition of type_conversion
2
3! =================================================================================================================================
4! MODULE       : stomate_glcc_bioe1
5!
6! CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
7!
8! LICENCE      : IPSL (2006)
9! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10!
11!>\BRIEF       This module implements gross land use change with age classes.
12!!
13!!\n DESCRIPTION: None
14!!
15!! RECENT CHANGE(S): Including permafrost carbon
16!!
17!! REFERENCE(S) : None
18!!
19!! SVN          :
20!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/perso/albert.jornet/ORCHIDEE-MICT/src_stomate/stomate_lcchange.f90 $
21!! $Date: 2015-07-30 15:38:45 +0200 (Thu, 30 Jul 2015) $
22!! $Revision: 2847 $
23!! \n
24!_ ================================================================================================================================
25
26
27MODULE stomate_glcc_bioe1
28
29  ! modules used:
30 
31  USE ioipsl_para
32  USE stomate_data
33  USE pft_parameters
34  USE constantes
35  USE constantes_soil_var
36  USE stomate_gluc_common
37  USE stomate_gluc_constants
38  USE stomate_glcchange_MulAgeC
39
40  IMPLICIT NONE
41 
42  PRIVATE
43  PUBLIC glcc_bioe1_firstday, glcc_bioe1 
44 
45CONTAINS
46
47! ================================================================================================================================
48!! SUBROUTINE   gross_lcchange
49!!
50!>\BRIEF       : Apply gross land cover change.
51!!
52!>\DESCRIPTION 
53!_ ================================================================================================================================
54  SUBROUTINE glcc_bioe1 (npts, dt_days, newvegfrac,  &
55               glccSecondShift,&
56               def_fuel_1hr_remain, def_fuel_10hr_remain,        &
57               def_fuel_100hr_remain, def_fuel_1000hr_remain,    &
58               deforest_litter_remain, deforest_biomass_remain,  &
59               convflux, cflux_prod10, cflux_prod100,                  &
60               glcc_pft, glcc_pftmtc,          &
61               veget_max, prod10, prod100, flux10, flux100,            &
62               PFTpresent, senescence, moiavail_month, moiavail_week,  &
63               gpp_week, ngd_minus5, resp_maint, resp_growth,          & 
64               resp_hetero, npp_daily, when_growthinit, npp_longterm,  &
65               ind, lm_lastyearmax, everywhere, age,                   &
66               co2_to_bm, gpp_daily, co2_fire,                         &
67               time_hum_min, gdd_midwinter, gdd_from_growthinit,       &
68               gdd_m5_dormance, ncd_dormance,                          &
69               lignin_struc, carbon, leaf_frac,                        &
70               deepC_a, deepC_s, deepC_p,                              &
71               leaf_age, bm_to_litter, biomass, litter,                &
72               fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
73 
74    IMPLICIT NONE
75
76    !! 0.1 Input variables
77
78    INTEGER, INTENT(in)                                  :: npts             !! Domain size - number of pixels (unitless)
79    REAL(r_std), INTENT(in)                              :: dt_days          !! Time step of vegetation dynamics for stomate
80    REAL(r_std), DIMENSION (npts,12),INTENT(in)        :: glccSecondShift     !! the land-cover-change (LCC) matrix in case a gross LCC is
81                                                                              !! used.
82    REAL(r_std), DIMENSION (npts,nvmap),INTENT(in)       :: newvegfrac             !!
83                                                                             !!
84
85    REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in)                 :: def_fuel_1hr_remain
86    REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in)                 :: def_fuel_10hr_remain
87    REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in)                 :: def_fuel_100hr_remain
88    REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(in)                 :: def_fuel_1000hr_remain
89    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(in) :: deforest_litter_remain   !! Vegetmax-weighted remaining litter on the ground for
90                                                                                                      !! deforestation region.
91    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in)      :: deforest_biomass_remain  !! Vegetmax-weighted remaining biomass on the ground for
92                                                                                                      !! deforestation region.
93
94
95    !! 0.2 Output variables
96    REAL(r_std), DIMENSION(npts,nwp), INTENT(out)            :: convflux         !! release during first year following land cover
97                                                                             !! change
98    REAL(r_std), DIMENSION(npts,nwp), INTENT(out)            :: cflux_prod10     !! total annual release from the 10 year-turnover
99                                                                             !! pool @tex ($gC m^{-2}$) @endtex
100    REAL(r_std), DIMENSION(npts,nwp), INTENT(out)            :: cflux_prod100    !! total annual release from the 100 year-
101    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)      :: glcc_pft         !! Loss of fraction in each PFT
102    REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout):: glcc_pftmtc      !! a temporary variable to hold the fractions each PFT is going to lose
103                                                                             !! i.e., the contribution of each PFT to the youngest age-class of MTC
104
105    !! 0.3 Modified variables
106    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)      :: veget_max        !! "maximal" coverage fraction of a PFT (LAI ->
107                                                                             !! infinity) on ground (unitless)
108    REAL(r_std), DIMENSION(npts,0:10,nwp), INTENT(inout)     :: prod10           !! products remaining in the 10 year-turnover
109                                                                             !! pool after the annual release for each
110                                                                             !! compartment (10 + 1 : input from year of land
111                                                                             !! cover change)
112    REAL(r_std), DIMENSION(npts,0:100,nwp), INTENT(inout)    :: prod100          !! products remaining in the 100 year-turnover
113                                                                             !! pool after the annual release for each
114                                                                             !! compartment (100 + 1 : input from year of land
115                                                                             !! cover change)
116    REAL(r_std), DIMENSION(npts,10,nwp), INTENT(inout)       :: flux10           !! annual release from the 10/100 year-turnover
117                                                                             !! pool compartments
118    REAL(r_std), DIMENSION(npts,100,nwp), INTENT(inout)      :: flux100          !! annual release from the 10/100 year-turnover
119                                                                             !! pool compartments
120    LOGICAL, DIMENSION(:,:), INTENT(inout)               :: PFTpresent       !! Tab indicating which PFTs are present in
121                                                                             !! each pixel
122    LOGICAL, DIMENSION(:,:), INTENT(inout)               :: senescence       !! Flag for setting senescence stage (only
123                                                                             !! for deciduous trees)
124    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: moiavail_month   !! "Monthly" moisture availability (0 to 1,
125                                                                             !! unitless)
126    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: moiavail_week    !! "Weekly" moisture availability
127                                                                             !! (0 to 1, unitless)
128    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: gpp_week         !! Mean weekly gross primary productivity
129                                                                             !! @tex $(gC m^{-2} day^{-1})$ @endtex
130    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: ngd_minus5       !! Number of growing days (days), threshold
131                                                                             !! -5 deg C (for phenology)   
132    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: resp_maint       !! Maintenance respiration 
133                                                                             !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
134    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: resp_growth      !! Growth respiration 
135                                                                             !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
136    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: resp_hetero      !! Heterotrophic respiration 
137                                                                             !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
138    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: npp_daily        !! Net primary productivity
139                                                                             !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
140    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: when_growthinit  !! How many days ago was the beginning of
141                                                                             !! the growing season (days)
142    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: npp_longterm     !! "Long term" mean yearly primary productivity
143    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: ind              !! Number of individuals at the stand level
144                                                                             !! @tex $(m^{-2})$ @endtex
145    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: lm_lastyearmax   !! last year's maximum leaf mass for each PFT
146                                                                             !! @tex ($gC m^{-2}$) @endtex
147    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: everywhere       !! is the PFT everywhere in the grid box or
148                                                                             !! very localized (after its introduction) (?)
149    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: age              !! mean age (years)
150    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: co2_to_bm        !! CO2 taken from the atmosphere to get C to create 
151                                                                             !! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
152    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: gpp_daily        !! Daily gross primary productivity 
153                                                                             !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
154    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: co2_fire         !! Fire carbon emissions
155                                                                             !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
156    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: time_hum_min     !! Time elapsed since strongest moisture
157                                                                             !! availability (days)
158    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: gdd_midwinter    !! Growing degree days (K), since midwinter
159                                                                             !! (for phenology) - this is written to the
160                                                                             !!  history files
161    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: gdd_from_growthinit !! growing degree days, since growthinit
162                                                                             !! for crops
163    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: gdd_m5_dormance  !! Growing degree days (K), threshold -5 deg
164                                                                             !! C (for phenology)
165    REAL(r_std), DIMENSION(:,:), INTENT(inout)           :: ncd_dormance     !! Number of chilling days (days), since
166                                                                             !! leaves were lost (for phenology)
167    REAL(r_std), DIMENSION(:,:,:), INTENT(inout)         :: lignin_struc     !! ratio Lignine/Carbon in structural litter,
168                                                                             !! above and below ground
169    REAL(r_std), DIMENSION(:,:,:), INTENT(inout)         :: carbon           !! carbon pool: active, slow, or passive
170                                                                             !! @tex ($gC m^{-2}$) @endtex
171    REAL(r_std), DIMENSION(:,:,:), INTENT(inout)         :: deepC_a          !! Permafrost soil carbon (g/m**3) active
172    REAL(r_std), DIMENSION(:,:,:), INTENT(inout)         :: deepC_s          !! Permafrost soil carbon (g/m**3) slow
173    REAL(r_std), DIMENSION(:,:,:), INTENT(inout)         :: deepC_p          !! Permafrost soil carbon (g/m**3) passive
174    REAL(r_std), DIMENSION(:,:,:), INTENT(inout)         :: leaf_frac        !! fraction of leaves in leaf age class (unitless;0-1)
175    REAL(r_std), DIMENSION(:,:,:), INTENT(inout)         :: leaf_age         !! Leaf age (days)
176    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)       :: bm_to_litter     !! Transfer of biomass to litter
177                                                                             !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
178    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)       :: biomass          !! Stand level biomass @tex $(gC.m^{-2})$ @endtex
179    REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout)     :: litter           !! metabolic and structural litter, above and
180                                                                             !! below ground @tex ($gC m^{-2}$) @endtex
181    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)       :: fuel_1hr
182    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)       :: fuel_10hr
183    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)       :: fuel_100hr
184    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)       :: fuel_1000hr
185
186    !! 0.4 Local variables
187    REAL(r_std), DIMENSION(nvmap,nparts,nelements)             :: bm_to_litter_pro !! conversion of biomass to litter
188                                                                             !! @tex ($gC m^{-2} day^{-1}$) @endtex
189    REAL(r_std), DIMENSION(nvmap,nparts,nelements)             :: biomass_pro      !! biomass @tex ($gC m^{-2}$) @endtex
190    REAL(r_std), DIMENSION(nvmap)                        :: veget_max_pro    !! "maximal" coverage fraction of a PFT (LAI ->
191                                                                             !! infinity) on ground (unitless)
192    REAL(r_std), DIMENSION(nvmap,ncarb)                        :: carbon_pro       !! carbon pool: active, slow, or passive
193                                                                             !! @tex ($gC m^{-2}$) @endtex
194    REAL(r_std), DIMENSION(nvmap,ndeep)                        :: deepC_a_pro      !! Permafrost carbon pool: active, slow, or passive
195                                                                             !! @tex ($gC m^{-3}$) @endtex
196    REAL(r_std), DIMENSION(nvmap,ndeep)                        :: deepC_s_pro      !! Permafrost carbon pool: active, slow, or passive
197                                                                             !! @tex ($gC m^{-3}$) @endtex
198    REAL(r_std), DIMENSION(nvmap,ndeep)                        :: deepC_p_pro      !! Permafrost carbon pool: active, slow, or passive
199                                                                             !! @tex ($gC m^{-3}$) @endtex
200    REAL(r_std), DIMENSION(nvmap,nlitt,nlevs,nelements)        :: litter_pro       !! metabolic and structural litter, above and
201                                                                             !! below ground @tex ($gC m^{-2}$) @endtex
202    REAL(r_std), DIMENSION(nvmap,nlitt,nelements)              :: fuel_1hr_pro
203    REAL(r_std), DIMENSION(nvmap,nlitt,nelements)              :: fuel_10hr_pro
204    REAL(r_std), DIMENSION(nvmap,nlitt,nelements)              :: fuel_100hr_pro
205    REAL(r_std), DIMENSION(nvmap,nlitt,nelements)              :: fuel_1000hr_pro
206    REAL(r_std), DIMENSION(nvmap,nlevs)                        :: lignin_struc_pro !! ratio Lignine/Carbon in structural litter
207                                                                             !! above and below ground
208    REAL(r_std), DIMENSION(nvmap,nleafages)                    :: leaf_frac_pro    !! fraction of leaves in leaf age class
209    REAL(r_std), DIMENSION(nvmap,nleafages)                    :: leaf_age_pro     !! fraction of leaves in leaf age class
210    LOGICAL, DIMENSION(nvmap)                :: PFTpresent_pro, senescence_pro                 !! Is pft there (unitless)
211    REAL(r_std), DIMENSION(nvmap)            :: ind_pro, age_pro, lm_lastyearmax_pro, npp_longterm_pro
212    REAL(r_std), DIMENSION(nvmap)            :: everywhere_pro
213    REAL(r_std), DIMENSION(nvmap)            :: gpp_daily_pro, npp_daily_pro, co2_to_bm_pro
214    REAL(r_std), DIMENSION(nvmap)            :: resp_maint_pro, resp_growth_pro
215    REAL(r_std), DIMENSION(nvmap)            :: resp_hetero_pro, co2_fire_pro
216 
217    INTEGER                :: ipts,ivm,ivma,l,m,ipft_young_agec
218    CHARACTER(LEN=10)      :: part_str                               !! string suffix indicating an index
219
220    REAL(r_std), DIMENSION(npts,nvmap)       :: glcc_mtc             !! Increase in fraction of each MTC in its youngest age-class
221    REAL(r_std), DIMENSION(npts,nvm)         :: glccReal_tmp         !! A temporary variable
222
223    WRITE(numout,*) 'Entering glcc_MulAgeC'
224    glccReal_tmp(:,:) = zero
225
226    CALL glcc_bioe1_firstday(npts,veget_max,newvegfrac,   &
227                          glccSecondShift,glcc_pft,glcc_pftmtc)
228
229    glcc_mtc(:,:) = SUM(glcc_pftmtc,DIM=2)
230
231    DO ipts=1,npts
232
233      !! Initialize the _pro variables
234      bm_to_litter_pro(:,:,:)=zero                                               
235      biomass_pro(:,:,:)=zero
236      veget_max_pro(:)=zero                                                       
237      carbon_pro(:,:)=zero                                                       
238      deepC_a_pro(:,:)=zero                                                       
239      deepC_s_pro(:,:)=zero                                                       
240      deepC_p_pro(:,:)=zero                                                       
241      litter_pro(:,:,:,:)=zero                                                   
242      fuel_1hr_pro(:,:,:)=zero                                                   
243      fuel_10hr_pro(:,:,:)=zero                                                   
244      fuel_100hr_pro(:,:,:)=zero                                                 
245      fuel_1000hr_pro(:,:,:)=zero                                                 
246      lignin_struc_pro(:,:)=zero                                                 
247
248      leaf_frac_pro = zero
249      leaf_age_pro = zero
250      PFTpresent_pro(:) = .FALSE.
251      senescence_pro(:) = .TRUE.
252      ind_pro = zero
253      age_pro = zero
254      lm_lastyearmax_pro = zero
255      npp_longterm_pro = zero
256      everywhere_pro = zero
257     
258      gpp_daily_pro=zero                                                       
259      npp_daily_pro=zero                                                       
260      co2_to_bm_pro=zero                                                       
261      resp_maint_pro=zero                                                     
262      resp_growth_pro=zero                                                     
263      resp_hetero_pro=zero                                                     
264      co2_fire_pro=zero                                                       
265
266      ! Note that we assume people don't intentionally change baresoil to
267      ! vegetated land.
268      DO ivma = 2,nvmap
269
270        ! here we set glcc_mtc(ipts,ivma) > min_stomate as a condition,
271        ! this is necessary because later on in the subroutine of
272        ! add_incoming_proxy_pft we have to merge the newly established
273        ! youngest proxy with potentially exisiting youngest receiving MTC,
274        ! thus have to devide a new fraction of (frac_proxy + frac_exist),
275        ! but in case frac_exist = zero, we risk deviding by a very small value
276        ! of frac_proxy and thus we want it to be bigger than min_stomate.
277        IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
278
279          ! 1. we accumulate the scalar variables that will be inherited.
280          !    note that in the subroutine of `collect_legacy_pft`, all
281          !    zero transitions will be automatically skipped.
282          CALL collect_legacy_pft(npts, ipts, ivma, glcc_pftmtc,    &
283                  biomass, bm_to_litter, carbon, litter,            &
284                  deepC_a, deepC_s, deepC_p,                        &
285                  fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr,     &
286                  lignin_struc, co2_to_bm, gpp_daily, npp_daily,    &
287                  resp_maint, resp_growth, resp_hetero, co2_fire,   &
288                  def_fuel_1hr_remain, def_fuel_10hr_remain,        &
289                  def_fuel_100hr_remain, def_fuel_1000hr_remain,    &
290                  deforest_litter_remain, deforest_biomass_remain,  &
291                  veget_max_pro(ivma), carbon_pro(ivma,:),          &
292                  lignin_struc_pro(ivma,:), litter_pro(ivma,:,:,:), &
293                  deepC_a_pro(ivma,:), deepC_s_pro(ivma,:), deepC_p_pro(ivma,:), &
294                  fuel_1hr_pro(ivma,:,:), fuel_10hr_pro(ivma,:,:),  &
295                  fuel_100hr_pro(ivma,:,:), fuel_1000hr_pro(ivma,:,:), &
296                  bm_to_litter_pro(ivma,:,:), co2_to_bm_pro(ivma),  &
297                  gpp_daily_pro(ivma), npp_daily_pro(ivma),         &
298                  resp_maint_pro(ivma), resp_growth_pro(ivma),      &
299                  resp_hetero_pro(ivma), co2_fire_pro(ivma),        &
300                  convflux,prod10,prod100)
301
302          !++TEMP++
303          ! Here we substract the outgoing fraction from the source PFT.
304          ! If a too small fraction remains in this source PFT, then it is
305          ! exhausted, we empty it. The subroutine 'empty_pft' might be
306          ! combined with 'collect_legacy_pft', but now we just put it here.
307          DO ivm = 1,nvm
308            IF( glcc_pftmtc(ipts,ivm,ivma)>zero ) THEN
309              veget_max(ipts,ivm) = veget_max(ipts,ivm)-glcc_pftmtc(ipts,ivm,ivma)
310              IF ( veget_max(ipts,ivm)<min_stomate ) THEN
311                CALL empty_pft(ipts, ivm, veget_max, biomass, ind,       &
312                       carbon, litter, lignin_struc, bm_to_litter,       &
313                       deepC_a, deepC_s, deepC_p,                        &
314                       fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr,     &
315                       gpp_daily, npp_daily, gpp_week, npp_longterm,     &
316                       co2_to_bm, resp_maint, resp_growth, resp_hetero,  &
317                       lm_lastyearmax, leaf_frac, leaf_age, age,         &
318                       everywhere, PFTpresent, when_growthinit,          &
319                       senescence, gdd_from_growthinit, gdd_midwinter,   &
320                       time_hum_min, gdd_m5_dormance, ncd_dormance,      &
321                       moiavail_month, moiavail_week, ngd_minus5)
322              ENDIF
323            ENDIF
324          ENDDO
325
326        ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
327      ENDDO ! (DO ivma = 2,nvmap)
328
329      ! We can only establish new youngest proxy and add it to the
330      ! existing youngest-age PFT after all wood product extraction is done, to
331      ! avoid the dilution of extractable biomass by the young proxy
332      ! and ensure consistency. Therefore now we have to loop again
333      ! over nvmap.
334      DO ivma = 2,nvmap
335        !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate ) THEN
336
337          ipft_young_agec = start_index(ivma)
338
339          ! 2. we establish a proxy PFT with the fraction of veget_max_pro,
340          !    which is going to be either merged with existing target
341          !    `ipft_young_agec` PFT, or fill the place if no existing target PFT
342          !    exits.
343          CALL initialize_proxy_pft(ipts,ipft_young_agec,veget_max_pro(ivma), &
344                 biomass_pro(ivma,:,:), co2_to_bm_pro(ivma), ind_pro(ivma),   &
345                 age_pro(ivma),                                               & 
346                 senescence_pro(ivma), PFTpresent_pro(ivma),                  &
347                 lm_lastyearmax_pro(ivma), everywhere_pro(ivma),              &
348                 npp_longterm_pro(ivma),                                      &
349                 leaf_frac_pro(ivma,:),leaf_age_pro(ivma,:))
350
351          CALL sap_take (ipts,ivma,veget_max,biomass_pro(ivma,:,:), &
352                         biomass,co2_to_bm_pro(ivma))
353
354          ! 3. we merge the newly initiazlized proxy PFT into existing one
355          !    or use it to fill an empty PFT slot.
356          CALL add_incoming_proxy_pft(npts, ipts, ipft_young_agec, veget_max_pro(ivma),&
357                 carbon_pro(ivma,:), litter_pro(ivma,:,:,:), lignin_struc_pro(ivma,:), &
358                 bm_to_litter_pro(ivma,:,:),    &
359                 deepC_a_pro(ivma,:), deepC_s_pro(ivma,:), deepC_p_pro(ivma,:), &
360                 fuel_1hr_pro(ivma,:,:), fuel_10hr_pro(ivma,:,:),               &
361                 fuel_100hr_pro(ivma,:,:), fuel_1000hr_pro(ivma,:,:),           &
362                 biomass_pro(ivma,:,:), co2_to_bm_pro(ivma),                    &
363                 npp_longterm_pro(ivma), ind_pro(ivma),                         &
364                 lm_lastyearmax_pro(ivma), age_pro(ivma), everywhere_pro(ivma), & 
365                 leaf_frac_pro(ivma,:), leaf_age_pro(ivma,:),                   &
366                 PFTpresent_pro(ivma), senescence_pro(ivma),                &
367                 gpp_daily_pro(ivma), npp_daily_pro(ivma),                      &
368                 resp_maint_pro(ivma), resp_growth_pro(ivma),                   &
369                 resp_hetero_pro(ivma), co2_fire_pro(ivma),                     &
370                 veget_max, carbon, litter, lignin_struc, bm_to_litter,         &
371                 deepC_a, deepC_s, deepC_p,                                     &
372                 fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr,                  &
373                 biomass, co2_to_bm, npp_longterm, ind,                         &
374                 lm_lastyearmax, age, everywhere,                               &
375                 leaf_frac, leaf_age, PFTpresent, senescence,                   &
376                 gpp_daily, npp_daily, resp_maint, resp_growth,                 &
377                 resp_hetero, co2_fire)
378         
379        !ENDIF !IF ( glcc_mtc(ipts,ivma) .GT. min_stomate )
380      ENDDO !(DO ivma=1,nvmap)
381
382    ENDDO !(DO ipts=1,npts)
383
384    !! Update 10 year-turnover pool content following flux emission
385    !!     (linear decay (10%) of the initial carbon input)
386    DO  l = 0, 8
387      m = 10 - l
388      cflux_prod10(:,:) =  cflux_prod10(:,:) + flux10(:,m,:)
389      prod10(:,m,:)     =  prod10(:,m-1,:)   - flux10(:,m-1,:)
390      flux10(:,m,:)     =  flux10(:,m-1,:)
391    ENDDO
392   
393    cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,1,:) 
394    flux10(:,1,:)     = 0.1 * prod10(:,0,:)
395    prod10(:,1,:)     = prod10(:,0,:)
396   
397    !! 2.4.3 update 100 year-turnover pool content following flux emission\n
398    DO l = 0, 98
399       m = 100 - l
400       cflux_prod100(:,:)  =  cflux_prod100(:,:) + flux100(:,m,:)
401       prod100(:,m,:)      =  prod100(:,m-1,:)   - flux100(:,m-1,:)
402       flux100(:,m,:)      =  flux100(:,m-1,:)
403    ENDDO
404   
405    cflux_prod100(:,:)  = cflux_prod100(:,:) + flux100(:,1,:) 
406    flux100(:,1,:)      = 0.01 * prod100(:,0,:)
407    prod100(:,1,:)      = prod100(:,0,:)
408    prod10(:,0,:)        = zero
409    prod100(:,0,:)       = zero 
410
411    ! Write out history files
412    CALL histwrite_p (hist_id_stomate, 'glcc_pft', itime, &
413         glcc_pft, npts*nvm, horipft_index)
414
415    DO ivma = 1, nvmap
416      WRITE(part_str,'(I2)') ivma
417      IF (ivma < 10) part_str(1:1) = '0'
418      CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_'//part_str(1:LEN_TRIM(part_str)), &
419           itime, glcc_pftmtc(:,:,ivma), npts*nvm, horipft_index)
420    ENDDO
421
422  END SUBROUTINE glcc_bioe1
423
424
425! ================================================================================================================================
426!! SUBROUTINE   : glcc_bioe1_firstday
427!!
428!>\BRIEF        : When necessary, adjust input glcc matrix, and allocate it
429!!                into different contributing age classes and receiving
430!!                youngest age classes.
431!! \n
432!_ ================================================================================================================================
433
434  ! Note: it has this name because this subroutine will also be called
435  ! the first day of each year to precalculate the forest loss for the
436  ! deforestation fire module.
437  SUBROUTINE glcc_bioe1_firstday(npts,veget_max_org,newvegfrac, &
438                          glccSecondShift, glcc_pft, glcc_pftmtc)
439
440    IMPLICIT NONE
441
442    !! 0.1 Input variables
443
444    INTEGER, INTENT(in)                                     :: npts           !! Domain size - number of pixels (unitless)
445    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: veget_max_org  !! "maximal" coverage fraction of a PFT on the ground
446                                                                              !! May sum to
447                                                                              !! less than unity if the pixel has
448                                                                              !! nobio area. (unitless, 0-1)
449    REAL(r_std), DIMENSION(npts,nvmap), INTENT(in)          :: newvegfrac     !! used to guid the allocation of new PFTs.
450                                                                              !!
451    REAL(r_std), DIMENSION (npts,12),INTENT(in)        :: glccSecondShift     !! the land-cover-change (LCC) matrix in case a gross LCC is
452                                                                              !! used.
453
454    !! 0.2 Output variables
455    REAL(r_std), DIMENSION(npts,nvm,nvmap), INTENT(inout)   :: glcc_pftmtc    !! a temporary variable to hold the fractions each PFT is going to lose
456    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)         :: glcc_pft       !! Loss of fraction in each PFT
457
458    !! 0.3 Modified variables
459   
460    !! 0.4 Local variables
461    REAL(r_std), DIMENSION(npts,nvmap)              :: veget_mtc           !! "maximal" coverage fraction of a PFT on the ground
462    REAL(r_std), DIMENSION(npts,nvmap)              :: veget_mtc_begin     !! "maximal" coverage fraction of a PFT on the ground
463    REAL(r_std), DIMENSION(npts,nagec_tree)         :: vegagec_tree        !! fraction of tree age-class groups, in sequence of old->young
464    REAL(r_std), DIMENSION(npts,nagec_herb)         :: vegagec_grass       !! fraction of grass age-class groups, in sequence of old->young
465    REAL(r_std), DIMENSION(npts,nagec_herb)         :: vegagec_pasture     !! fraction of pasture age-class groups, in sequence of old->young
466    REAL(r_std), DIMENSION(npts,nagec_herb)         :: vegagec_crop        !! fraction of crop age-class groups, in sequence of old->young
467    REAL(r_std), DIMENSION(npts,nagec_bioe1)        :: vegagec_bioe1       !! fraction of crop age-class groups, in sequence of old->young
468
469   
470    REAL(r_std), DIMENSION(npts,4)                  :: veget_4veg      !! "maximal" coverage fraction of a PFT on the ground
471    REAL(r_std), DIMENSION(npts)                    :: veget_tree      !! "maximal" coverage fraction of a PFT on the ground
472    REAL(r_std), DIMENSION(npts)                    :: veget_grass     !! "maximal" coverage fraction of a PFT on the ground
473    REAL(r_std), DIMENSION(npts)                    :: veget_pasture   !! "maximal" coverage fraction of a PFT on the ground
474    REAL(r_std), DIMENSION(npts)                    :: veget_crop      !! "maximal" coverage fraction of a PFT on the ground
475    REAL(r_std), DIMENSION(npts)                    :: veget_bioe1     !! "maximal" coverage fraction of a PFT on the ground
476
477    REAL(r_std), DIMENSION(npts,nvm)         :: veget_max              !! "maximal" coverage fraction of a PFT on the ground
478    REAL(r_std), DIMENSION(npts,nvm)         :: veget_max_old          !! "maximal" coverage fraction of a PFT on the ground
479    REAL(r_std), DIMENSION(npts,nvm)         :: glcc_pft_tmp           !! Loss of fraction in each PFT
480
481    REAL(r_std), DIMENSION(npts,nvm,nvmap)  :: glcc_pftmtc_SecShift   !! a temporary variable to hold the fractions each PFT is going to lose
482    REAL(r_std), DIMENSION(npts,12)      :: glccRealSecShift   !! real matrix applied for secondary shifting cultivation.
483    REAL(r_std), DIMENSION(npts,12)      :: glccDefSecShift    !! deficit for the glccSecondShift
484    REAL(r_std), DIMENSION(npts,12)         :: glccRemain      !!
485    REAL(r_std), DIMENSION(npts,12)         :: glccSecondShift_remain      !!
486
487    REAL(r_std), DIMENSION(npts,nvm)         :: glccReal_tmp         !! A temporary variable
488
489    LOGICAL, SAVE  :: glcc_bioe1_firstday_done = .FALSE.
490
491    ! Different indexes for convenient local uses
492    INTEGER :: f_to_bioe1=1, g_to_bioe1=2, p_to_bioe1=3, c_to_bioe1=4, & 
493               bioe1_to_f=5, bioe1_to_g=6, bioe1_to_p=7, bioe1_to_c=8
494
495    INTEGER :: ivma
496
497    INTEGER :: ipts,IndStart_f,IndEnd_f
498    CHARACTER(LEN=10)      :: part_str                               !! string suffix indicating an index
499
500    !Some more local configurations
501    LOGICAL                                 :: allow_youngest_forest_convert = .TRUE.
502   
503    ! Initialization
504    glcc_pftmtc = zero
505    glcc_pft = zero
506    glcc_pft_tmp = zero
507
508    !!! ** Land cover change processes start here ** !!!
509    ! we make copies of original input veget_max (which is veget_max_org
510    ! in the subroutine parameter list).
511    ! veget_max will be modified through different operations in order to
512    ! check various purposes, e.g., whether input glcc matrix
513    ! is compatible with existing veget_max and how to allocate it etc.
514    ! veget_max_old will not be modified
515    veget_max(:,:) = veget_max_org(:,:)
516    veget_max_old(:,:) = veget_max_org(:,:)
517
518    !! 3. Treat secondary-agriculture shifting cultivation transition matrix
519    !! [The primary-agriculture shifting cultivation will be treated together
520    !!  with the netLCC transitions, with the conversion sequence of oldest->
521    !!  youngest is applied.]
522    ! When we prepare the driving data, secondary-agriculture shifting cultivation
523    ! is intended to include the "constant transitions" over time. Ideally, we
524    ! should start applying this secondary-agriculture shifting cultivation with
525    ! the "secondary forest" in the model. Here we tentatively start with the 3rd
526    ! youngest age class and move to the 2ne youngest age class. But if the prescribed
527    ! transition fraction is not met, we then move further to 4th youngest age class
528    ! and then move to the oldest age class sequentially.
529
530    CALL calc_cover_bioe1(npts,veget_max,veget_mtc_begin,vegagec_tree,vegagec_grass, &
531           vegagec_pasture,vegagec_crop,vegagec_bioe1)
532    veget_mtc = veget_mtc_begin
533 
534    !! 3.1 We start treating secondary-agriculture cultivation from the 3rd youngest
535    !! age class and then move to the younger age class.
536    ! Because it's rather complicated to calculate which transtion fraction between
537    ! which vegetation types should occur in case there is deficit occuring
538    ! for the overall donation vegetation type, we will just start from some
539    ! priority and leave the unrealized parts into the latter section.
540
541    ! For this purpose, we should first make a copy of glccSecondShift into
542    ! glccRemain. glccRemain will tell us the transition fractions that have to
543    ! be treated starting from `IndStart_f+1` oldest age class and moving torward older
544    ! age class.
545    glccRemain(:,:) = glccSecondShift(:,:)
546
547    ! Now we will call type_conversion for each of the 12 transitions, starting
548    ! from `IndStart_f` age class moving to the 2nd youngest age class. We use glccRemain
549    ! to track the transtion fractions we should leave for the second case.
550    ! To make the code more flexible, we will store the start and end indecies
551    ! in variables.
552
553    !*[Note: we do above process only for forest now, as we assume the conversion
554    !  of crop/pasture/grass to other types will start always from the oldest
555    !  age class]
556
557    IndStart_f = nagec_tree-1  ! note the indecies and vegetfrac for tree age class
558                               ! is from old to young
559    IndEnd_f = nagec_tree    ! nagec_tree-2: The 3rd youngest age class
560                               ! nagec_tree-1: The 2nd youngest age class
561                               ! nagec_tree: The youngest age class
562
563    IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
564      write(numout,*) 'glcc_bioe1_firstday: Age class index cannot be negative or zero!'
565      STOP
566    ENDIF
567
568    DO ipts=1,npts
569      !f_to_bioe1
570      CALL type_conversion(ipts,f_to_bioe1,glccSecondShift,veget_mtc,newvegfrac,       &
571                       indold_tree,indagec_tree,indagec_bioe1,num_bioe1_mulagec,     &
572                       IndEnd_f,nagec_bioe1,                    &
573                       vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
574                       glccRemain, &
575                       .TRUE., iagec_start=IndStart_f)
576      !g_to_bioe1
577      CALL type_conversion(ipts,g_to_bioe1,glccSecondShift,veget_mtc,newvegfrac,       &
578                       indold_grass,indagec_grass,indagec_bioe1,num_bioe1_mulagec,     &
579                       nagec_herb,nagec_bioe1,                    &
580                       vegagec_grass,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
581                       glccRemain, &
582                       .TRUE.)
583      !p_to_bioe1
584      CALL type_conversion(ipts,p_to_bioe1,glccSecondShift,veget_mtc,newvegfrac,       &
585                       indold_pasture,indagec_pasture,indagec_bioe1,num_bioe1_mulagec,     &
586                       nagec_herb,nagec_bioe1,                    &
587                       vegagec_pasture,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
588                       glccRemain, &
589                       .TRUE.)
590      !c_to_bioe1
591      CALL type_conversion(ipts,c_to_bioe1,glccSecondShift,veget_mtc,newvegfrac,       &
592                       indold_crop,indagec_crop,indagec_bioe1,num_bioe1_mulagec,     &
593                       nagec_herb,nagec_bioe1,                    &
594                       vegagec_crop,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
595                       glccRemain, &
596                       .TRUE.)
597      !bioe1_to_f
598      CALL type_conversion(ipts,bioe1_to_f,glccSecondShift,veget_mtc,newvegfrac,       &
599                       indold_bioe1,indagec_bioe1,indagec_tree,num_tree_mulagec,     &
600                       nagec_bioe1,nagec_tree,                    &
601                       vegagec_bioe1,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
602                       glccRemain, &
603                       .TRUE.)
604      !bioe1_to_g
605      CALL type_conversion(ipts,bioe1_to_g,glccSecondShift,veget_mtc,newvegfrac,       &
606                       indold_bioe1,indagec_bioe1,indagec_grass,num_grass_mulagec,     &
607                       nagec_bioe1,nagec_herb,                    &
608                       vegagec_bioe1,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
609                       glccRemain, &
610                       .TRUE.)
611      !bioe1_to_p
612      CALL type_conversion(ipts,bioe1_to_p,glccSecondShift,veget_mtc,newvegfrac,       &
613                       indold_bioe1,indagec_bioe1,indagec_pasture,num_pasture_mulagec,     &
614                       nagec_bioe1,nagec_herb,                    &
615                       vegagec_bioe1,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
616                       glccRemain, &
617                       .TRUE.)
618      !bioe1_to_c
619      CALL type_conversion(ipts,bioe1_to_c,glccSecondShift,veget_mtc,newvegfrac,       &
620                       indold_bioe1,indagec_bioe1,indagec_crop,num_crop_mulagec,     &
621                       nagec_bioe1,nagec_herb,                    &
622                       vegagec_bioe1,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp,&
623                       glccRemain, &
624                       .TRUE.)
625    ENDDO
626    glccSecondShift_remain(:,:) = glccRemain(:,:)
627
628    !! 3.2 We treat the remaing unrealized transtions from forest. Now we will
629    !! start with the 3rd oldest age class and then move to the oldest age class.
630
631    CALL calc_cover_bioe1(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
632           vegagec_pasture,vegagec_crop,vegagec_bioe1)
633    veget_mtc = veget_mtc_begin
634 
635    IndStart_f = nagec_tree  ! note the indecies and vegetfrac for tree age class
636                               ! is from old to young.
637                               ! nagec_tree -3: The 4th youngest age class.
638
639    IndEnd_f = 1               ! oldest-age class forest.
640
641    IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
642      write(numout,*) 'glcc_bioe1_firstday: Age class index cannot be negative or zero!'
643      STOP
644    ENDIF
645
646    ! we start with the 3rd youngest age class and move up to the oldest age
647    ! class in the sequence of young->old, as indicated by the .FALSE. parameter
648    ! when calling the subroutine type_conversion.
649    DO ipts=1,npts
650      !f_to_bioe1
651      CALL type_conversion(ipts,f_to_bioe1,glccSecondShift_remain,veget_mtc,newvegfrac,       &
652                       indold_tree,indagec_tree,indagec_bioe1,num_bioe1_mulagec,     &
653                       IndEnd_f,nagec_bioe1,                    &
654                       vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
655                       glccRemain, &
656                       .FALSE., iagec_start=IndStart_f)
657    ENDDO
658
659    IF (allow_youngest_forest_convert) THEN
660      !!++Temp++!!
661      !! this block of 3.3 could be commented to remove this process as desribed
662      !! below.
663
664      ! [2016-04-20] This is temporarily added: Normally we assume the youngest
665      ! forest age cohort will not be cut because in a shifting cultivation, they
666      ! are grown to let the land recover from agricultural process. (Or at least)
667      ! we can set the threshold of youngest age cohort to be very small. But there
668      ! are two reasons we allow the youngest forest cohort to be cut for shifting
669      ! cultivation purpose: a). Farmers may decide to harvest the wood of a forest
670      ! and then convert to crop. We don't simulate explicitly this process because
671      ! this will depend on input land change matrix and land use data assumptions.
672      ! However,we can implicitly account for this by assuming "farmers plant young
673      ! trees after harvesting the wood, and immediately convert this young trees
674      ! to crops. b). For the sake of conserving the total sum of veget_max before
675      ! and after the transition as one, we need to allow the youngest forest cohort
676      ! eligible for cutting.
677
678      !! 3.3 We treat the remaing unrealized transtions from forest, allowing
679      !! the youngest forest cohort to be cut. For this purpose, we will
680      !! start with the 2nd youngest age class and then move to the youngest one.
681
682      glccSecondShift_remain(:,:) = glccRemain(:,:)
683
684      CALL calc_cover_bioe1(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
685           vegagec_pasture,vegagec_crop,vegagec_bioe1)
686      veget_mtc = veget_mtc_begin
687 
688      ! Note: the setting of index here must be consistent with those of 3.1 and 3.2
689      IndStart_f = nagec_tree-1  ! note the indecies and vegetfrac for tree age class
690                                 ! is from old to young.
691                                 ! nagec_tree -1: The 2nd youngest age class.
692
693      IndEnd_f = nagec_tree      ! youngest class forest.
694
695      IF (IndStart_f .LE. 0 .OR. IndEnd_f .LE. 0) THEN
696        write(numout,*) 'glcc_bioe1_firstday: Age class index cannot be negative or zero!'
697        STOP
698      ENDIF
699
700      ! we start with the 3rd youngest age class and move up to the oldest age
701      ! class in the sequence of young->old, as indicated by the .FALSE. parameter
702      ! when calling the subroutine type_conversion.
703      DO ipts=1,npts
704        !f_to_bioe1
705        CALL type_conversion(ipts,f_to_bioe1,glccSecondShift_remain,veget_mtc,newvegfrac,       &
706                         indold_tree,indagec_tree,indagec_bioe1,num_bioe1_mulagec,     &
707                         IndEnd_f,nagec_bioe1,                    &
708                         vegagec_tree,veget_max,glcc_pft,glcc_pftmtc,glcc_pft_tmp, &
709                         glccRemain, &
710                         .FALSE., iagec_start=IndStart_f)
711      ENDDO
712      !! End of ++Temp++ Section 3.3
713    ENDIF
714
715    ! Final handling of some output variables.
716    ! we separate the glcc_pftmtc_SecShift
717    glcc_pftmtc_SecShift = glcc_pftmtc
718
719    ! we put the remaining glccRemain into the deficit
720    glccDefSecShift = -1 * glccRemain
721    glccRealSecShift = glccSecondShift - glccRemain
722
723    !*****end block to handle glcc involving bioenergy vegtation type *******
724
725    IF (.NOT. glcc_bioe1_firstday_done) THEN
726
727      glccReal_tmp = zero
728
729      glccReal_tmp(:,1:12) = glccRealSecShift
730      CALL histwrite_p (hist_id_stomate, 'glccRealSecShift', itime, &
731           glccReal_tmp, npts*nvm, horipft_index)
732
733      glccReal_tmp(:,1:12) = glccDefSecShift
734      CALL histwrite_p (hist_id_stomate, 'glccDefSecShift', itime, &
735           glccReal_tmp, npts*nvm, horipft_index)
736
737      DO ivma = 1, nvmap
738        WRITE(part_str,'(I2)') ivma
739        IF (ivma < 10) part_str(1:1) = '0'
740        CALL histwrite_p (hist_id_stomate, 'glcc_pftmtc_SF_'//part_str(1:LEN_TRIM(part_str)), &
741             itime, glcc_pftmtc_SecShift(:,:,ivma), npts*nvm, horipft_index)
742      ENDDO
743
744      glcc_bioe1_firstday_done = .TRUE.
745    ENDIF
746
747  END SUBROUTINE glcc_bioe1_firstday
748
749  SUBROUTINE calc_cover_bioe1(npts,veget_max,veget_mtc,vegagec_tree,vegagec_grass, &
750                 vegagec_pasture,vegagec_crop,vegagec_bioe1)
751
752   
753    IMPLICIT NONE
754
755    !! Input variables
756    INTEGER, INTENT(in)                                       :: npts             !! Domain size - number of pixels (unitless)
757    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)         :: veget_max           !! "maximal" coverage fraction of a PFT on the ground
758
759    !! Output variables
760    REAL(r_std), DIMENSION(npts,nvmap), INTENT(inout)         :: veget_mtc        !! "maximal" coverage fraction of a PFT on the ground
761    REAL(r_std), DIMENSION(npts,nagec_tree), INTENT(inout)    :: vegagec_tree     !! fraction of tree age-class groups, in sequence of old->young
762    REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout)    :: vegagec_grass    !! fraction of grass age-class groups, in sequence of old->young
763    REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout)    :: vegagec_pasture  !! fraction of pasture age-class groups, in sequence of old->young
764    REAL(r_std), DIMENSION(npts,nagec_herb), INTENT(inout)    :: vegagec_crop     !! fraction of crop age-class groups, in sequence of old->young
765    REAL(r_std), DIMENSION(npts,nagec_bioe1), INTENT(inout)    :: vegagec_bioe1     !! fraction of bioenergy tree age-class groups, in sequence of old->young
766
767    !! Local variables
768    INTEGER(i_std)                                          :: ivma,staind,endind,j    !! indices (unitless)
769
770    veget_mtc(:,:) = 0.
771    vegagec_tree(:,:) = 0.
772    vegagec_grass(:,:) = 0.
773    vegagec_pasture(:,:) = 0.
774    vegagec_crop(:,:) = 0.
775    vegagec_bioe1(:,:) = 0.
776
777    ! Calculate veget_max for MTCs
778    DO ivma = 1,nvmap
779      staind = start_index(ivma)
780      IF (nagec_pft(ivma) == 1) THEN
781        veget_mtc(:,ivma) = veget_max(:,staind)
782      ELSE
783        veget_mtc(:,ivma) = \
784          SUM(veget_max(:,staind:staind+nagec_pft(ivma)-1),DIM=2)
785      ENDIF
786    ENDDO
787
788    ! Calculate veget_max for each age class
789    DO ivma = 2,nvmap  !here we start with 2 to exclude baresoil (always PFT1)
790      staind = start_index(ivma)
791      endind = staind+nagec_pft(ivma)-1
792
793      ! Single-age-class MTC goest to oldest age class.
794      IF (nagec_pft(ivma) == 1) THEN
795        WRITE(numout,*) "Error: metaclass has only a single age group: ",ivma
796        STOP
797      ELSE
798        IF (is_tree(staind)) THEN
799          DO j=1,nagec_tree
800            vegagec_tree(:,j) = vegagec_tree(:,j)+veget_max(:,endind-j+1)
801          ENDDO
802        ELSE IF (is_bioe1(staind)) THEN
803          DO j=1,nagec_bioe1
804            vegagec_bioe1(:,j) = vegagec_bioe1(:,j)+veget_max(:,endind-j+1)
805          ENDDO
806        ELSE IF (is_grassland_manag(staind)) THEN
807          DO j=1,nagec_herb
808            vegagec_pasture(:,j) = vegagec_pasture(:,j)+veget_max(:,endind-j+1)
809          ENDDO
810        ELSE IF (natural(staind)) THEN
811          DO j=1,nagec_herb
812            vegagec_grass(:,j) = vegagec_grass(:,j)+veget_max(:,endind-j+1)
813          ENDDO
814        ELSE
815          DO j=1,nagec_herb
816            vegagec_crop(:,j) = vegagec_crop(:,j)+veget_max(:,endind-j+1)
817          ENDDO
818        ENDIF
819      ENDIF
820    ENDDO
821
822  END SUBROUTINE calc_cover_bioe1
823
824END MODULE stomate_glcc_bioe1
Note: See TracBrowser for help on using the repository browser.