source: branches/publications/ORCHIDEE_gmd-2018-261/src_stomate/stomate_lpj.f90 @ 8793

Last change on this file since 8793 was 4998, checked in by nicolas.vuichard, 7 years ago

rev29012018

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