source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/stomate_lpj.f90 @ 5816

Last change on this file since 5816 was 5816, checked in by jinfeng.chang, 5 years ago

copy ORCHIDEE-GMv3.2 for publication

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