source: branches/publications/ORCHIDEE_gmd-2018-182/src_stomate/stomate_lpj.f90

Last change on this file was 1170, checked in by josefine.ghattas, 11 years ago

Merge revision 1137 from branches/MERGE-OCN into the trunk.

  • Added supplementary dimension for carbon to following variables : bm_sapl, litter, bm_to_litter, turnover_daily, turnover_littercalc, turnover_longterm, bm_to_litttercalc and biomass. No change in results.

by Didier

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 69.4 KB
Line 
1! ================================================================================================================================
2! MODULE       : stomate_lpj
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       Main entry point for daily processes in STOMATE and LPJ (phenology,
10!! allocation, npp_calc, kill, turn, light, establish, crown, cover, lcchange)
11!!
12!!\n DESCRIPTION: None
13!!
14!! RECENT CHANGE(S): None
15!!
16!! REFERENCE(S) : None
17!!
18!! SVN          :
19!! $HeadURL$
20!! $Date$
21!! $Revision$
22!! \n
23!_ ================================================================================================================================
24
25MODULE stomate_lpj
26
27  ! modules used:
28
29  USE ioipsl_para
30  USE grid
31  USE stomate_data
32  USE constantes
33  USE constantes_soil
34  USE pft_parameters
35  USE lpj_constraints
36  USE lpj_pftinout
37  USE lpj_kill
38  USE lpj_crown
39  USE lpj_fire
40  USE lpj_gap
41  USE lpj_light
42  USE lpj_establish
43  USE lpj_cover
44  USE stomate_prescribe
45  USE stomate_phenology
46  USE stomate_alloc
47  USE stomate_npp
48  USE stomate_turnover
49  USE stomate_litter
50  USE stomate_soilcarbon
51  USE stomate_vmax
52  USE stomate_assimtemp
53  USE stomate_lcchange
54!  USE Orch_Write_field_p
55
56
57  IMPLICIT NONE
58
59  ! private & public routines
60
61  PRIVATE
62  PUBLIC StomateLpj,StomateLpj_clear
63
64  LOGICAL, SAVE                         :: firstcall = .TRUE.             !! first call
65!$OMP THREADPRIVATE(firstcall)
66
67CONTAINS
68
69
70!! ================================================================================================================================
71!! SUBROUTINE   : StomateLpj_clear
72!!
73!>\BRIEF        Re-initialisation of variable
74!!
75!! DESCRIPTION  : This subroutine reinitializes variables. To be used if we want to relaunch
76!! ORCHIDEE but the routine is not used in current version.
77!!
78!! RECENT CHANGE(S) : None
79!!
80!! MAIN OUTPUT VARIABLE(S): None
81!!
82!! REFERENCE(S) : None
83!!
84!! FLOWCHART    : None
85!! \n
86!_ ================================================================================================================================
87
88  SUBROUTINE StomateLpj_clear
89
90    CALL prescribe_clear
91    CALL phenology_clear
92    CALL npp_calc_clear
93    CALL turn_clear
94    CALL soilcarbon_clear
95    CALL constraints_clear
96    CALL establish_clear
97    CALL fire_clear
98    CALL gap_clear
99    CALL light_clear
100    CALL pftinout_clear
101    CALL alloc_clear
102  END SUBROUTINE StomateLpj_clear
103
104
105!! ================================================================================================================================
106!! SUBROUTINE   : StomateLPJ
107!!
108!>\BRIEF        Main entry point for daily processes in STOMATE and LPJ, structures the call sequence
109!!              to the different processes such as dispersion, establishment, competition and mortality of PFT's.
110!!
111!! DESCRIPTION  : This routine is the main entry point to all processes calculated on a
112!! daily time step. Is mainly devoted to call the different STOMATE and LPJ routines
113!! depending of the ok_dgvm (is dynamic veg used) and lpj_constant_mortality (is background mortality used).
114!! It also prepares the cumulative
115!! fluxes or pools (e.g TOTAL_M TOTAL_BM_LITTER etc...)
116!!
117!! This routine makes frequent use of "weekly", "monthly" and "long term" variables. Quotion is used because
118!! by default "weekly" denotes 7 days, by default "monthly" denotes 20 days and by default "Long term" denotes
119!! 3 years. dtslow refers to 24 hours (1 day).
120!!
121!!
122!! RECENT CHANGE(S) : None
123!!
124!! MAIN OUTPUT VARIABLE(S): All variables related to stomate and required for LPJ dynamic vegetation mode.
125!!
126!! REFERENCE(S) :
127!! - Krinner, G., N. Viovy, N. de Noblet-Ducoudré, J. Ogeé, J. Polcher, P. Friedlingstein, P. Ciais, S. Sitch,
128!! and I. C. Prentice. 2005. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere
129!! system. Global Biogeochemical Cycles 19:GB1015, doi:1010.1029/2003GB002199.
130!! - Sitch, S., B. Smith, I. C. Prentice, A. Arneth, A. Bondeau, W. Cramer, J. O. Kaplan, S. Levis, W. Lucht,
131!! M. T. Sykes, K. Thonicke, and S. Venevsky. 2003. Evaluation of ecosystem dynamics, plant geography and
132!! terrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biology 9:161-185.
133!!
134!! FLOWCHART    : Update with existing flowchart from N Viovy (Jan 19, 2012)
135!! \n
136!_ ================================================================================================================================
137 
138  SUBROUTINE StomateLpj (npts, dt_days, EndOfYear, EndOfMonth, &
139       neighbours, resolution, &
140       clay, herbivores, &
141       tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
142       litterhum_daily, soilhum_daily, &
143       maxmoiavail_lastyear, minmoiavail_lastyear, &
144       gdd0_lastyear, precip_lastyear, &
145       moiavail_month, moiavail_week, tlong_ref, t2m_month, t2m_week, &
146       tsoil_month, soilhum_month, &
147       gdd_m5_dormance, gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
148       turnover_longterm, gpp_daily, &
149       time_hum_min, maxfpc_lastyear, resp_maint_part, &
150       PFTpresent, age, fireindex, firelitter, &
151       leaf_age, leaf_frac, biomass, ind, adapted, regenerate, &
152       senescence, when_growthinit, &
153       litterpart, litter, dead_leaves, carbon, black_carbon, lignin_struc, &
154       veget_max, npp_longterm, lm_lastyearmax, veget_lastlight, &
155       everywhere, need_adjacent, RIP_time, &
156       lai, rprof,npp_daily, turnover_daily, turnover_time,&
157       control_moist, control_temp, soilcarbon_input, &
158       co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, &
159       height, deadleaf_cover, vcmax, vjmax, &
160       t_photo_min, t_photo_opt, t_photo_max,bm_to_litter, &
161       prod10,prod100,flux10, flux100, veget_max_new, &
162       convflux,cflux_prod10,cflux_prod100, harvest_above, carb_mass_total, lcchange, &
163       fpc_max, MatrixA)
164   
165  !! 0. Variable and parameter declaration
166
167    !! 0.1 input
168
169    INTEGER(i_std), INTENT(in)                                 :: npts                 !! Domain size (unitless)
170    REAL(r_std), INTENT(in)                                    :: dt_days              !! Time step of Stomate (days)
171    INTEGER(i_std), DIMENSION(npts,8), INTENT(in)              :: neighbours           !! Indices of the 8 neighbours of each grid
172                                                                                       !! point [1=N, 2=NE, 3=E, 4=SE,
173                                                                                       !!  5=S, 6=SW, 7=W, 8=NW]
174    REAL(r_std), DIMENSION(npts,2), INTENT(in)                 :: resolution           !! Resolution at each grid point (m) 
175                                                                                       !! [1=E-W, 2=N-S]
176    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: clay                 !! Clay fraction (0 to 1, unitless)
177    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: herbivores           !! Time constant of probability of a leaf to
178                                                                                       !! be eaten by a herbivore (days)
179    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: tsurf_daily          !! Daily surface temperatures (K)
180    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: tsoil_daily          !! Daily soil temperatures (K)
181    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_daily            !! Daily 2 meter temperatures (K)
182    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_min_daily        !! Daily minimum 2 meter temperatures (K)
183    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: litterhum_daily      !! Daily litter humidity (0 to 1, unitless)
184    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: soilhum_daily        !! Daily soil humidity (0 to 1, unitless)
185    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: maxmoiavail_lastyear !! Last year's maximum moisture availability
186                                                                                       !! (0 to 1, unitless)
187    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: minmoiavail_lastyear !! Last year's minimum moisture availability
188                                                                                       !! (0 to 1, unitless)
189    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: gdd0_lastyear        !! Last year's GDD0 (K)
190    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: precip_lastyear      !! Lastyear's precipitation
191                                                                                       !! @tex $(mm year^{-1})$ @endtex
192                                                                                       !! to determine if establishment possible
193    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: moiavail_month       !! "Monthly" moisture availability (0 to 1,
194                                                                                       !! unitless)
195    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: moiavail_week        !! "Weekly" moisture availability
196                                                                                       !! (0 to 1, unitless)
197    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: tlong_ref            !! "Long term" 2 meter reference
198                                                                                       !! temperatures (K)
199    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_month            !! "Monthly" 2-meter temperatures (K)
200    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_week             !! "Weekly" 2-meter temperatures (K)
201    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: tsoil_month          !! "Monthly" soil temperatures (K)
202    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: soilhum_month        !! "Monthly" soil humidity
203                                                                                       !! (0 to 1, unitless)
204    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gdd_m5_dormance      !! Growing degree days (K), threshold -5 deg
205                                                                                       !! C (for phenology)
206    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gdd_from_growthinit  !! growing degree days, since growthinit for crops
207    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: gdd_midwinter        !! Growing degree days (K), since midwinter
208                                                                                       !! (for phenology) - this is written to the history files
209    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: ncd_dormance         !! Number of chilling days (days), since
210                                                                                       !! leaves were lost (for phenology)
211    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: ngd_minus5           !! Number of growing days (days), threshold
212                                                                                       !! -5 deg C (for phenology)
213    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)        :: turnover_longterm    !! "Long term" turnover rate 
214                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
215    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gpp_daily            !! Daily gross primary productivity 
216                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
217    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: time_hum_min         !! Time elapsed since strongest moisture
218                                                                                       !! availability (days)
219    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: maxfpc_lastyear      !! Last year's maximum foliage projected
220                                                                                       !! coverage for each natural PFT,
221                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
222    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)        :: resp_maint_part      !! Maintenance respiration of different
223                                                                                       !! plant parts 
224                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
225    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: fpc_max              !! "Maximal" coverage fraction of a PFT (LAI
226                                                                                       !! -> infinity) on ground 
227                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
228    LOGICAL, INTENT(in)                                        :: lcchange             !! Land cover change flag
229    LOGICAL, INTENT(in)                                        :: EndOfYear            !! Flag set on the last day of the year used
230                                                                                       !! to update "yearly variables". This
231                                                                                       !! variable must be .TRUE. once a year
232    LOGICAL, INTENT(in)                                        :: EndOfMonth           !! Flag set at end of each month to update
233                                                                                       !! monthly variable
234
235  !! 0.2 Output variables
236   
237    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: npp_daily            !! Net primary productivity
238                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
239    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out) :: turnover_daily       !! Turnover rates
240                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
241    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: co2_to_bm            !! CO2 taken up from atmosphere when
242                                                                                       !! introducing a new PFT (introduced for
243                                                                                       !! carbon balance closure)
244                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
245    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: co2_fire             !! Carbon emitted into the atmosphere by
246                                                                                       !! fire (living and dead biomass) 
247                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
248    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: resp_hetero          !! Heterotrophic respiration
249                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
250    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: resp_maint           !! Maintenance respiration 
251                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
252    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: resp_growth          !! Growth respiration 
253                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
254   
255    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: deadleaf_cover       !! Fraction of soil covered by dead leaves
256                                                                                       !! (0 to 1, unitless)
257    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: vcmax                !! Maximum rate of carboxylation
258    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: vjmax                !! Maximum rate of RUbp regeneration 
259                                                                                       !! @tex $(\mumol CO_2 m^{-2}s^{-1})$ @endtex
260    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: t_photo_min          !! Minimum temperature for photosynthesis
261                                                                                       !! (C) This is not K because it is a scalar.
262    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: t_photo_opt          !! Optimum temperature for photosynthesis
263                                                                                       !! (C) This is not K because it is a scalar.
264    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: t_photo_max          !! Maximum temperature for photosynthesis
265                                                                                       !! (C) This is not K because it is a scalar.
266    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out):: bm_to_litter      !! Conversion of biomass to litter
267                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
268
269    !! 0.3 Modified variables
270   
271    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: height               !! Height of vegetation (m)
272    REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout)          :: control_moist        !! Moisture control of heterotrophic
273                                                                                       !! respiration (0 to 1, unitless)
274    REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout)          :: control_temp         !! Temperature control of heterotrophic
275                                                                                       !! respiration, above and below
276                                                                                       !! (0 to 1, unitless)
277    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout)      :: soilcarbon_input     !! Quantity of carbon going into carbon
278                                                                                       !! pools from litter decomposition 
279                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
280    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: lai                  !! Leaf area index OF AN INDIVIDUAL PLANT,
281                                                                                       !! where a PFT contains n indentical plants
282                                                                                       !! i.e., using the mean individual approach
283                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
284    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: rprof                !! Prescribed root depth (m)
285    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: PFTpresent           !! Tab indicating which PFTs are present in
286                                                                                       !! each pixel
287    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: age                  !! Age (years)   
288    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: fireindex            !! Probability of fire (0 to 1, unitless)
289    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: firelitter           !! Longer term litter above the ground that
290                                                                                       !! can be burned, @tex $(gC m^{-2})$ @endtex
291    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_age             !! Leaf age (days)
292    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_frac            !! Fraction of leaves in leaf age class,
293                                                                                       !! (0 to 1, unitless)
294    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass        !! Biomass @tex $(gC m^{-2})$ @endtex
295    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: ind                  !! Density of individuals
296                                                                                       !! @tex $(m^{-2})$ @endtex
297    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: adapted              !! Adaptation of PFT (killed if too cold)
298                                                                                       !! (0 to 1, unitless)
299    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: regenerate           !! "Fitness": Winter sufficiently cold for
300                                                                                       !! PFT regeneration ? (0 to 1, unitless)
301    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: senescence           !! Flag for setting senescence stage (only
302                                                                                       !! for deciduous trees)
303    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: when_growthinit      !! How many days ago was the beginning of
304                                                                                       !! the growing season (days)
305    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)      :: litterpart           !! Fraction of litter above the ground
306                                                                                       !! belonging to different PFTs
307                                                                                       !! (0 to 1, unitless)
308    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout):: litter     !! Metabolic and structural litter, above
309                                                                                       !! and below ground
310                                                                                       !! @tex $(gC m^{-2})$ @endtex
311    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)      :: dead_leaves          !! Dead leaves on ground, per PFT, metabolic
312                                                                                       !! and structural, 
313                                                                                       !! @tex $(gC m^{-2})$ @endtex
314    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout)      :: carbon               !! Carbon pool: active, slow, or passive,
315                                                                                       !! @tex $(gC m^{-2})$ @endtex 
316    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: black_carbon         !! Black carbon on the ground
317                                                                                       !! @tex $(gC m^{-2})$ @endtex
318    REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout)      :: lignin_struc         !! Ratio of Lignin/Carbon in structural
319                                                                                       !! litter, above and below ground, 
320                                                                                       !! @tex $(gC m^{-2})$ @endtex
321    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: veget_max            !! "Maximal" coverage fraction of a PFT (LAI
322                                                                                       !! -> infinity) on ground
323    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: npp_longterm         !! "Long term" mean yearly primary
324                                                                                       !! productivity
325                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
326    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: lm_lastyearmax       !! Last year's maximum leaf mass, for each
327                                                                                       !! PFT @tex $(gC m^{-2})$ @endtex
328    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: veget_lastlight      !! Vegetation fractions (on ground) after
329                                                                                       !! last light competition 
330                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
331    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: everywhere           !! Is the PFT everywhere in the grid box or
332                                                                                       !! very localized (after its introduction)
333                                                                                       !! (unitless)
334    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: need_adjacent        !! In order for this PFT to be introduced,
335                                                                                       !! does it have to be present in an
336                                                                                       !! adjacent grid box?
337    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: RIP_time             !! How much time ago was the PFT eliminated
338                                                                                       !! for the last time (y)
339    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: turnover_time        !! Turnover_time of leaves for grasses
340                                                                                       !! (days)
341    REAL(r_std), DIMENSION(npts,nvm),INTENT(inout)             :: veget_max_new        !! New "maximal" coverage fraction of a PFT
342                                                                                       !! (LAI -> infinity) (unitless)
343    REAL(r_std),DIMENSION(npts,0:10), INTENT(inout)            :: prod10               !! Products remaining in the 10
344                                                                                       !! year-turnover pool after the annual
345                                                                                       !! release for each compartment (10
346                                                                                       !! + 1 : input from year of land cover
347                                                                                       !! change) @tex $(gC m^{-2})$ @endtex
348    REAL(r_std),DIMENSION(npts,0:100), INTENT(inout)           :: prod100              !! Products remaining in the 100
349                                                                                       !! year-turnover pool after the annual
350                                                                                       !! release for each compartment (100
351                                                                                       !! + 1 : input from year of land cover
352                                                                                       !! change) @tex $(gC m^{-2})$ @endtex
353    REAL(r_std),DIMENSION(npts,10), INTENT(inout)              :: flux10               !! Annual release from the 10
354                                                                                       !! year-turnover pool compartments 
355                                                                                       !! @tex $(gC m^{-2})$ @endtex
356    REAL(r_std),DIMENSION(npts,100), INTENT(inout)             :: flux100              !! Annual release from the 100
357                                                                                       !! year-turnover pool compartments 
358                                                                                       !! @tex $(gC m^{-2})$ @endtex
359    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: convflux             !! Release during first year following land
360                                                                                       !! cover change @tex $(gC m^{-2})$ @endtex
361    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: cflux_prod10         !! Total annual release from the 10
362                                                                                       !! year-turnover pool
363                                                                                       !! @tex $(gC m^{-2})$ @endtex
364    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: cflux_prod100        !! Total annual release from the 100
365                                                                                       !! year-turnover pool
366                                                                                       !! @tex $(gC m^{-2})$ @endtex
367    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: harvest_above        !! Harvest above ground biomass for
368                                                                                       !! agriculture @tex $(gC m^{-2})$ @endtex
369    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: carb_mass_total      !! Carbon Mass total (soil, litter, veg)
370                                                                                       !! @tex $(gC m^{-2})$ @endtex 
371    REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA         !! Matrix containing the fluxes 
372                                                                                       !! between the carbon pools
373                                                                                       !! per sechiba time step
374                                                                                       !! @tex $(gC.m^2.day^{-1})$ @endtex
375
376    !! 0.4 Local variables
377
378    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_bm_to_litter    !! Total conversion of biomass to litter
379                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
380    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_live_biomass    !! Total living biomass 
381                                                                                       !! @tex $(gC m{-2})$ @endtex
382    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements)           :: bm_alloc            !! Biomass increase, i.e. NPP per plant part
383                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
384    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_turnover        !! Total turnover rate 
385                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
386    REAL(r_std), DIMENSION(npts,nvm)                            :: tot_litter_soil_carb!! Total soil and litter carbon 
387                                                                                       !! @tex $(gC m^{-2})$ @endtex
388    REAL(r_std), DIMENSION(npts,nvm)                            :: tot_litter_carb     !! Total litter carbon
389                                                                                       !! @tex $(gC m^{-2})$ @endtex
390    REAL(r_std), DIMENSION(npts,nvm)                            :: tot_soil_carb       !! Total soil carbon 
391                                                                                       !! @tex $(gC m^{-2})$ @endtex
392    REAL(r_std), DIMENSION(npts)                                :: carb_mass_variation !! Carbon Mass variation 
393                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
394    REAL(r_std), DIMENSION(npts,nvm)                            :: cn_ind              !! Crown area of individuals
395                                                                                       !! @tex $(m^{2})$ @endtex
396    REAL(r_std), DIMENSION(npts,nvm)                            :: woodmass_ind        !! Woodmass of individuals (gC)
397    REAL(r_std), DIMENSION(npts,nvm,nparts)                     :: f_alloc             !! Fraction that goes into plant part
398                                                                                       !! (0 to 1, unitless)
399    REAL(r_std), DIMENSION(npts)                                :: avail_tree          !! Space availability for trees
400                                                                                       !! (0 to 1, unitless)
401    REAL(r_std), DIMENSION(npts)                                :: avail_grass         !! Space availability for grasses
402                                                                                       !! (0 to 1, unitless)
403    INTEGER                                                     :: j
404    REAL(r_std),DIMENSION(npts)                                 :: prod10_total        !! Total products remaining in the pool
405                                                                                       !! after the annual release
406                                                                                       !! @tex $(gC m^{-2})$ @endtex
407    REAL(r_std),DIMENSION(npts)                                 :: prod100_total       !! Total products remaining in the pool
408                                                                                       !! after the annual release
409                                                                                       !! @tex $(gC m^{-2})$ @endtex
410    REAL(r_std),DIMENSION(npts)                                 :: cflux_prod_total    !! Total flux from conflux and the 10/100
411                                                                                       !! year-turnover pool
412                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
413    REAL(r_std),DIMENSION(npts,nvm)                             :: veget_max_old       !! "Maximal" coverage fraction of a PFT 
414                                                                                       !! (LAI-> infinity) on ground (unitless)
415    REAL(r_std), DIMENSION(npts,nvm)                            :: mortality           !! Fraction of individual dying this time
416                                                                                       !! step (0 to 1, unitless)
417    REAL(r_std), DIMENSION(npts)                                :: vartmp              !! Temporary variable used to add history
418    REAL(r_std), DIMENSION(npts,nvm)                            :: histvar             !! History variables
419!_ ================================================================================================================================
420
421    IF (bavard.GE.3) WRITE(numout,*) 'Entering stomate_lpj'
422
423 
424  !! 1. Initializations
425   
426    !! 1.1 Initialize variables to zero
427    co2_to_bm(:,:) = zero
428    co2_fire(:,:) = zero
429    npp_daily(:,:) = zero
430    resp_maint(:,:) = zero
431    resp_growth(:,:) = zero
432    harvest_above(:) = zero
433    bm_to_litter(:,:,:,:) = zero
434    cn_ind(:,:) = zero
435    woodmass_ind(:,:) = zero
436    turnover_daily(:,:,:,:) = zero
437   
438    !! 1.2  Initialize variables to veget_max
439    veget_max_old(:,:) = veget_max(:,:)
440
441    !! 1.3 Calculate some vegetation characteristics
442   
443    !! 1.3.1 Calculate some vegetation characteristics
444    !        Calculate cn_ind (individual crown mass) and individual height from
445    !        state variables if running DGVM or dynamic mortality in static cover mode
446    !??        Explain (maybe in the header once) why you mulitply with veget_max in the DGVM
447    !??        and why you don't multiply with veget_max in stomate.
448    IF ( control%ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
449       IF(control%ok_dgvm) THEN
450          WHERE (ind(:,:).GT.min_stomate)
451             woodmass_ind(:,:) = &
452                  ((biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
453                  +biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)) & 
454                  *veget_max(:,:))/ind(:,:)
455          ENDWHERE
456       ELSE
457          WHERE (ind(:,:).GT.min_stomate)
458             woodmass_ind(:,:) = &
459                  (biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
460                  +biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon))/ind(:,:)
461          ENDWHERE
462       ENDIF
463
464       CALL crown (npts,  PFTpresent, &
465            ind, biomass, woodmass_ind, &
466            veget_max, cn_ind, height)
467    ENDIF
468
469    !! 1.3.2 Prescribe characteristics if the vegetation is not dynamic
470    !        IF the DGVM is not activated, the density of individuals and their crown
471    !        areas don't matter, but they should be defined for the case we switch on
472    !        the DGVM afterwards. At the first call, if the DGVM is not activated,
473    !        impose a minimum biomass for prescribed PFTs and declare them present.
474    CALL prescribe (npts, &
475         veget_max, dt_days, PFTpresent, everywhere, when_growthinit, &
476         biomass, leaf_frac, ind, cn_ind, co2_to_bm)
477
478
479  !! 2. Climatic constraints for PFT presence and regenerativeness
480
481    !   Call this even when DGVM is not activated so that "adapted" and "regenerate"
482    !   are kept up to date for the moment when the DGVM is activated.
483    CALL constraints (npts, dt_days, &
484         t2m_month, t2m_min_daily,when_growthinit, &
485         adapted, regenerate)
486
487   
488  !! 3. Determine introduction and elimination of PTS based on climate criteria
489 
490    IF ( control%ok_dgvm ) THEN
491     
492       !! 3.1 Calculate introduction and elimination
493       CALL pftinout (npts, dt_days, adapted, regenerate, &
494            neighbours, veget_max, &
495            biomass, ind, cn_ind, age, leaf_frac, npp_longterm, lm_lastyearmax, senescence, &
496            PFTpresent, everywhere, when_growthinit, need_adjacent, RIP_time, &
497            co2_to_bm, &
498            avail_tree, avail_grass)
499
500       !! 3.2 Reset attributes for eliminated PFTs.
501       !     This also kills PFTs that had 0 leafmass during the last year. The message
502       !     "... after pftinout" is misleading in this case.
503       CALL kill (npts, 'pftinout  ', lm_lastyearmax, &
504            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
505            lai, age, leaf_age, leaf_frac, npp_longterm, &
506            when_growthinit, everywhere, veget_max, bm_to_litter)
507
508       
509       !! 3.3 Calculate new crown area and diameter
510       !      Calculate new crown area, diameter and maximum vegetation cover**[No longer used in the subroutine]
511       !      unsure whether this is really required
512       !      - in theory this could ONLY be done at the END of stomate_lpj
513       !      calculate woodmass of individual tree
514       WHERE ((ind(:,:).GT.min_stomate))
515          WHERE  ( veget_max(:,:) .GT. min_stomate)
516             woodmass_ind(:,:) = &
517                  ((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
518                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))*veget_max(:,:))/ind(:,:)
519          ELSEWHERE
520             woodmass_ind(:,:) =(biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
521                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
522          ENDWHERE
523
524       ENDWHERE
525       
526       ! Calculate crown area and diameter for all PFTs (including the newly established)
527       CALL crown (npts, PFTpresent, &
528            ind, biomass, woodmass_ind, &
529            veget_max, cn_ind, height)
530
531    ENDIF
532   
533  !! 4. Phenology
534
535    !! 4.1 Write values to history file
536    !      Current values for ::when_growthinit
537    CALL histwrite_p (hist_id_stomate, 'WHEN_GROWTHINIT', itime, when_growthinit, npts*nvm, horipft_index)
538
539    ! Set and write values for ::PFTpresent
540    WHERE(PFTpresent)
541       histvar=un
542    ELSEWHERE
543       histvar=zero
544    ENDWHERE
545    CALL histwrite_p (hist_id_stomate, 'PFTPRESENT', itime, histvar, npts*nvm, horipft_index)
546
547    ! Set and write values for gdd_midwinter
548    WHERE(gdd_midwinter.EQ.undef)
549       histvar=val_exp
550    ELSEWHERE
551       histvar=gdd_midwinter
552    ENDWHERE
553    CALL histwrite_p (hist_id_stomate, 'GDD_MIDWINTER', itime, histvar, npts*nvm, horipft_index)
554
555    ! Set and write values for gdd_m5_dormance
556    WHERE(gdd_m5_dormance.EQ.undef)
557       histvar=val_exp
558    ELSEWHERE
559       histvar=gdd_m5_dormance
560    ENDWHERE
561    CALL histwrite (hist_id_stomate, 'GDD_M5_DORMANCE', itime, histvar, npts*nvm, horipft_index)
562
563    ! Set and write values for ncd_dormance
564    WHERE(ncd_dormance.EQ.undef)
565       histvar=val_exp
566    ELSEWHERE
567       histvar=ncd_dormance
568    ENDWHERE
569    CALL histwrite_p (hist_id_stomate, 'NCD_DORMANCE', itime, histvar, npts*nvm, horipft_index)
570
571    !! 4.2 Calculate phenology
572    CALL phenology (npts, dt_days, PFTpresent, &
573         veget_max, &
574         tlong_ref, t2m_month, t2m_week, gpp_daily, &
575         maxmoiavail_lastyear, minmoiavail_lastyear, &
576         moiavail_month, moiavail_week, &
577         gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
578         senescence, time_hum_min, &
579         biomass, leaf_frac, leaf_age, &
580         when_growthinit, co2_to_bm)
581   
582  !! 5. Allocate C to different plant parts
583   
584    CALL alloc (npts, dt_days, &
585         lai, veget_max, senescence, when_growthinit, &
586         moiavail_week, tsoil_month, soilhum_month, &
587         biomass, age, leaf_age, leaf_frac, rprof, f_alloc)
588
589  !! 6. NPP, maintenance and growth respiration
590
591    !! 6.1 Calculate NPP and respiration terms
592    CALL npp_calc (npts, dt_days, &
593         PFTpresent, &
594         tlong_ref, t2m_daily, tsoil_daily, lai, rprof, &
595         gpp_daily, f_alloc, bm_alloc, resp_maint_part,&
596         biomass, leaf_age, leaf_frac, age, &
597         resp_maint, resp_growth, npp_daily)
598
599    !! 6.2 Kill slow growing PFTs in DGVM or STOMATE with constant mortality
600    IF ( control%ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
601       CALL kill (npts, 'npp       ', lm_lastyearmax,  &
602            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
603            lai, age, leaf_age, leaf_frac, npp_longterm, &
604            when_growthinit, everywhere, veget_max, bm_to_litter)
605
606       !! 6.2.1 Update wood biomass     
607       !        For the DGVM
608       IF(control%ok_dgvm) THEN
609          WHERE (ind(:,:).GT.min_stomate)
610             woodmass_ind(:,:) = &
611                  ((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
612                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon)) & 
613                  *veget_max(:,:))/ind(:,:)
614          ENDWHERE
615
616       ! For all pixels with individuals
617       ELSE
618          WHERE (ind(:,:).GT.min_stomate)
619             woodmass_ind(:,:) = &
620                  (biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
621                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
622          ENDWHERE
623       ENDIF ! control%ok_dgvm
624
625       !! 6.2.2 New crown area and maximum vegetation cover after growth
626       CALL crown (npts, PFTpresent, &
627            ind, biomass, woodmass_ind,&
628            veget_max, cn_ind, height)
629
630    ENDIF ! control%ok_dgvm
631   
632  !! 7. fire
633
634    !! 7.1. Burn PFTs
635    CALL fire (npts, dt_days, litterpart, &
636         litterhum_daily, t2m_daily, lignin_struc, veget_max, &
637         fireindex, firelitter, biomass, ind, &
638         litter, dead_leaves, bm_to_litter, black_carbon, &
639         co2_fire, MatrixA)
640
641    !! 7.2 Kill PFTs in DGVM
642    IF ( control%ok_dgvm ) THEN
643
644       ! reset attributes for eliminated PFTs
645       CALL kill (npts, 'fire      ', lm_lastyearmax, &
646            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
647            lai, age, leaf_age, leaf_frac, npp_longterm, &
648            when_growthinit, everywhere, veget_max, bm_to_litter)
649
650    ENDIF ! control%ok_dgvm
651 
652  !! 8. Tree mortality
653
654    ! Does not depend on age, therefore does not change crown area.
655    CALL gap (npts, dt_days, &
656         npp_longterm, turnover_longterm, lm_lastyearmax, &
657         PFTpresent, biomass, ind, bm_to_litter, mortality)
658
659    IF ( control%ok_dgvm ) THEN
660
661       ! reset attributes for eliminated PFTs
662       CALL kill (npts, 'gap       ', lm_lastyearmax, &
663            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
664            lai, age, leaf_age, leaf_frac, npp_longterm, &
665            when_growthinit, everywhere, veget_max, bm_to_litter)
666
667    ENDIF
668
669  !! 9. Calculate vcmax, vjmax and photosynthesis temperatures
670
671    CALL vmax (npts, dt_days, &
672         leaf_age, leaf_frac, &
673         vcmax, vjmax)
674
675    CALL assim_temp (npts, tlong_ref, t2m_month, &
676         t_photo_min, t_photo_opt, t_photo_max)
677
678  !! 10. Leaf senescence, new lai and other turnover processes
679
680    CALL turn (npts, dt_days, PFTpresent, &
681         herbivores, &
682         maxmoiavail_lastyear, minmoiavail_lastyear, &
683         moiavail_week,  moiavail_month,tlong_ref, t2m_month, t2m_week, veget_max, &
684         gdd_from_growthinit, leaf_age, leaf_frac, age, lai, biomass, &
685         turnover_daily, senescence,turnover_time)
686
687    !! 11. Light competition
688   
689    !! If not using constant mortality then kill with light competition
690!    IF ( control%ok_dgvm .OR. .NOT.(lpj_gap_const_mort) ) THEN
691    IF ( control%ok_dgvm ) THEN
692 
693       !! 11.1 Light competition
694       CALL light (npts, dt_days, &
695            veget_max, fpc_max, PFTpresent, cn_ind, lai, maxfpc_lastyear, &
696            lm_lastyearmax, ind, biomass, veget_lastlight, bm_to_litter, mortality)
697       
698       !! 11.2 Reset attributes for eliminated PFTs
699       CALL kill (npts, 'light     ', lm_lastyearmax, &
700            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
701            lai, age, leaf_age, leaf_frac, npp_longterm, &
702            when_growthinit, everywhere, veget_max, bm_to_litter)
703
704    ENDIF
705
706   
707  !! 12. Establishment of saplings
708   
709    IF ( control%ok_dgvm .OR. .NOT.lpj_gap_const_mort ) THEN
710
711       !! 12.1 Establish new plants
712       CALL establish (npts, dt_days, PFTpresent, regenerate, &
713            neighbours, resolution, need_adjacent, herbivores, &
714            precip_lastyear, gdd0_lastyear, lm_lastyearmax, &
715            cn_ind, lai, avail_tree, avail_grass, npp_longterm, &
716            leaf_age, leaf_frac, &
717            ind, biomass, age, everywhere, co2_to_bm, veget_max, woodmass_ind)
718
719       !! 12.2 Calculate new crown area (and maximum vegetation cover)
720       CALL crown (npts, PFTpresent, &
721            ind, biomass, woodmass_ind, &
722            veget_max, cn_ind, height)
723
724    ENDIF
725
726  !! 13. Calculate final LAI and vegetation cover
727   
728    CALL cover (npts, cn_ind, ind, biomass, &
729         veget_max, veget_max_old, lai, &
730         litter, carbon, turnover_daily, bm_to_litter)
731
732  !! 14. Update litter pools to account for harvest
733 
734    ! the whole litter stuff:
735    !    litter update, lignin content, PFT parts, litter decay,
736    !    litter heterotrophic respiration, dead leaf soil cover.
737    !    No vertical discretisation in the soil for litter decay.\n
738    ! added by shilong for harvest
739    IF(harvest_agri) THEN
740       CALL harvest(npts, dt_days, veget_max, &
741            bm_to_litter, turnover_daily, &
742            harvest_above)
743    ENDIF
744
745  !! 15. Land cover change
746
747    !shilong adde turnover_daily
748    IF(EndOfYear) THEN
749       IF (lcchange) THEN
750          CALL lcchange_main (npts, dt_days, veget_max, veget_max_new, &
751               biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
752               co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,flux10,flux100, &
753!!$               prod10,prod100,prod10_total,prod100_total,&
754!!$               convflux,cflux_prod_total,cflux_prod10,cflux_prod100,leaf_frac,&
755               prod10,prod100,convflux,cflux_prod10,cflux_prod100,leaf_frac,&
756               npp_longterm, lm_lastyearmax, litter, carbon)
757       ENDIF
758    ENDIF
759    !MM déplacement pour initialisation correcte des grandeurs cumulées :
760    cflux_prod_total(:) = convflux(:) + cflux_prod10(:) + cflux_prod100(:)
761    prod10_total(:)=SUM(prod10,dim=2)
762    prod100_total(:)=SUM(prod100,dim=2)
763   
764  !! 16. Total heterotrophic respiration
765
766    tot_soil_carb(:,:) = zero
767    tot_litter_carb(:,:) = zero
768    DO j=2,nvm
769
770       tot_litter_carb(:,j) = tot_litter_carb(:,j) + (litter(:,istructural,j,iabove,icarbon) + &
771            &          litter(:,imetabolic,j,iabove,icarbon) + &
772            &          litter(:,istructural,j,ibelow,icarbon) + litter(:,imetabolic,j,ibelow,icarbon))
773
774       tot_soil_carb(:,j) = tot_soil_carb(:,j) + (carbon(:,iactive,j) + &
775            &          carbon(:,islow,j)+  carbon(:,ipassive,j))
776
777    ENDDO
778    tot_litter_soil_carb(:,:) = tot_litter_carb(:,:) + tot_soil_carb(:,:)
779
780!!$     DO k = 1, nelements ! Loop over # elements
781!!$        tot_live_biomass(:,:,k) = biomass(:,:,ileaf,k) + biomass(:,:,isapabove,k) + biomass(:,:,isapbelow,k) +&
782!!$             &                    biomass(:,:,iheartabove,k) + biomass(:,:,iheartbelow,k) + &
783!!$             &                    biomass(:,:,iroot,k) + biomass(:,:,ifruit,k) + biomass(:,:,icarbres,k)
784!!$    END DO ! Loop over # elements
785
786    tot_live_biomass(:,:,:) = biomass(:,:,ileaf,:) + biomass(:,:,isapabove,:) + biomass(:,:,isapbelow,:) +&
787             &                    biomass(:,:,iheartabove,:) + biomass(:,:,iheartbelow,:) + &
788             &                    biomass(:,:,iroot,:) + biomass(:,:,ifruit,:) + biomass(:,:,icarbres,:)
789
790
791    tot_turnover(:,:,:) = turnover_daily(:,:,ileaf,:) + turnover_daily(:,:,isapabove,:) + &
792         &         turnover_daily(:,:,isapbelow,:) + turnover_daily(:,:,iheartabove,:) + &
793         &         turnover_daily(:,:,iheartbelow,:) + turnover_daily(:,:,iroot,:) + &
794         &         turnover_daily(:,:,ifruit,:) + turnover_daily(:,:,icarbres,:)
795
796    tot_bm_to_litter(:,:,:) = bm_to_litter(:,:,ileaf,:) + bm_to_litter(:,:,isapabove,:) +&
797         &             bm_to_litter(:,:,isapbelow,:) + bm_to_litter(:,:,iheartbelow,:) +&
798         &             bm_to_litter(:,:,iheartabove,:) + bm_to_litter(:,:,iroot,:) + &
799         &             bm_to_litter(:,:,ifruit,:) + bm_to_litter(:,:,icarbres,:)
800
801    carb_mass_variation(:)=-carb_mass_total(:)
802    carb_mass_total(:)=SUM((tot_live_biomass(:,:,icarbon)+tot_litter_carb+tot_soil_carb)*veget_max,dim=2) + &
803         &                 (prod10_total + prod100_total)
804    carb_mass_variation(:)=carb_mass_total(:)+carb_mass_variation(:)
805   
806  !! 17. Write history
807
808    CALL histwrite_p (hist_id_stomate, 'RESOLUTION_X', itime, &
809         resolution(:,1), npts, hori_index)
810    CALL histwrite_p (hist_id_stomate, 'RESOLUTION_Y', itime, &
811         resolution(:,2), npts, hori_index)
812    CALL histwrite_p (hist_id_stomate, 'CONTFRAC', itime, &
813         contfrac(:), npts, hori_index)
814
815    CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AB', itime, &
816         litter(:,istructural,:,iabove,icarbon), npts*nvm, horipft_index)
817    CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AB', itime, &
818         litter(:,imetabolic,:,iabove,icarbon), npts*nvm, horipft_index)
819    CALL histwrite_p (hist_id_stomate, 'LITTER_STR_BE', itime, &
820         litter(:,istructural,:,ibelow,icarbon), npts*nvm, horipft_index)
821    CALL histwrite_p (hist_id_stomate, 'LITTER_MET_BE', itime, &
822         litter(:,imetabolic,:,ibelow,icarbon), npts*nvm, horipft_index)
823
824    CALL histwrite_p (hist_id_stomate, 'DEADLEAF_COVER', itime, &
825         deadleaf_cover, npts, hori_index)
826
827    CALL histwrite_p (hist_id_stomate, 'TOTAL_SOIL_CARB', itime, &
828         tot_litter_soil_carb, npts*nvm, horipft_index)
829    CALL histwrite_p (hist_id_stomate, 'CARBON_ACTIVE', itime, &
830         carbon(:,iactive,:), npts*nvm, horipft_index)
831    CALL histwrite_p (hist_id_stomate, 'CARBON_SLOW', itime, &
832         carbon(:,islow,:), npts*nvm, horipft_index)
833    CALL histwrite_p (hist_id_stomate, 'CARBON_PASSIVE', itime, &
834         carbon(:,ipassive,:), npts*nvm, horipft_index)
835
836    CALL histwrite_p (hist_id_stomate, 'T2M_MONTH', itime, &
837         t2m_month, npts, hori_index)
838    CALL histwrite_p (hist_id_stomate, 'T2M_WEEK', itime, &
839         t2m_week, npts, hori_index)
840
841    CALL histwrite_p (hist_id_stomate, 'HET_RESP', itime, &
842         resp_hetero(:,:), npts*nvm, horipft_index)
843
844    CALL histwrite_p (hist_id_stomate, 'BLACK_CARBON', itime, &
845         black_carbon, npts, hori_index)
846
847    CALL histwrite_p (hist_id_stomate, 'FIREINDEX', itime, &
848         fireindex(:,:), npts*nvm, horipft_index)
849    CALL histwrite_p (hist_id_stomate, 'LITTERHUM', itime, &
850         litterhum_daily, npts, hori_index)
851    CALL histwrite_p (hist_id_stomate, 'CO2_FIRE', itime, &
852         co2_fire, npts*nvm, horipft_index)
853    CALL histwrite_p (hist_id_stomate, 'CO2_TAKEN', itime, &
854         co2_to_bm, npts*nvm, horipft_index)
855    ! land cover change
856    CALL histwrite_p (hist_id_stomate, 'CONVFLUX', itime, &
857         convflux, npts, hori_index)
858    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10', itime, &
859         cflux_prod10, npts, hori_index)
860    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100', itime, &
861         cflux_prod100, npts, hori_index)
862    CALL histwrite_p (hist_id_stomate, 'HARVEST_ABOVE', itime, &
863         harvest_above, npts, hori_index)
864
865    CALL histwrite_p (hist_id_stomate, 'LAI', itime, &
866         lai, npts*nvm, horipft_index)
867    CALL histwrite_p (hist_id_stomate, 'VEGET_MAX', itime, &
868         veget_max, npts*nvm, horipft_index)
869    CALL histwrite_p (hist_id_stomate, 'NPP', itime, &
870         npp_daily, npts*nvm, horipft_index)
871    CALL histwrite_p (hist_id_stomate, 'GPP', itime, &
872         gpp_daily, npts*nvm, horipft_index)
873    CALL histwrite_p (hist_id_stomate, 'IND', itime, &
874         ind, npts*nvm, horipft_index)
875    CALL histwrite_p (hist_id_stomate, 'CN_IND', itime, &
876         cn_ind, npts*nvm, horipft_index)
877    CALL histwrite_p (hist_id_stomate, 'WOODMASS_IND', itime, &
878         woodmass_ind, npts*nvm, horipft_index)
879    CALL histwrite_p (hist_id_stomate, 'TOTAL_M', itime, &
880         tot_live_biomass(:,:,icarbon), npts*nvm, horipft_index)
881    CALL histwrite_p (hist_id_stomate, 'LEAF_M', itime, &
882         biomass(:,:,ileaf,icarbon), npts*nvm, horipft_index)
883    CALL histwrite_p (hist_id_stomate, 'SAP_M_AB', itime, &
884         biomass(:,:,isapabove,icarbon), npts*nvm, horipft_index)
885    CALL histwrite_p (hist_id_stomate, 'SAP_M_BE', itime, &
886         biomass(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
887    CALL histwrite_p (hist_id_stomate, 'HEART_M_AB', itime, &
888         biomass(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
889    CALL histwrite_p (hist_id_stomate, 'HEART_M_BE', itime, &
890         biomass(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
891    CALL histwrite_p (hist_id_stomate, 'ROOT_M', itime, &
892         biomass(:,:,iroot,icarbon), npts*nvm, horipft_index)
893    CALL histwrite_p (hist_id_stomate, 'FRUIT_M', itime, &
894         biomass(:,:,ifruit,icarbon), npts*nvm, horipft_index)
895    CALL histwrite_p (hist_id_stomate, 'RESERVE_M', itime, &
896         biomass(:,:,icarbres,icarbon), npts*nvm, horipft_index)
897    CALL histwrite_p (hist_id_stomate, 'TOTAL_TURN', itime, &
898         tot_turnover(:,:,icarbon), npts*nvm, horipft_index)
899    CALL histwrite_p (hist_id_stomate, 'LEAF_TURN', itime, &
900         turnover_daily(:,:,ileaf,icarbon), npts*nvm, horipft_index)
901    CALL histwrite_p (hist_id_stomate, 'SAP_AB_TURN', itime, &
902         turnover_daily(:,:,isapabove,icarbon), npts*nvm, horipft_index)
903    CALL histwrite_p (hist_id_stomate, 'ROOT_TURN', itime, &
904         turnover_daily(:,:,iroot,icarbon), npts*nvm, horipft_index)
905    CALL histwrite_p (hist_id_stomate, 'FRUIT_TURN', itime, &
906         turnover_daily(:,:,ifruit,icarbon), npts*nvm, horipft_index)
907    CALL histwrite_p (hist_id_stomate, 'TOTAL_BM_LITTER', itime, &
908         tot_bm_to_litter(:,:,icarbon), npts*nvm, horipft_index)
909    CALL histwrite_p (hist_id_stomate, 'LEAF_BM_LITTER', itime, &
910         bm_to_litter(:,:,ileaf,icarbon), npts*nvm, horipft_index)
911    CALL histwrite_p (hist_id_stomate, 'SAP_AB_BM_LITTER', itime, &
912         bm_to_litter(:,:,isapabove,icarbon), npts*nvm, horipft_index)
913    CALL histwrite_p (hist_id_stomate, 'SAP_BE_BM_LITTER', itime, &
914         bm_to_litter(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
915    CALL histwrite_p (hist_id_stomate, 'HEART_AB_BM_LITTER', itime, &
916         bm_to_litter(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
917    CALL histwrite_p (hist_id_stomate, 'HEART_BE_BM_LITTER', itime, &
918         bm_to_litter(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
919    CALL histwrite_p (hist_id_stomate, 'ROOT_BM_LITTER', itime, &
920         bm_to_litter(:,:,iroot,icarbon), npts*nvm, horipft_index)
921    CALL histwrite_p (hist_id_stomate, 'FRUIT_BM_LITTER', itime, &
922         bm_to_litter(:,:,ifruit,icarbon), npts*nvm, horipft_index)
923    CALL histwrite_p (hist_id_stomate, 'RESERVE_BM_LITTER', itime, &
924         bm_to_litter(:,:,icarbres,icarbon), npts*nvm, horipft_index)
925    CALL histwrite_p (hist_id_stomate, 'MAINT_RESP', itime, &
926         resp_maint, npts*nvm, horipft_index)
927    CALL histwrite_p (hist_id_stomate, 'GROWTH_RESP', itime, &
928         resp_growth, npts*nvm, horipft_index)
929    CALL histwrite_p (hist_id_stomate, 'AGE', itime, &
930         age, npts*nvm, horipft_index)
931    CALL histwrite_p (hist_id_stomate, 'HEIGHT', itime, &
932         height, npts*nvm, horipft_index)
933    CALL histwrite_p (hist_id_stomate, 'MOISTRESS', itime, &
934         moiavail_week, npts*nvm, horipft_index)
935    CALL histwrite_p (hist_id_stomate, 'VCMAX', itime, &
936         vcmax, npts*nvm, horipft_index)
937    CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME', itime, &
938         turnover_time, npts*nvm, horipft_index)
939    ! land cover change
940    CALL histwrite_p (hist_id_stomate, 'PROD10', itime, &
941         prod10, npts*11, horip11_index)
942    CALL histwrite_p (hist_id_stomate, 'PROD100', itime, &
943         prod100, npts*101, horip101_index)
944    CALL histwrite_p (hist_id_stomate, 'FLUX10', itime, &
945         flux10, npts*10, horip10_index)
946    CALL histwrite_p (hist_id_stomate, 'FLUX100', itime, &
947         flux100, npts*100, horip100_index)
948
949    IF ( hist_id_stomate_IPCC > 0 ) THEN
950       vartmp(:)=SUM(tot_live_biomass(:,:,icarbon)*veget_max,dim=2)/1e3*contfrac
951       CALL histwrite_p (hist_id_stomate_IPCC, "cVeg", itime, &
952            vartmp, npts, hori_index)
953       vartmp(:)=SUM(tot_litter_carb*veget_max,dim=2)/1e3*contfrac
954       CALL histwrite_p (hist_id_stomate_IPCC, "cLitter", itime, &
955            vartmp, npts, hori_index)
956       vartmp(:)=SUM(tot_soil_carb*veget_max,dim=2)/1e3*contfrac
957       CALL histwrite_p (hist_id_stomate_IPCC, "cSoil", itime, &
958            vartmp, npts, hori_index)
959       vartmp(:)=(prod10_total + prod100_total)/1e3
960       CALL histwrite_p (hist_id_stomate_IPCC, "cProduct", itime, &
961            vartmp, npts, hori_index)
962       vartmp(:)=carb_mass_variation/1e3/one_day*contfrac
963       CALL histwrite_p (hist_id_stomate_IPCC, "cMassVariation", itime, &
964            vartmp, npts, hori_index)
965       vartmp(:)=SUM(lai*veget_max,dim=2)*contfrac
966       CALL histwrite_p (hist_id_stomate_IPCC, "lai", itime, &
967            vartmp, npts, hori_index)
968       vartmp(:)=SUM(gpp_daily*veget_max,dim=2)/1e3/one_day*contfrac
969       CALL histwrite_p (hist_id_stomate_IPCC, "gpp", itime, &
970            vartmp, npts, hori_index)
971       vartmp(:)=SUM((resp_maint+resp_growth)*veget_max,dim=2)/1e3/one_day*contfrac
972       CALL histwrite_p (hist_id_stomate_IPCC, "ra", itime, &
973            vartmp, npts, hori_index)
974       vartmp(:)=SUM(npp_daily*veget_max,dim=2)/1e3/one_day*contfrac
975       CALL histwrite_p (hist_id_stomate_IPCC, "npp", itime, &
976            vartmp, npts, hori_index)
977       vartmp(:)=SUM(resp_hetero*veget_max,dim=2)/1e3/one_day*contfrac
978       CALL histwrite_p (hist_id_stomate_IPCC, "rh", itime, &
979            vartmp, npts, hori_index)
980       vartmp(:)=SUM(co2_fire*veget_max,dim=2)/1e3/one_day*contfrac
981       CALL histwrite_p (hist_id_stomate_IPCC, "fFire", itime, &
982            vartmp, npts, hori_index)
983       vartmp(:)=harvest_above/1e3/one_day*contfrac
984       CALL histwrite_p (hist_id_stomate_IPCC, "fHarvest", itime, &
985            vartmp, npts, hori_index)
986       vartmp(:)=cflux_prod_total/1e3/one_day*contfrac
987       CALL histwrite_p (hist_id_stomate_IPCC, "fLuc", itime, &
988            vartmp, npts, hori_index)
989       vartmp(:)=(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
990            &        *veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_day*contfrac
991       CALL histwrite_p (hist_id_stomate_IPCC, "nbp", itime, &
992            vartmp, npts, hori_index)
993       vartmp(:)=SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*veget_max,dim=2)/1e3/one_day*contfrac
994       CALL histwrite_p (hist_id_stomate_IPCC, "fVegLitter", itime, &
995            vartmp, npts, hori_index)
996       vartmp(:)=SUM(SUM(soilcarbon_input,dim=2)*veget_max,dim=2)/1e3/one_day*contfrac
997       CALL histwrite_p (hist_id_stomate_IPCC, "fLitterSoil", itime, &
998            vartmp, npts, hori_index)
999       vartmp(:)=SUM(biomass(:,:,ileaf,icarbon)*veget_max,dim=2)/1e3*contfrac
1000       CALL histwrite_p (hist_id_stomate_IPCC, "cLeaf", itime, &
1001            vartmp, npts, hori_index)
1002       vartmp(:)=SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*veget_max,dim=2)/1e3*contfrac
1003       CALL histwrite_p (hist_id_stomate_IPCC, "cWood", itime, &
1004            vartmp, npts, hori_index)
1005       vartmp(:)=SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + biomass(:,:,iheartbelow,icarbon) ) &
1006            &        *veget_max,dim=2)/1e3*contfrac
1007       CALL histwrite_p (hist_id_stomate_IPCC, "cRoot", itime, &
1008            vartmp, npts, hori_index)
1009       vartmp(:)=SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*veget_max,dim=2)/1e3*contfrac
1010       CALL histwrite_p (hist_id_stomate_IPCC, "cMisc", itime, &
1011            vartmp, npts, hori_index)
1012       vartmp(:)=SUM((litter(:,istructural,:,iabove,icarbon)+litter(:,imetabolic,:,iabove,icarbon))*veget_max,dim=2)/1e3*contfrac
1013       CALL histwrite_p (hist_id_stomate_IPCC, "cLitterAbove", itime, &
1014            vartmp, npts, hori_index)
1015       vartmp(:)=SUM((litter(:,istructural,:,ibelow,icarbon)+litter(:,imetabolic,:,ibelow,icarbon))*veget_max,dim=2)/1e3*contfrac
1016       CALL histwrite_p (hist_id_stomate_IPCC, "cLitterBelow", itime, &
1017            vartmp, npts, hori_index)
1018       vartmp(:)=SUM(carbon(:,iactive,:)*veget_max,dim=2)/1e3*contfrac
1019       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilFast", itime, &
1020            vartmp, npts, hori_index)
1021       vartmp(:)=SUM(carbon(:,islow,:)*veget_max,dim=2)/1e3*contfrac
1022       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilMedium", itime, &
1023            vartmp, npts, hori_index)
1024       vartmp(:)=SUM(carbon(:,ipassive,:)*veget_max,dim=2)/1e3*contfrac
1025       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilSlow", itime, &
1026            vartmp, npts, hori_index)
1027       DO j=1,nvm
1028          histvar(:,j)=veget_max(:,j)*contfrac(:)*100
1029       ENDDO
1030       CALL histwrite_p (hist_id_stomate_IPCC, "landCoverFrac", itime, &
1031            histvar, npts*nvm, horipft_index)
1032       !-
1033       vartmp(:)=zero
1034       DO j = 2,nvm
1035          IF (is_deciduous(j)) THEN
1036             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
1037          ENDIF
1038       ENDDO
1039       CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimDec", itime, &
1040            vartmp, npts, hori_index)
1041       !-
1042       vartmp(:)=zero
1043       DO j = 2,nvm
1044          IF (is_evergreen(j)) THEN
1045             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
1046          ENDIF
1047       ENDDO
1048       CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimEver", itime, &
1049            vartmp, npts, hori_index)
1050       !-
1051       vartmp(:)=zero
1052       DO j = 2,nvm
1053          IF ( .NOT.(is_c4(j)) ) THEN
1054             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
1055          ENDIF
1056       ENDDO
1057       CALL histwrite_p (hist_id_stomate_IPCC, "c3PftFrac", itime, &
1058            vartmp, npts, hori_index)
1059       !-
1060       vartmp(:)=zero
1061       DO j = 2,nvm
1062          IF ( is_c4(j) ) THEN
1063             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
1064          ENDIF
1065       ENDDO
1066       CALL histwrite_p (hist_id_stomate_IPCC, "c4PftFrac", itime, &
1067            vartmp, npts, hori_index)
1068       !-
1069       vartmp(:)=SUM(resp_growth*veget_max,dim=2)/1e3/one_day*contfrac
1070       CALL histwrite_p (hist_id_stomate_IPCC, "rGrowth", itime, &
1071            vartmp, npts, hori_index)
1072       vartmp(:)=SUM(resp_maint*veget_max,dim=2)/1e3/one_day*contfrac
1073       CALL histwrite_p (hist_id_stomate_IPCC, "rMaint", itime, &
1074            vartmp, npts, hori_index)
1075       vartmp(:)=SUM(bm_alloc(:,:,ileaf,icarbon)*veget_max,dim=2)/1e3/one_day*contfrac
1076       CALL histwrite_p (hist_id_stomate_IPCC, "nppLeaf", itime, &
1077            vartmp, npts, hori_index)
1078       vartmp(:)=SUM(bm_alloc(:,:,isapabove,icarbon)*veget_max,dim=2)/1e3/one_day*contfrac
1079       CALL histwrite_p (hist_id_stomate_IPCC, "nppWood", itime, &
1080            vartmp, npts, hori_index)
1081       vartmp(:)=SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )*veget_max,dim=2)/1e3/one_day*contfrac
1082       CALL histwrite_p (hist_id_stomate_IPCC, "nppRoot", itime, &
1083            vartmp, npts, hori_index)
1084
1085       CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_X', itime, &
1086            resolution(:,1), npts, hori_index)
1087       CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_Y', itime, &
1088            resolution(:,2), npts, hori_index)
1089       CALL histwrite_p (hist_id_stomate_IPCC, 'CONTFRAC', itime, &
1090            contfrac(:), npts, hori_index)
1091
1092    ENDIF
1093
1094    IF (bavard.GE.4) WRITE(numout,*) 'Leaving stomate_lpj'
1095
1096  END SUBROUTINE StomateLpj
1097
1098
1099!! ================================================================================================================================
1100!! SUBROUTINE   : harvest
1101!!
1102!>\BRIEF        Harvest of croplands
1103!!
1104!! DESCRIPTION  : To take into account biomass harvest from crop (mainly to take
1105!! into account for the reduced litter input and then decreased soil carbon. it is a
1106!! constant (40\%) fraction of above ground biomass.
1107!!
1108!! RECENT CHANGE(S) : None
1109!!
1110!! MAIN OUTPUT VARIABLE(S): ::harvest_above the harvested biomass
1111!!
1112!! REFERENCE(S) :
1113!! - Piao, S., P. Ciais, P. Friedlingstein, N. de Noblet-Ducoudre, P. Cadule, N. Viovy, and T. Wang. 2009.
1114!!   Spatiotemporal patterns of terrestrial carbon cycle during the 20th century. Global Biogeochemical
1115!!   Cycles 23:doi:10.1029/2008GB003339.
1116!!
1117!! FLOWCHART    : None
1118!! \n
1119!_ ================================================================================================================================
1120
1121  SUBROUTINE harvest(npts, dt_days, veget_max, &
1122       bm_to_litter, turnover_daily, &
1123       harvest_above)
1124
1125  !! 0. Variable and parameter declaration
1126
1127    !! 0.1 Input variables
1128
1129    INTEGER, INTENT(in)                                    :: npts            !! Domain size (unitless)
1130    REAL(r_std), INTENT(in)                                :: dt_days         !! Time step (days)                               
1131    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: veget_max       !! new "maximal" coverage fraction of a PFT (LAI ->
1132                                                                              !! infinity) on ground @tex $(m^2 m^{-2})$ @endtex
1133   
1134   !! 0.2 Output variables
1135   
1136   !! 0.3 Modified variables
1137
1138    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! [DISPENSABLE] conversion of biomass to litter
1139                                                                                     !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1140    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily   !! Turnover rates
1141                                                                               !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1142    REAL(r_std), DIMENSION(npts), INTENT(inout)            :: harvest_above    !! harvest above ground biomass for agriculture
1143                                                                               !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1144    !! 0.4 Local variables
1145
1146    INTEGER(i_std)                                         :: i, j, k, l, m    !! indices                       
1147    REAL(r_std)                                            :: above_old        !! biomass of previous time step
1148                                                                               !! @tex $(gC m^{-2})$ @endtex
1149!_ ================================================================================================================================
1150
1151  !! 1. Yearly initialisation
1152
1153    above_old             = zero
1154    harvest_above         = zero
1155
1156    DO i = 1, npts
1157       DO j = 1,nvm
1158          IF (.NOT. natural(j)) THEN
1159             above_old = turnover_daily(i,j,ileaf,icarbon) + turnover_daily(i,j,isapabove,icarbon) + &
1160                  &       turnover_daily(i,j,iheartabove,icarbon) + turnover_daily(i,j,ifruit,icarbon) + &
1161                  &       turnover_daily(i,j,icarbres,icarbon) + turnover_daily(i,j,isapbelow,icarbon) + &
1162                  &       turnover_daily(i,j,iheartbelow,icarbon) + turnover_daily(i,j,iroot,icarbon)
1163
1164             turnover_daily(i,j,ileaf,icarbon) = turnover_daily(i,j,ileaf,icarbon)*frac_turnover_daily
1165             turnover_daily(i,j,isapabove,icarbon) = turnover_daily(i,j,isapabove,icarbon)*frac_turnover_daily
1166             turnover_daily(i,j,isapbelow,icarbon) = turnover_daily(i,j,isapbelow,icarbon)*frac_turnover_daily
1167             turnover_daily(i,j,iheartabove,icarbon) = turnover_daily(i,j,iheartabove,icarbon)*frac_turnover_daily
1168             turnover_daily(i,j,iheartbelow,icarbon) = turnover_daily(i,j,iheartbelow,icarbon)*frac_turnover_daily
1169             turnover_daily(i,j,iroot,icarbon) = turnover_daily(i,j,iroot,icarbon)*frac_turnover_daily
1170             turnover_daily(i,j,ifruit,icarbon) = turnover_daily(i,j,ifruit,icarbon)*frac_turnover_daily
1171             turnover_daily(i,j,icarbres,icarbon) = turnover_daily(i,j,icarbres,icarbon)*frac_turnover_daily
1172             harvest_above(i)  = harvest_above(i) + veget_max(i,j) * above_old *(un - frac_turnover_daily)
1173          ENDIF
1174       ENDDO
1175    ENDDO
1176
1177!!$    harvest_above = harvest_above
1178  END SUBROUTINE harvest
1179END MODULE stomate_lpj
Note: See TracBrowser for help on using the repository browser.