source: branches/publications/ORCHIDEE-Clateral/src_stomate/stomate_lpj.f90 @ 7346

Last change on this file since 7346 was 7191, checked in by haicheng.zhang, 3 years ago

Implementing latral transfers of sediment and POC from land to ocean through river in ORCHILEAK
  • Property svn:keywords set to HeadURL Date Author Revision
File size: 109.3 KB
Line 
1! ================================================================================================================================
2! MODULE       : stomate_lpj
3!
4! CONTACT      : orchidee-help _at_ listes.ipsl.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  USE stomate_woodharvest
55
56
57  IMPLICIT NONE
58
59  ! private & public routines
60
61  PRIVATE
62  PUBLIC StomateLpj,StomateLpj_clear
63
64CONTAINS
65
66
67!! ================================================================================================================================
68!! SUBROUTINE   : StomateLpj_clear
69!!
70!>\BRIEF        Re-initialisation of variable
71!!
72!! DESCRIPTION  : This subroutine reinitializes variables. To be used if we want to relaunch
73!! ORCHIDEE but the routine is not used in current version.
74!!
75!! RECENT CHANGE(S) : None
76!!
77!! MAIN OUTPUT VARIABLE(S): None
78!!
79!! REFERENCE(S) : None
80!!
81!! FLOWCHART    : None
82!! \n
83!_ ================================================================================================================================
84
85  SUBROUTINE StomateLpj_clear
86
87    CALL prescribe_clear
88    CALL phenology_clear
89    CALL npp_calc_clear
90    CALL turn_clear
91    CALL soilcarbon_clear
92    CALL constraints_clear
93    CALL establish_clear
94    CALL fire_clear
95    CALL gap_clear
96    CALL light_clear
97    CALL pftinout_clear
98    CALL alloc_clear
99  END SUBROUTINE StomateLpj_clear
100
101
102!! ================================================================================================================================
103!! SUBROUTINE   : StomateLPJ
104!!
105!>\BRIEF        Main entry point for daily processes in STOMATE and LPJ, structures the call sequence
106!!              to the different processes such as dispersion, establishment, competition and mortality of PFT's.
107!!
108!! DESCRIPTION  : This routine is the main entry point to all processes calculated on a
109!! daily time step. Is mainly devoted to call the different STOMATE and LPJ routines
110!! depending of the ok_dgvm (is dynamic veg used) and lpj_constant_mortality (is background mortality used).
111!! It also prepares the cumulative
112!! fluxes or pools (e.g TOTAL_M TOTAL_BM_LITTER etc...)
113!!
114!! This routine makes frequent use of "weekly", "monthly" and "long term" variables. Quotion is used because
115!! by default "weekly" denotes 7 days, by default "monthly" denotes 20 days and by default "Long term" denotes
116!! 3 years. dtslow refers to 24 hours (1 day).
117!!
118!!
119!! RECENT CHANGE(S) : None
120!!
121!! MAIN OUTPUT VARIABLE(S): All variables related to stomate and required for LPJ dynamic vegetation mode.
122!!
123!! REFERENCE(S) :
124!! - Krinner, G., N. Viovy, N. de Noblet-Ducoudré, J. Ogeé, J. Polcher, P. Friedlingstein, P. Ciais, S. Sitch,
125!! and I. C. Prentice. 2005. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere
126!! system. Global Biogeochemical Cycles 19:GB1015, doi:1010.1029/2003GB002199.
127!! - Sitch, S., B. Smith, I. C. Prentice, A. Arneth, A. Bondeau, W. Cramer, J. O. Kaplan, S. Levis, W. Lucht,
128!! M. T. Sykes, K. Thonicke, and S. Venevsky. 2003. Evaluation of ecosystem dynamics, plant geography and
129!! terrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biology 9:161-185.
130!!
131!! FLOWCHART    : Update with existing flowchart from N Viovy (Jan 19, 2012)
132!! \n
133!_ ================================================================================================================================
134 
135  SUBROUTINE StomateLpj (npts, dt_days, &
136       neighbours, resolution, &
137       clay, herbivores, &
138       tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
139       litterhum_daily, soilhum_daily, &
140       maxmoiavail_lastyear, minmoiavail_lastyear, &
141       gdd0_lastyear, precip_lastyear, &
142       moiavail_month, moiavail_week, t2m_longterm, t2m_month, t2m_week, &
143       tsoil_month, soilhum_month, &
144       gdd_m5_dormance, gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
145       turnover_longterm, gpp_daily, &
146       time_hum_min, maxfpc_lastyear, resp_maint_part, &
147       PFTpresent, age, fireindex, firelitter, &
148       leaf_age, leaf_frac, biomass, ind, adapted, regenerate, &
149       senescence, when_growthinit, &
150       litterpart, litter_above, litter_below, depth_deepsoil, &
151       dead_leaves, carbon, DOC, DOC_EXP, lignin_struc_above,&
152       veget_cov_max, veget_cov_max_new, woodharvest, fraclut, npp_longterm, lm_lastyearmax, veget_lastlight, &
153       everywhere, need_adjacent, RIP_time, &
154       lai, rprof,npp_daily, turnover_daily, turnover_time,&
155       soilcarbon_input, &
156       co2_to_bm, co2_fire, resp_hetero, tot_soil_resp, resp_maint, resp_growth, &
157       height, deadleaf_cover, vcmax, &
158       bm_to_litter, &
159       prod10,prod100,flux10, flux100, &
160       convflux,cflux_prod10,cflux_prod100, &
161       prod10_harvest,prod100_harvest,flux10_harvest, flux100_harvest, &
162       convflux_harvest,cflux_prod10_harvest,cflux_prod100_harvest, woodharvestpft, &
163       convfluxpft, fDeforestToProduct, fLulccResidue,fHarvestToProduct, &
164       harvest_above, carb_mass_total, fpc_max, MatrixA, &
165       Tseason, Tmin_spring_time, begin_leaves, onset_date, moist_soil)
166   
167  !! 0. Variable and parameter declaration
168
169    !! 0.1 input
170
171    INTEGER(i_std), INTENT(in)                                 :: npts                 !! Domain size (unitless)
172    REAL(r_std), INTENT(in)                                    :: dt_days              !! Time step of Stomate (days)
173    INTEGER(i_std), DIMENSION(npts,NbNeighb), INTENT(in)       :: neighbours           !! Indices of the 8 neighbours of each grid
174                                                                                       !! point [1=North and then clockwise]
175    REAL(r_std), DIMENSION(npts,2), INTENT(in)                 :: resolution           !! Resolution at each grid point (m) 
176                                                                                       !! [1=E-W, 2=N-S]
177    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: clay                 !! Clay fraction (0 to 1, unitless)
178    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: herbivores           !! Time constant of probability of a leaf to
179                                                                                       !! be eaten by a herbivore (days)
180    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: tsurf_daily          !! Daily surface temperatures (K)
181    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)              :: tsoil_daily          !! Daily soil temperatures (K)
182    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_daily            !! Daily 2 meter temperatures (K)
183    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_min_daily        !! Daily minimum 2 meter temperatures (K)
184    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: litterhum_daily      !! Daily litter humidity (0 to 1, unitless)
185    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)              :: soilhum_daily        !! Daily soil humidity (0 to 1, unitless)
186    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: maxmoiavail_lastyear !! Last year's maximum moisture availability
187                                                                                       !! (0 to 1, unitless)
188    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: minmoiavail_lastyear !! Last year's minimum moisture availability
189                                                                                       !! (0 to 1, unitless)
190    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: gdd0_lastyear        !! Last year's GDD0 (K)
191    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: precip_lastyear      !! Lastyear's precipitation
192                                                                                       !! @tex $(mm year^{-1})$ @endtex
193                                                                                       !! to determine if establishment possible
194    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: moiavail_month       !! "Monthly" moisture availability (0 to 1,
195                                                                                       !! unitless)
196    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: moiavail_week        !! "Weekly" moisture availability
197                                                                                       !! (0 to 1, unitless)
198    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_longterm         !! "Long term" 2 meter reference
199                                                                                       !! temperatures (K)
200    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_month            !! "Monthly" 2-meter temperatures (K)
201    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_week             !! "Weekly" 2-meter temperatures (K)
202    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: Tseason              !! "seasonal" 2-meter temperatures (K)
203    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: Tmin_spring_time     !! Number of days after begin_leaves (leaf onset)
204    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: onset_date           !! Date in the year at when the leaves started to grow(begin_leaves)
205    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)              :: tsoil_month          !! "Monthly" soil temperatures (K)
206    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)              :: soilhum_month        !! "Monthly" soil humidity
207                                                                                       !! (0 to 1, unitless)
208    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gdd_m5_dormance      !! Growing degree days (K), threshold -5 deg
209                                                                                       !! C (for phenology)
210    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gdd_from_growthinit  !! growing degree days, since growthinit for crops
211    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: gdd_midwinter        !! Growing degree days (K), since midwinter
212                                                                                       !! (for phenology) - this is written to the history files
213    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: ncd_dormance         !! Number of chilling days (days), since
214                                                                                       !! leaves were lost (for phenology)
215    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: ngd_minus5           !! Number of growing days (days), threshold
216                                                                                       !! -5 deg C (for phenology)
217    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in):: turnover_longterm    !! "Long term" turnover rate 
218                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
219    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: gpp_daily            !! Daily gross primary productivity 
220                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
221    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: time_hum_min         !! Time elapsed since strongest moisture
222                                                                                       !! availability (days)
223    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: maxfpc_lastyear      !! Last year's maximum foliage projected
224                                                                                       !! coverage for each natural PFT,
225                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
226    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)        :: resp_maint_part      !! Maintenance respiration of different
227                                                                                       !! plant parts 
228                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
229    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: fpc_max              !! "Maximal" coverage fraction of a PFT (LAI
230                                                                                       !! -> infinity) on ground 
231                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
232    REAL(r_std), DIMENSION(npts,nvm),INTENT(in)                :: veget_cov_max_new    !! New "maximal" coverage fraction of a PFT
233    REAL(r_std), DIMENSION(npts),INTENT(in)                    :: woodharvest          !! Harvested wood biomass (gC m-2 yr-1)
234    REAL(r_std),DIMENSION(npts, nlut),INTENT(in)               :: fraclut              !! Fraction of landuse tile
235    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)              :: moist_soil            !! moisture in soil for one pixel (m3/m3)
236    REAL(r_std), DIMENSION(npts,nvm,nslmd,npool,nelements),INTENT(in) :: soilcarbon_input    !! Amount of carbon going into the DOC pools from litter
237                                                                                                                                                                                        !! decomposition \f$(gC m^{-2} day^{-1})$\f     
238
239  !! 0.2 Output variables
240   
241    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: npp_daily            !! Net primary productivity
242                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
243    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out) :: turnover_daily       !! Turnover rates
244                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
245    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: co2_to_bm            !! CO2 taken up from atmosphere when
246                                                                                       !! introducing a new PFT (introduced for
247                                                                                       !! carbon balance closure)
248                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
249    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: co2_fire             !! Carbon emitted into the atmosphere by
250                                                                                       !! fire (living and dead biomass) 
251                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
252    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: resp_hetero          !! Heterotrophic respiration
253                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
254    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: resp_maint           !! Maintenance respiration 
255                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
256    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: resp_growth          !! Growth respiration 
257                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
258   
259    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: deadleaf_cover       !! Fraction of soil covered by dead leaves
260                                                                                       !! (0 to 1, unitless)
261    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: vcmax                !! Maximum rate of carboxylation
262    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out):: bm_to_litter      !! Conversion of biomass to litter
263                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
264    LOGICAL, DIMENSION(npts,nvm), INTENT(out)                  :: begin_leaves         !! signal to start putting leaves on (true/false)
265
266    !! 0.3 Modified variables
267   
268    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: height               !! Height of vegetation (m)
269    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: lai                  !! Leaf area index OF AN INDIVIDUAL PLANT,
270                                                                                       !! where a PFT contains n indentical plants
271                                                                                       !! i.e., using the mean individual approach
272                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
273    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: rprof                !! Prescribed root depth (m)
274    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: PFTpresent           !! Tab indicating which PFTs are present in
275                                                                                       !! each pixel
276    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: age                  !! Age (years)   
277    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: fireindex            !! Probability of fire (0 to 1, unitless)
278    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: firelitter           !! Longer term litter above the ground that
279                                                                                       !! can be burned, @tex $(gC m^{-2})$ @endtex
280    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_age             !! Leaf age (days)
281    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_frac            !! Fraction of leaves in leaf age class,
282                                                                                       !! (0 to 1, unitless)
283    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass        !! Biomass @tex $(gC m^{-2})$ @endtex
284    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: ind                  !! Density of individuals
285                                                                                       !! @tex $(m^{-2})$ @endtex
286    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: adapted              !! Adaptation of PFT (killed if too cold)
287                                                                                       !! (0 to 1, unitless)
288    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: regenerate           !! "Fitness": Winter sufficiently cold for
289                                                                                       !! PFT regeneration ? (0 to 1, unitless)
290    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: senescence           !! Flag for setting senescence stage (only
291                                                                                       !! for deciduous trees)
292    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: when_growthinit      !! How many days ago was the beginning of
293                                                                                       !! the growing season (days)
294    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)      :: litterpart           !! Fraction of litter above the ground
295                                                                                       !! belonging to different PFTs
296                                                                                       !! (0 to 1, unitless)
297    REAL(r_std), DIMENSION(npts,nlitt,nvm, nelements), INTENT(inout):: litter_above    !! Metabolic and structural litter, above ground
298                                                                                       !! @tex $(gC m^{-2})$ @endtex
299    REAL(r_std), DIMENSION(npts,nlitt,nvm,nslmd,nelements), INTENT(inout):: litter_below !! Metabolic and structural litter, below ground
300                                                                                       !! @tex $(gC m^{-2})$ @endtex
301    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)      :: dead_leaves          !! Dead leaves on ground, per PFT, metabolic
302                                                                                       !! and structural, 
303                                                                                       !! @tex $(gC m^{-2})$ @endtex
304    REAL(r_std), DIMENSION(npts,ncarb,nvm,nslmd), INTENT(inout) :: carbon               !! Carbon pool: active, slow, or passive,
305                                                                                       !! @tex $(gC m^{-2})$ @endtex
306    REAL(r_std), DIMENSION(npts,nvm,nslmd,ndoc,npool,nelements), INTENT(inout) :: DOC   !! Dissolved Organic Carbon in soil
307                                                                                       !! The unit is given by m^2 of ground
308                                                                                       !! @tex $(gC m{-2} of ground)$ @endtex
309    REAL(r_std), DIMENSION(npts,nvm,nexp,nflow,nelements), INTENT(in) :: DOC_EXP       !! Exported Dissolved Organic Carbon from soil
310                                                                                       !! The unit is given by m^2 of ground
311                                                                                       !! @tex $(gC m{-2} of ground) $ @endtex
312    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: lignin_struc_above   !! Ratio of Lignin/Carbon in structural
313                                                                                       !! litter, above ground, 
314                                                                                       !! @tex $(gC m^{-2})$ @endtex
315        REAL(r_std),DIMENSION (npts,nvm), INTENT (inout)           :: depth_deepsoil       !! Depth of the soil layer deeper than 2 m.
316                                                                                       !! When sediment deposition occuring, the original surface (0-2)
317                                                                                                                                                                   !! soil DOC, and SOC will enther into this layer.   
318    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: tot_soil_resp        !! Total soil respiration
319                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
320    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: veget_cov_max        !! "Maximal" coverage fraction of a PFT (LAI
321                                                                                       !! -> infinity) on ground
322    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: npp_longterm         !! "Long term" mean yearly primary
323                                                                                       !! productivity
324                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
325    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: lm_lastyearmax       !! Last year's maximum leaf mass, for each
326                                                                                       !! PFT @tex $(gC m^{-2})$ @endtex
327    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: veget_lastlight      !! Vegetation fractions (on ground) after
328                                                                                       !! last light competition 
329                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
330    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: everywhere           !! Is the PFT everywhere in the grid box or
331                                                                                       !! very localized (after its introduction)
332                                                                                       !! (unitless)
333    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: need_adjacent        !! In order for this PFT to be introduced,
334                                                                                       !! does it have to be present in an
335                                                                                       !! adjacent grid box?
336    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: RIP_time             !! How much time ago was the PFT eliminated
337                                                                                       !! for the last time (y)
338    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: turnover_time        !! Turnover_time of leaves for grasses
339                                                                                       !! (days)
340    REAL(r_std),DIMENSION(npts,0:10), INTENT(inout)            :: prod10               !! Products remaining in the 10
341                                                                                       !! year-turnover pool after the annual
342                                                                                       !! release for each compartment (10
343                                                                                       !! + 1 : input from year of land cover
344                                                                                       !! change) @tex $(gC m^{-2})$ @endtex
345    REAL(r_std),DIMENSION(npts,0:100), INTENT(inout)           :: prod100              !! Products remaining in the 100
346                                                                                       !! year-turnover pool after the annual
347                                                                                       !! release for each compartment (100
348                                                                                       !! + 1 : input from year of land cover
349                                                                                       !! change) @tex $(gC m^{-2})$ @endtex
350    REAL(r_std),DIMENSION(npts,10), INTENT(inout)              :: flux10               !! Annual release from the 10
351                                                                                       !! year-turnover pool compartments 
352                                                                                       !! @tex $(gC m^{-2})$ @endtex
353    REAL(r_std),DIMENSION(npts,100), INTENT(inout)             :: flux100              !! Annual release from the 100
354                                                                                       !! year-turnover pool compartments 
355                                                                                       !! @tex $(gC m^{-2})$ @endtex
356    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: convflux             !! Release during first year following land
357                                                                                       !! cover change @tex $(gC m^{-2})$ @endtex
358    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: cflux_prod10         !! Total annual release from the 10
359                                                                                       !! year-turnover pool
360                                                                                       !! @tex $(gC m^{-2})$ @endtex
361    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: cflux_prod100        !! Total annual release from the 100
362                                                                                       !! year-turnover pool
363                                                                                       !! @tex $(gC m^{-2})$ @endtex
364    REAL(r_std),DIMENSION(npts,0:10), INTENT(inout)            :: prod10_harvest       !! Products remaining in the 10
365                                                                                       !! year-turnover pool after the annual
366                                                                                       !! release for each compartment (10
367                                                                                       !! + 1 : input from year of wood harvest)
368                                                                                       !! @tex $(gC m^{-2})$ @endtex
369    REAL(r_std),DIMENSION(npts,0:100), INTENT(inout)           :: prod100_harvest      !! Products remaining in the 100
370                                                                                       !! year-turnover pool after the annual
371                                                                                       !! release for each compartment (100
372                                                                                       !! + 1 : input from year of wood harvest)
373                                                                                       !! @tex $(gC m^{-2})$ @endtex
374    REAL(r_std),DIMENSION(npts,10), INTENT(inout)              :: flux10_harvest       !! Annual release from the 10
375                                                                                       !! year-turnover pool compartments 
376                                                                                       !! @tex $(gC m^{-2})$ @endtex
377    REAL(r_std),DIMENSION(npts,100), INTENT(inout)             :: flux100_harvest      !! Annual release from the 100
378                                                                                       !! year-turnover pool compartments 
379                                                                                       !! @tex $(gC m^{-2})$ @endtex
380    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: convflux_harvest     !! Release during first year following wood
381                                                                                       !! harvest @tex $(gC m^{-2})$ @endtex
382    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: cflux_prod10_harvest !! Total annual release from the 10
383                                                                                       !! year-turnover pool
384                                                                                       !! @tex $(gC m^{-2})$ @endtex
385    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: cflux_prod100_harvest!! Total annual release from the 100
386                                                                                       !! year-turnover pool
387                                                                                       !! @tex $(gC m^{-2})$ @endtex
388    REAL(r_std),DIMENSION(npts,nvm), INTENT(inout)             :: woodharvestpft       !! Harvested wood biomass (gC m-2 dt_stomate-1)
389    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)             :: convfluxpft         !! Convflux per PFT
390    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)             :: fDeforestToProduct  !! Deforested biomass into product pool due to anthropogenic
391                                                                                       !! land use change
392    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)             :: fLulccResidue       !! Carbon mass flux into soil and litter due to anthropogenic land use or land cover change
393    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)             :: fHarvestToProduct   !! Deforested biomass into product pool due to anthropogenic
394                                                                                       !! land use
395
396    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: harvest_above        !! Harvest above ground biomass for
397                                                                                       !! agriculture @tex $(gC m^{-2})$ @endtex
398    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: carb_mass_total      !! Carbon Mass total (soil, litter, veg)
399                                                                                       !! @tex $(gC m^{-2})$ @endtex 
400    REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA         !! Matrix containing the fluxes 
401                                                                                       !! between the carbon pools
402                                                                                       !! per sechiba time step
403                                                                                       !! @tex $(gC.m^2.day^{-1})$ @endtex
404
405    !! 0.4 Local variables
406
407    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_bm_to_litter    !! Total conversion of biomass to litter
408                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
409    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_live_biomass    !! Total living biomass 
410                                                                                       !! @tex $(gC m{-2})$ @endtex
411    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: cOther              !! Carbon Mass in Vegetation Components
412                                                                                       !! other than Leaves, Stems and Roots
413                                                                                       !! @tex $(gC m{-2})$ @endtex
414
415    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements)           :: bm_alloc            !! Biomass increase, i.e. NPP per plant part
416                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
417    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_turnover        !! Total turnover rate 
418                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
419    REAL(r_std), DIMENSION(npts,nvm,nslmd)                      :: tot_litter_soil_carb!! Total soil and litter carbon 
420                                                                                       !! @tex $(gC m^{-2})$ @endtex
421    REAL(r_std), DIMENSION(npts,nvm,nslmd)                      :: tot_litter_carb     !! Total litter carbon
422                                                                                       !! @tex $(gC m^{-2})$ @endtex
423    REAL(r_std), DIMENSION(npts,nvm,nslmd)                      :: tot_soil_carb       !! Total soil carbon 
424                                                                                       !! @tex $(gC m^{-2})$ @endtex
425    REAL(r_std), DIMENSION(npts)                                :: sum_cLitterGrass    !! Carbon mass in litter on grass tiles
426                                                                                       !! @tex $(kgC !m^{-2})$ @endtex
427    REAL(r_std), DIMENSION(npts)                                :: sum_cLitterCrop     !! Carbon mass in litter on crop tiles
428                                                                                       !! @tex $(kgC m^{-2})$ @endtex
429    REAL(r_std), DIMENSION(npts)                                :: sum_cSoilGrass      !! Carbon mass in soil on grass tiles
430                                                                                       !! @tex $(kgC !m^{-2})$ @endtex
431    REAL(r_std), DIMENSION(npts)                                :: sum_cSoilCrop       !! Carbon mass in soil on crop tiles
432                                                                                       !! @tex $(kgC m^{-2})$ @endtex
433    REAL(r_std), DIMENSION(npts)                                :: sum_cVegGrass       !! Carbon mass in vegetation on grass tiles
434                                                                                       !! @tex $(kgC !m^{-2})$ @endtex
435    REAL(r_std), DIMENSION(npts)                                :: sum_cVegCrop        !! Carbon mass in vegetation on crop tiles
436                                                                                       !! @tex $(kgC m^{-2})$ @endtex
437    REAL(r_std), DIMENSION(npts)                                :: sum_cLitterTree     !! Carbon mass in litter on tree tiles
438                                                                                       !! @tex $(kgC !!m^{-2})$ @endtex
439    REAL(r_std), DIMENSION(npts)                                :: sum_cSoilTree       !! Carbon mass in soil on tree tiles
440                                                                                       !! @tex $(kgC !m^{-2})$ @endtex
441    REAL(r_std), DIMENSION(npts)                                :: sum_cVegTree        !! Carbon mass in vegetation on tree tiles
442                                                                                       !! @tex $(kgC !m^{-2})$ @endtex
443    REAL(r_std), DIMENSION(npts)                                :: carb_mass_variation !! Carbon Mass variation 
444                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
445    REAL(r_std), DIMENSION(npts,nvm)                            :: cn_ind              !! Crown area of individuals
446                                                                                       !! @tex $(m^{2})$ @endtex
447    REAL(r_std), DIMENSION(npts,nvm)                            :: woodmass_ind        !! Woodmass of individuals (gC)
448    REAL(r_std), DIMENSION(npts,nvm,nparts)                     :: f_alloc             !! Fraction that goes into plant part
449                                                                                       !! (0 to 1, unitless)
450    REAL(r_std), DIMENSION(npts)                                :: avail_tree          !! Space availability for trees
451                                                                                       !! (0 to 1, unitless)
452    REAL(r_std), DIMENSION(npts)                                :: avail_grass         !! Space availability for grasses
453                                                                                       !! (0 to 1, unitless)
454    INTEGER                                                     :: ig,j,k,l,m
455    REAL(r_std),DIMENSION(npts)                                 :: prod10_total        !! Total products remaining in the pool
456                                                                                       !! after the annual release
457                                                                                       !! @tex $(gC m^{-2})$ @endtex
458    REAL(r_std),DIMENSION(npts)                                 :: prod100_total       !! Total products remaining in the pool
459                                                                                       !! after the annual release
460                                                                                       !! @tex $(gC m^{-2})$ @endtex
461    REAL(r_std),DIMENSION(npts)                                 :: cflux_prod_total    !! Total flux from conflux and the 10/100
462                                                                                       !! year-turnover pool
463                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
464    REAL(r_std),DIMENSION(npts)                                 :: prod10_harvest_total!! Total products remaining in the pool
465                                                                                       !! after the annual release
466                                                                                       !! @tex $(gC m^{-2})$ @endtex
467    REAL(r_std),DIMENSION(npts)                                 :: prod100_harvest_total!! Total products remaining in the pool
468                                                                                       !! after the annual release
469                                                                                       !! @tex $(gC m^{-2})$ @endtex
470    REAL(r_std),DIMENSION(npts)                                 :: cflux_prod_harvest_total!! Total flux from conflux and the 10/100
471                                                                                       !! year-turnover pool
472                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
473    REAL(r_std),DIMENSION(npts,nvm)                             :: veget_cov_max_tmp   !! "Maximal" coverage fraction of a PFT 
474                                                                                       !! (LAI-> infinity) on ground (unitless)
475    REAL(r_std), DIMENSION(npts,nvm)                            :: mortality           !! Fraction of individual dying this time
476                                                                                       !! step (0 to 1, unitless)
477    REAL(r_std), DIMENSION(npts)                                :: vartmp              !! Temporary variable used to add history
478    REAL(r_std), DIMENSION(npts,nslmd)                          :: var_tmp             !! Temporary variable used to add history
479    REAL(r_std), DIMENSION(npts,nvm)                            :: tmp_var             !! Temporary variable used to add history       
480    REAL(r_std), DIMENSION(npts,nvm)                            :: histvar             !! History variables
481    REAL(r_std), DIMENSION(0:nslm)                              :: z_soil              !! Variable to store depth of the different
482                                                                                       !! soil layers (m)
483                                                                                                                                                                           
484    REAL(r_std), DIMENSION(npts,nlut)                           :: clitterlut          !! Litter carbon on landusetype4 (nlut)
485    REAL(r_std), DIMENSION(npts,nlut)                           :: csoillut            !! Soil carbon on landusetype4 (nlut)
486    REAL(r_std), DIMENSION(npts,nlut)                           :: cveglut             !! Carbon in vegetation on landusetype4 (nlut)
487    REAL(r_std), DIMENSION(npts,nlut)                           :: lailut              !! LAI on landusetype4 (nlut)
488    REAL(r_std), DIMENSION(npts,nlut)                           :: ralut               !! Autotrophic respiration on landusetype4 (nlut)
489    REAL(r_std), DIMENSION(npts,nlut)                           :: rhlut               !! Heterotrophic respiration on landusetype4 (nlut)
490    REAL(r_std), DIMENSION(npts,nlut)                           :: npplut              !! Net Primary Productivity on landusetype4 (nlut)   
491    REAL(r_std), DIMENSION(npts,nlut)                           :: ctotfirelut         !! Fire CO2 emission on landusetype4 (nlut)
492    REAL(r_std), DIMENSION(npts,nlut)                           :: cproductlut
493    REAL(r_std), DIMENSION(npts,nlut)                           :: flulccatmlut
494    REAL(r_std), DIMENSION(npts,nlut)                           :: flulccproductlut
495    REAL(r_std), DIMENSION(npts,nlut)                           :: flulccresiduelut
496    REAL(r_std), DIMENSION(npts,ncarb)                          :: csoilpools          !! Diagnostics for carbon in soil pools
497!_ ================================================================================================================================
498
499    IF (printlev>=3) WRITE(numout,*) 'Entering stomate_lpj'
500
501 
502  !! 1. Initializations
503    !WRITE(numout,*) 'STOMATE_LPJ1'
504        !WRITE(numout,*) 'litter_above: ', MAXVAL(litter_above(:,:,:,icarbon))
505        !WRITE(numout,*) 'litter_below: ', MAXVAL(litter_below(:,:,:,:,icarbon))
506        !WRITE(numout,*) 'SOC_max',MAXVAL(carbon),MAXVAL(DOC)
507        !WRITE(numout,*) 'SOC_min',MINVAL(carbon),MINVAL(DOC)
508   
509    !! 1.1 Initialize variables to zero
510    co2_to_bm(:,:) = zero
511    co2_fire(:,:) = zero
512    npp_daily(:,:) = zero
513    resp_maint(:,:) = zero
514    resp_growth(:,:) = zero
515    harvest_above(:) = zero
516    bm_to_litter(:,:,:,:) = zero
517    cn_ind(:,:) = zero
518    woodmass_ind(:,:) = zero
519    turnover_daily(:,:,:,:) = zero
520       
521        z_soil(0) = zero
522    z_soil(1:nslm) = zlt(1:nslm)
523   
524    !! 1.2  Initialize variables to veget_cov_max
525    veget_cov_max_tmp(:,:) = veget_cov_max(:,:)
526
527    !! 1.3 Calculate some vegetation characteristics
528   
529    !! 1.3.1 Calculate some vegetation characteristics
530    !        Calculate cn_ind (individual crown mass) and individual height from
531    !        state variables if running DGVM or dynamic mortality in static cover mode
532    !??        Explain (maybe in the header once) why you mulitply with veget_cov_max in the DGVM
533    !??        and why you don't multiply with veget_cov_max in stomate.
534    IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
535       IF(ok_dgvm) THEN
536          WHERE (ind(:,:).GT.min_stomate)
537             woodmass_ind(:,:) = &
538                  ((biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
539                  +biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)) & 
540                  *veget_cov_max(:,:))/ind(:,:)
541          ENDWHERE
542       ELSE
543          WHERE (ind(:,:).GT.min_stomate)
544             woodmass_ind(:,:) = &
545                  (biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
546                  +biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon))/ind(:,:)
547          ENDWHERE
548       ENDIF
549
550       CALL crown (npts,  PFTpresent, &
551            ind, biomass, woodmass_ind, &
552            veget_cov_max, cn_ind, height)
553    ENDIF
554
555    !! 1.3.2 Prescribe characteristics if the vegetation is not dynamic
556    !        IF the DGVM is not activated, the density of individuals and their crown
557    !        areas don't matter, but they should be defined for the case we switch on
558    !        the DGVM afterwards. At the first call, if the DGVM is not activated,
559    !        impose a minimum biomass for prescribed PFTs and declare them present.
560    CALL prescribe (npts, &
561         veget_cov_max, dt_days, PFTpresent, everywhere, when_growthinit, &
562         biomass, leaf_frac, ind, cn_ind, co2_to_bm)
563
564
565  !! 2. Climatic constraints for PFT presence and regenerativeness
566
567    !   Call this even when DGVM is not activated so that "adapted" and "regenerate"
568    !   are kept up to date for the moment when the DGVM is activated.
569    CALL constraints (npts, dt_days, &
570         t2m_month, t2m_min_daily,when_growthinit, Tseason, &
571         adapted, regenerate)
572
573   
574  !! 3. Determine introduction and elimination of PTS based on climate criteria
575 
576    IF ( ok_dgvm ) THEN
577     
578       !! 3.1 Calculate introduction and elimination
579       CALL pftinout (npts, dt_days, adapted, regenerate, &
580            neighbours, veget_cov_max, &
581            biomass, ind, cn_ind, age, leaf_frac, npp_longterm, lm_lastyearmax, senescence, &
582            PFTpresent, everywhere, when_growthinit, need_adjacent, RIP_time, &
583            co2_to_bm, &
584            avail_tree, avail_grass)
585
586       !! 3.2 Reset attributes for eliminated PFTs.
587       !     This also kills PFTs that had 0 leafmass during the last year. The message
588       !     "... after pftinout" is misleading in this case.
589       CALL kill (npts, 'pftinout  ', lm_lastyearmax, &
590            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
591            lai, age, leaf_age, leaf_frac, npp_longterm, &
592            when_growthinit, everywhere, veget_cov_max, bm_to_litter)
593
594       
595       !! 3.3 Calculate woodmass of individual tree
596       IF(ok_dgvm) THEN
597          WHERE ((ind(:,:).GT.min_stomate))
598             woodmass_ind(:,:) = &
599                  ((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
600                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))*veget_cov_max(:,:))/ind(:,:)
601          ENDWHERE
602       ELSE
603          WHERE ((ind(:,:).GT.min_stomate))
604             woodmass_ind(:,:) =(biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
605                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
606          ENDWHERE
607       ENDIF
608       
609       ! Calculate crown area and diameter for all PFTs (including the newly established)
610       CALL crown (npts, PFTpresent, &
611            ind, biomass, woodmass_ind, &
612            veget_cov_max, cn_ind, height)
613
614    ENDIF
615   
616  !! 4. Phenology
617
618    !! 4.1 Write values to history file
619    !      Current values for ::when_growthinit
620    CALL xios_orchidee_send_field("WHEN_GROWTHINIT",when_growthinit)
621
622    CALL histwrite_p (hist_id_stomate, 'WHEN_GROWTHINIT', itime, when_growthinit, npts*nvm, horipft_index)
623
624    ! Set and write values for ::PFTpresent
625    WHERE(PFTpresent)
626       histvar=un
627    ELSEWHERE
628       histvar=zero
629    ENDWHERE
630
631    CALL xios_orchidee_send_field("PFTPRESENT",histvar)
632
633    CALL histwrite_p (hist_id_stomate, 'PFTPRESENT', itime, histvar, npts*nvm, horipft_index)
634
635    ! Set and write values for gdd_midwinter
636    WHERE(gdd_midwinter.EQ.undef)
637       histvar=val_exp
638    ELSEWHERE
639       histvar=gdd_midwinter
640    ENDWHERE
641
642    CALL xios_orchidee_send_field("GDD_MIDWINTER",histvar)
643
644    CALL histwrite_p (hist_id_stomate, 'GDD_MIDWINTER', itime, histvar, npts*nvm, horipft_index)
645
646    ! Set and write values for gdd_m5_dormance
647    WHERE(gdd_m5_dormance.EQ.undef)
648       histvar=val_exp
649    ELSEWHERE
650       histvar=gdd_m5_dormance
651    ENDWHERE
652   
653    CALL xios_orchidee_send_field('GDD_M5_DORMANCE',histvar)
654    CALL histwrite_p (hist_id_stomate, 'GDD_M5_DORMANCE', itime, histvar, npts*nvm, horipft_index)
655
656    ! Set and write values for ncd_dormance
657    WHERE(ncd_dormance.EQ.undef)
658       histvar=val_exp
659    ELSEWHERE
660       histvar=ncd_dormance
661    ENDWHERE
662
663    CALL xios_orchidee_send_field("NCD_DORMANCE",histvar)
664
665    CALL histwrite_p (hist_id_stomate, 'NCD_DORMANCE', itime, histvar, npts*nvm, horipft_index)
666
667    !! 4.2 Calculate phenology
668    CALL phenology (npts, dt_days, PFTpresent, &
669         veget_cov_max, &
670         t2m_longterm, t2m_month, t2m_week, gpp_daily, &
671         maxmoiavail_lastyear, minmoiavail_lastyear, &
672         moiavail_month, moiavail_week, &
673         gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
674         senescence, time_hum_min, &
675         biomass, leaf_frac, leaf_age, &
676         when_growthinit, co2_to_bm, &
677         begin_leaves)
678   
679  !! 5. Allocate C to different plant parts
680   
681    CALL alloc (npts, dt_days, &
682         lai, veget_cov_max, senescence, when_growthinit, &
683         moiavail_week, tsoil_month, soilhum_month, &
684         biomass, age, leaf_age, leaf_frac, rprof, f_alloc)
685
686  !! 6. NPP, maintenance and growth respiration
687
688    !! 6.1 Calculate NPP and respiration terms
689    CALL npp_calc (npts, dt_days, &
690         PFTpresent, &
691         t2m_daily, tsoil_daily, lai, rprof, &
692         gpp_daily, f_alloc, bm_alloc, resp_maint_part,&
693         biomass, leaf_age, leaf_frac, age, &
694         resp_maint, resp_growth, npp_daily, co2_to_bm)
695
696    !! 6.2 Kill slow growing PFTs in DGVM or STOMATE with constant mortality
697    IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
698       CALL kill (npts, 'npp       ', lm_lastyearmax,  &
699            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
700            lai, age, leaf_age, leaf_frac, npp_longterm, &
701            when_growthinit, everywhere, veget_cov_max, bm_to_litter)
702
703       !! 6.2.1 Update wood biomass     
704       !        For the DGVM
705       IF(ok_dgvm) THEN
706          WHERE (ind(:,:).GT.min_stomate)
707             woodmass_ind(:,:) = &
708                  ((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
709                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon)) & 
710                  *veget_cov_max(:,:))/ind(:,:)
711          ENDWHERE
712
713       ! For all pixels with individuals
714       ELSE
715          WHERE (ind(:,:).GT.min_stomate)
716             woodmass_ind(:,:) = &
717                  (biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
718                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
719          ENDWHERE
720       ENDIF ! ok_dgvm
721
722       !! 6.2.2 New crown area and maximum vegetation cover after growth
723       CALL crown (npts, PFTpresent, &
724            ind, biomass, woodmass_ind,&
725            veget_cov_max, cn_ind, height)
726
727    ENDIF ! ok_dgvm
728   
729  !! 7. fire
730
731    !! 7.1. Burn PFTs
732    CALL fire (npts, dt_days, litterpart, &
733         litterhum_daily, t2m_daily, lignin_struc_above, veget_cov_max, &
734         fireindex, firelitter, biomass, ind, &
735         litter_above, dead_leaves, bm_to_litter, &
736         co2_fire, MatrixA)
737
738    !! 7.2 Kill PFTs in DGVM
739    IF ( ok_dgvm ) THEN
740
741       ! reset attributes for eliminated PFTs
742       CALL kill (npts, 'fire      ', lm_lastyearmax, &
743            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
744            lai, age, leaf_age, leaf_frac, npp_longterm, &
745            when_growthinit, everywhere, veget_cov_max, bm_to_litter)
746
747    ENDIF ! ok_dgvm
748 
749  !! 8. Tree mortality
750
751    ! Does not depend on age, therefore does not change crown area.
752    CALL gap (npts, dt_days, &
753         npp_longterm, turnover_longterm, lm_lastyearmax, &
754         PFTpresent, t2m_min_daily, Tmin_spring_time, &
755         biomass, ind, bm_to_litter, mortality)
756
757
758    IF ( ok_dgvm ) THEN
759
760       ! reset attributes for eliminated PFTs
761       CALL kill (npts, 'gap       ', lm_lastyearmax, &
762            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
763            lai, age, leaf_age, leaf_frac, npp_longterm, &
764            when_growthinit, everywhere, veget_cov_max, bm_to_litter)
765
766    ENDIF
767
768  !! 9. Leaf senescence, new lai and other turnover processes
769
770    CALL turn (npts, dt_days, PFTpresent, &
771         herbivores, &
772         maxmoiavail_lastyear, minmoiavail_lastyear, &
773         moiavail_week,  moiavail_month,t2m_longterm, t2m_month, t2m_week, veget_cov_max, &
774         gdd_from_growthinit, leaf_age, leaf_frac, age, lai, biomass, &
775         turnover_daily, senescence,turnover_time)
776
777    !! 10. Light competition
778   
779    !! If not using constant mortality then kill with light competition
780!    IF ( ok_dgvm .OR. .NOT.(lpj_gap_const_mort) ) THEN
781    IF ( ok_dgvm ) THEN
782 
783       !! 10.1 Light competition
784       CALL light (npts, dt_days, &
785            veget_cov_max, fpc_max, PFTpresent, cn_ind, lai, maxfpc_lastyear, &
786            lm_lastyearmax, ind, biomass, veget_lastlight, bm_to_litter, mortality)
787       
788       !! 10.2 Reset attributes for eliminated PFTs
789       CALL kill (npts, 'light     ', lm_lastyearmax, &
790            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
791            lai, age, leaf_age, leaf_frac, npp_longterm, &
792            when_growthinit, everywhere, veget_cov_max, bm_to_litter)
793    ENDIF
794
795   
796  !! 11. Establishment of saplings
797   
798    IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort ) THEN
799
800       !! 11.1 Establish new plants
801       CALL establish (npts, dt_days, PFTpresent, regenerate, &
802            neighbours, resolution, need_adjacent, herbivores, &
803            precip_lastyear, gdd0_lastyear, lm_lastyearmax, &
804            cn_ind, lai, avail_tree, avail_grass, npp_longterm, &
805            leaf_age, leaf_frac, &
806            ind, biomass, age, everywhere, co2_to_bm, veget_cov_max, woodmass_ind, &
807            mortality, bm_to_litter)
808
809       !! 11.2 Calculate new crown area (and maximum vegetation cover)
810       CALL crown (npts, PFTpresent, &
811            ind, biomass, woodmass_ind, &
812            veget_cov_max, cn_ind, height)
813
814    ENDIF
815
816  !! 12. Calculate final LAI and vegetation cover
817               
818    CALL cover (npts, cn_ind, ind, biomass, &
819         veget_cov_max, veget_cov_max_tmp, lai, &
820         litter_above, litter_below, carbon, DOC, turnover_daily, bm_to_litter, &
821         co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily)
822               
823  !! 13. Update litter pools to account for harvest
824    !WRITE(numout,*) 'STOMATE_LPJ2'
825        !WRITE(numout,*) 'litter_above: ', MAXVAL(litter_above(:,:,:,icarbon))
826        !WRITE(numout,*) 'litter_below: ', MAXVAL(litter_below(:,:,:,:,icarbon))
827        !WRITE(numout,*) 'SOC_max',MAXVAL(carbon),MAXVAL(DOC)
828        !WRITE(numout,*) 'SOC_min',MINVAL(carbon),MINVAL(DOC)
829    ! the whole litter stuff:
830    !    litter update, lignin content, PFT parts, litter decay,
831    !    litter heterotrophic respiration, dead leaf soil cover.
832    !    No vertical discretisation in the soil for litter decay.\n
833    ! added by shilong for harvest
834    IF(harvest_agri) THEN
835       CALL harvest(npts, dt_days, veget_cov_max, &
836            bm_to_litter, turnover_daily, &
837            harvest_above)
838    ENDIF
839
840    !! 14. Land cover change if it is time to do so
841    !! The flag do_now_lcchange is set in slowproc_main at the same time as the vegetation is read from file.
842    !! The vegetation fractions are not updated yet and will be updated in the end of sechiba_main.
843    IF (do_now_stomate_woodharvest) THEN
844       CALL woodharvest_main(npts, dt_days, veget_cov_max, &
845            biomass, &
846            flux10_harvest,flux100_harvest, prod10_harvest,prod100_harvest,&
847            convflux_harvest,cflux_prod10_harvest,cflux_prod100_harvest,&
848            woodharvest,woodharvestpft,fHarvestToProduct)
849       do_now_stomate_woodharvest=.FALSE.
850    ENDIF
851
852
853    IF (do_now_stomate_lcchange) THEN
854       CALL lcchange_main (npts, dt_days, veget_cov_max, veget_cov_max_new, &
855            biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
856            co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,flux10,flux100, &
857            prod10,prod100,convflux,cflux_prod10,cflux_prod100,leaf_frac,&
858            npp_longterm, lm_lastyearmax, litter_above, litter_below, carbon, DOC, &
859            depth_deepsoil,convfluxpft,  fDeforestToProduct, fLulccResidue)
860       do_now_stomate_lcchange=.FALSE.
861
862       ! Set the flag done_stomate_lcchange to be used in the end of sechiba_main to update the fractions.
863       done_stomate_lcchange=.TRUE.
864    ENDIF
865    !WRITE(numout,*) 'STOMATE_LPJ3'
866        !WRITE(numout,*) 'litter_above: ', MAXVAL(litter_above(:,:,:,icarbon))
867        !WRITE(numout,*) 'litter_below: ', MAXVAL(litter_below(:,:,:,:,icarbon))
868        !WRITE(numout,*) 'SOC_max',MAXVAL(carbon),MAXVAL(DOC)
869        !WRITE(numout,*) 'SOC_min',MINVAL(carbon),MINVAL(DOC)
870    !! 15. Calculate vcmax
871
872    CALL vmax (npts, dt_days, &
873         leaf_age, leaf_frac, &
874         vcmax)
875
876    !MM déplacement pour initialisation correcte des grandeurs cumulées :
877    cflux_prod_total(:) = convflux(:) + cflux_prod10(:) + cflux_prod100(:)
878    prod10_total(:)=SUM(prod10,dim=2)
879    prod100_total(:)=SUM(prod100,dim=2)
880
881    cflux_prod_harvest_total(:) = convflux_harvest(:) + cflux_prod10_harvest(:) + cflux_prod100_harvest(:)
882    prod10_harvest_total(:)=SUM(prod10_harvest,dim=2)
883    prod100_harvest_total(:)=SUM(prod100_harvest,dim=2)
884   
885  !! 16. Total heterotrophic respiration
886
887    tot_soil_carb(:,:,:) = zero
888    tot_litter_carb(:,:,:) = zero
889    sum_cLitterGrass = zero
890    sum_cLitterCrop = zero
891    sum_cSoilGrass = zero
892    sum_cSoilCrop = zero
893    sum_cVegGrass = zero
894    sum_cVegCrop = zero
895    sum_cLitterTree = zero
896    sum_cSoilTree = zero
897    sum_cVegTree = zero
898       
899        DO j=2,nvm
900          !tot_litter_carb(:,j,1) = tot_litter_carb(:,j,1) + &
901          !     &          litter_above(:,istructural,j,icarbon) + litter_above(:,imetabolic,j,icarbon)
902       DO l=1,nslmd
903          tot_litter_carb(:,j,l) = tot_litter_carb(:,j,l) + &
904               &          litter_below(:,istructural,j,l,icarbon) + litter_below(:,imetabolic,j,l,icarbon)
905       ENDDO
906
907       DO l=1,nslmd
908
909          tot_soil_carb(:,j,l) = tot_soil_carb(:,j,l) + (carbon(:,iactive,j,l) + &
910               &          carbon(:,islow,j,l) + carbon(:,ipassive,j,l)) 
911       ENDDO
912
913      DO k =1,npool
914        DO l=1,nslmd
915
916           tot_soil_carb(:,j,l) = tot_soil_carb(:,j,l) + &
917                &          (DOC(:,j,l,ifree,k,icarbon) + DOC(:,j,l,iadsorbed,k,icarbon))
918        ENDDO
919
920      ENDDO
921
922    ENDDO 
923
924    tot_litter_soil_carb(:,:,:) = tot_litter_carb(:,:,:) + tot_soil_carb(:,:,:)
925       
926        WHERE (ISNAN(tot_litter_soil_carb(:,:,:)))
927                tot_litter_soil_carb(:,:,:)=zero
928        ENDWHERE
929
930 !   DO j=2,nvm
931
932!      IF ((.NOT. is_tree(j))  .AND. natural(j)) THEN
933!         sum_cLitterGrass(:) = sum_cLitterGrass(:) + tot_litter_carb(:,j)*veget_cov_max(:,j)
934!      ELSE IF ((.NOT. is_tree(j))  .AND. (.NOT. natural(j)) ) THEN
935!         sum_cLitterCrop(:) = sum_cLitterCrop(:) + tot_litter_carb(:,j)*veget_cov_max(:,j)
936!      ELSE IF (is_tree(j)) THEN
937!         sum_cLitterTree(:) = sum_cLitterTree(:) + tot_litter_carb(:,j)*veget_cov_max(:,j)
938!      ENDIF
939!
940!      tot_soil_carb(:,j) = tot_soil_carb(:,j) + (carbon(:,iactive,j) + &
941!           &          carbon(:,islow,j)+  carbon(:,ipassive,j))
942
943!      IF ((.NOT. is_tree(j)) .AND. natural(j)) THEN
944!         sum_cSoilGrass(:) = sum_cSoilGrass(:) + tot_soil_carb(:,j)*veget_cov_max(:,j)
945!      ELSE IF ((.NOT. is_tree(j))  .AND. (.NOT. natural(j)) ) THEN
946!         sum_cSoilCrop(:) = sum_cSoilCrop(:) + tot_soil_carb(:,j)*veget_cov_max(:,j)
947!      ELSE IF (is_tree(j)) THEN
948!         sum_cSoilTree(:) = sum_cSoilTree(:) + tot_soil_carb(:,j)*veget_cov_max(:,j)
949!      END IF
950!    ENDDO
951
952!!$     DO k = 1, nelements ! Loop over # elements
953!!$        tot_live_biomass(:,:,k) = biomass(:,:,ileaf,k) + biomass(:,:,isapabove,k) + biomass(:,:,isapbelow,k) +&
954!!$             &                    biomass(:,:,iheartabove,k) + biomass(:,:,iheartbelow,k) + &
955!!$             &                    biomass(:,:,iroot,k) + biomass(:,:,ifruit,k) + biomass(:,:,icarbres,k)
956!!$    END DO ! Loop over # elements
957
958    tot_live_biomass(:,:,:) = biomass(:,:,ileaf,:) + biomass(:,:,isapabove,:) + biomass(:,:,isapbelow,:) +&
959             &                    biomass(:,:,iheartabove,:) + biomass(:,:,iheartbelow,:) + &
960             &                    biomass(:,:,iroot,:) + biomass(:,:,ifruit,:) + biomass(:,:,icarbres,:)
961
962    DO j= 1,nvm
963       cOther(:,j,:) = biomass(:,j,ifruit,:) + biomass(:,j,icarbres,:)
964
965       IF ((.NOT. is_tree(j))  .AND. natural(j)) THEN
966          sum_cVegGrass(:) = sum_cVegGrass(:) + tot_live_biomass(:,j,icarbon)*veget_cov_max(:,j)
967       ELSE IF ((.NOT. is_tree(j))  .AND. (.NOT. natural(j)) ) THEN
968          sum_cVegCrop(:) = sum_cVegCrop(:) + tot_live_biomass(:,j,icarbon)*veget_cov_max(:,j)
969       ELSE IF (is_tree(j)) THEN
970          sum_cVegTree(:) = sum_cVegTree(:) + tot_live_biomass(:,j,icarbon)*veget_cov_max(:,j)
971       ENDIF
972    END DO
973
974
975    tot_turnover(:,:,:) = turnover_daily(:,:,ileaf,:) + turnover_daily(:,:,isapabove,:) + &
976         &         turnover_daily(:,:,isapbelow,:) + turnover_daily(:,:,iheartabove,:) + &
977         &         turnover_daily(:,:,iheartbelow,:) + turnover_daily(:,:,iroot,:) + &
978         &         turnover_daily(:,:,ifruit,:) + turnover_daily(:,:,icarbres,:)
979
980    tot_bm_to_litter(:,:,:) = bm_to_litter(:,:,ileaf,:) + bm_to_litter(:,:,isapabove,:) +&
981         &             bm_to_litter(:,:,isapbelow,:) + bm_to_litter(:,:,iheartbelow,:) +&
982         &             bm_to_litter(:,:,iheartabove,:) + bm_to_litter(:,:,iroot,:) + &
983         &             bm_to_litter(:,:,ifruit,:) + bm_to_litter(:,:,icarbres,:)
984
985    carb_mass_variation(:)=-carb_mass_total(:) 
986    DO k=1,nvm
987           carb_mass_total(:)=carb_mass_total(:)+ tot_live_biomass(:,k,icarbon)*veget_cov_max(:,k)+ &
988           &                    (prod10_total + prod100_total) +  (prod10_harvest_total + prod100_harvest_total)
989       DO l= 1,nslmd
990         carb_mass_total(:)=carb_mass_total(:)+ &
991              &                 (tot_litter_carb(:,k,l)+tot_soil_carb(:,k,l))*veget_cov_max(:,k) ! H. Zhang: should not multipied by *(z_soil(l)-z_soil(l-1))
992       ENDDO
993           !carb_mass_total(:)=carb_mass_total(:)+(tot_litter_carb(:,k,nslmd)+tot_soil_carb(:,k,nslmd))*veget_cov_max(:,k)
994    ENDDO               
995    carb_mass_variation(:)=carb_mass_total(:)+carb_mass_variation(:)
996   
997  !! 17. Write history
998       
999    CALL xios_orchidee_send_field("RESOLUTION_X",resolution(:,1))
1000    CALL xios_orchidee_send_field("RESOLUTION_Y",resolution(:,2))
1001    CALL xios_orchidee_send_field("CONTFRAC_STOMATE",contfrac(:))
1002    CALL xios_orchidee_send_field("T2M_MONTH",t2m_month)
1003    CALL xios_orchidee_send_field("T2M_WEEK",t2m_week)
1004    CALL xios_orchidee_send_field("TSEASON",Tseason)
1005    CALL xios_orchidee_send_field("TMIN_SPRING_TIME", Tmin_spring_time)
1006    CALL xios_orchidee_send_field("ONSET_DATE",onset_date)
1007    CALL xios_orchidee_send_field("FPC_MAX",fpc_max)
1008    CALL xios_orchidee_send_field("MAXFPC_LASTYEAR",maxfpc_lastyear)
1009    CALL xios_orchidee_send_field("HET_RESP",resp_hetero(:,:))
1010    CALL xios_orchidee_send_field("TOT_SOIL_RESP",tot_soil_resp(:,:))   
1011    CALL xios_orchidee_send_field("CO2_FIRE",co2_fire)
1012    CALL xios_orchidee_send_field("CO2_TAKEN",co2_to_bm)
1013    CALL xios_orchidee_send_field("LAI",lai)
1014    CALL xios_orchidee_send_field("VEGET_COV_MAX",veget_cov_max)
1015    CALL xios_orchidee_send_field("NPP_STOMATE",npp_daily)
1016    CALL xios_orchidee_send_field("GPP",gpp_daily)
1017    CALL xios_orchidee_send_field("IND",ind)
1018    CALL xios_orchidee_send_field("CN_IND",cn_ind)
1019    CALL xios_orchidee_send_field("WOODMASS_IND",woodmass_ind)
1020    CALL xios_orchidee_send_field("TOTAL_M",tot_live_biomass)
1021    CALL xios_orchidee_send_field("MOISTRESS",moiavail_week)
1022    CALL xios_orchidee_send_field("LEAF_M",biomass(:,:,ileaf,icarbon))
1023    CALL xios_orchidee_send_field("SAP_M_AB",biomass(:,:,isapabove,icarbon))
1024    CALL xios_orchidee_send_field("SAP_M_BE",biomass(:,:,isapbelow,icarbon))
1025    CALL xios_orchidee_send_field("HEART_M_AB",biomass(:,:,iheartabove,icarbon))
1026    CALL xios_orchidee_send_field("HEART_M_BE",biomass(:,:,iheartbelow,icarbon))
1027    CALL xios_orchidee_send_field("ROOT_M",biomass(:,:,iroot,icarbon))
1028    CALL xios_orchidee_send_field("FRUIT_M",biomass(:,:,ifruit,icarbon))
1029    CALL xios_orchidee_send_field("RESERVE_M",biomass(:,:,icarbres,icarbon))
1030    CALL xios_orchidee_send_field("TOTAL_TURN",tot_turnover)
1031    CALL xios_orchidee_send_field("LEAF_TURN",turnover_daily(:,:,ileaf,icarbon))
1032    CALL xios_orchidee_send_field("MAINT_RESP",resp_maint)
1033    CALL xios_orchidee_send_field("GROWTH_RESP",resp_growth)
1034    CALL xios_orchidee_send_field("SAP_AB_TURN",turnover_daily(:,:,isapabove,icarbon))
1035    CALL xios_orchidee_send_field("ROOT_TURN",turnover_daily(:,:,iroot,icarbon))
1036    CALL xios_orchidee_send_field("FRUIT_TURN",turnover_daily(:,:,ifruit,icarbon))
1037    CALL xios_orchidee_send_field("TOTAL_BM_LITTER",tot_bm_to_litter(:,:,icarbon))
1038    CALL xios_orchidee_send_field("LEAF_BM_LITTER",bm_to_litter(:,:,ileaf,icarbon))
1039    CALL xios_orchidee_send_field("SAP_AB_BM_LITTER",bm_to_litter(:,:,isapabove,icarbon))
1040    CALL xios_orchidee_send_field("SAP_BE_BM_LITTER",bm_to_litter(:,:,isapbelow,icarbon))
1041    CALL xios_orchidee_send_field("HEART_AB_BM_LITTER",bm_to_litter(:,:,iheartabove,icarbon))
1042    CALL xios_orchidee_send_field("HEART_BE_BM_LITTER",bm_to_litter(:,:,iheartbelow,icarbon))
1043    CALL xios_orchidee_send_field("ROOT_BM_LITTER",bm_to_litter(:,:,iroot,icarbon))
1044    CALL xios_orchidee_send_field("FRUIT_BM_LITTER",bm_to_litter(:,:,ifruit,icarbon))
1045    CALL xios_orchidee_send_field("RESERVE_BM_LITTER",bm_to_litter(:,:,icarbres,icarbon))
1046    CALL xios_orchidee_send_field("LITTER_STR_AB",litter_above(:,istructural,:,icarbon))
1047    CALL xios_orchidee_send_field("LITTER_MET_AB",litter_above(:,imetabolic,:,icarbon))
1048    CALL xios_orchidee_send_field("DEADLEAF_COVER",deadleaf_cover)
1049       
1050 var_tmp(:,:)=zero
1051 DO l= 1,nvm
1052   var_tmp(:,1)=var_tmp(:,1)+litter_below(:,istructural,l,1,icarbon)*veget_cov_max(:,l)/z_soil(1)
1053   WHERE (depth_deepsoil(:,l).GT.min_sechiba)
1054       var_tmp(:,nslmd)=var_tmp(:,nslmd)+litter_below(:,istructural,l,nslmd,icarbon)*veget_cov_max(:,l)/depth_deepsoil(:,l)
1055   ENDWHERE
1056 ENDDO
1057 DO j= 2,nslm
1058   DO l= 1,nvm
1059    var_tmp(:,j)=var_tmp(:,j)+litter_below(:,istructural,l,j,icarbon)*veget_cov_max(:,l)/(z_soil(j)-z_soil(j-1))
1060   ENDDO
1061 ENDDO
1062 CALL xios_orchidee_send_field("LITTER_STR_BE",var_tmp(:,:))
1063       
1064 var_tmp(:,:)=zero
1065 DO l= 1,nvm
1066   var_tmp(:,1)=var_tmp(:,1)+litter_below(:,imetabolic,l,1,icarbon)*veget_cov_max(:,l)/z_soil(1)
1067   WHERE (depth_deepsoil(:,l).GT.min_sechiba)
1068       var_tmp(:,nslmd)=var_tmp(:,nslmd)+litter_below(:,imetabolic,l,nslmd,icarbon)*veget_cov_max(:,l)/depth_deepsoil(:,l)
1069   ENDWHERE
1070 ENDDO
1071 DO j=2,nslm
1072   DO l= 1,nvm
1073    var_tmp(:,j)=var_tmp(:,j)+litter_below(:,imetabolic,l,j,icarbon)*veget_cov_max(:,l)/(z_soil(j)-z_soil(j-1))
1074   ENDDO
1075 ENDDO
1076 CALL xios_orchidee_send_field("LITTER_MET_BE",var_tmp(:,:))
1077 
1078 var_tmp(:,:)=zero
1079 DO l= 1,nvm
1080   var_tmp(:,1)=var_tmp(:,1)+tot_litter_soil_carb(:,l,1)*veget_cov_max(:,l)/z_soil(1)
1081   WHERE (depth_deepsoil(:,l).GT.min_sechiba)
1082       var_tmp(:,nslmd)=var_tmp(:,nslmd)+tot_litter_soil_carb(:,l,nslmd)*veget_cov_max(:,l)/depth_deepsoil(:,l)
1083   ENDWHERE
1084 ENDDO
1085 
1086 DO j=2,nslm
1087   DO l= 1,nvm
1088      var_tmp(:,j)=var_tmp(:,j)+tot_litter_soil_carb(:,l,j)*veget_cov_max(:,l)/(z_soil(j)-z_soil(j-1))
1089   ENDDO
1090 ENDDO
1091 
1092 !WRITE(numout,*) 'STOMATE_LPJ0'
1093 !WRITE(numout,*) 'litter_above: ', MAXVAL(litter_above),MINVAL(litter_above)
1094 !WRITE(numout,*) 'litter_below: ', MAXVAL(litter_below),MINVAL(litter_below)
1095 !WRITE(numout,*) 'SOC',MAXVAL(carbon),MINVAL(carbon)
1096 !WRITE(numout,*) 'DOC',MAXVAL(DOC),MINVAL(DOC)
1097 !WRITE(numout,*) 'Ctot',MAXVAL(var_tmp),MINVAL(var_tmp)
1098 !DO ig=1,npts
1099!       DO j=1,nslmd
1100!               IF (var_tmp(ig,j).GT.1.0e+10 .OR.var_tmp(ig,j).LT.zero .OR. ISNAN(var_tmp(ig,j))) THEN
1101!!                      WRITE(numout,*) 'Ctotig: ', ig,j,var_tmp(ig,j),MAXVAL(carbon),MAXVAL(DOC),MAXVAL(litter_below),MAXVAL(litter_above)
1102!                       WRITE(numout,*) 'depth_deepsoil: ', depth_deepsoil(ig,:)
1103!                       WRITE(numout,*) 'tot_litter_carb: ',tot_litter_carb(ig,:,j)
1104!                       WRITE(numout,*) 'tot_soil_carb: ',tot_soil_carb(ig,:,j)
1105!                       WRITE(numout,*) 'carbon: ', carbon(ig,1,:,nslmd)
1106!                       WRITE(numout,*) 'carbon: ', carbon(ig,2,:,nslmd)
1107!                       WRITE(numout,*) 'carbon: ', carbon(ig,3,:,nslmd)
1108!               ENDIF
1109!       ENDDO
1110! ENDDO
1111 
1112 CALL xios_orchidee_send_field("TOTAL_SOIL_CARB",var_tmp(:,:))
1113       
1114 var_tmp(:,:)=zero
1115 DO l= 1,nvm
1116    var_tmp(:,1)=var_tmp(:,1)+carbon(:,iactive,l,1)*veget_cov_max(:,l)/z_soil(1)
1117        WHERE (depth_deepsoil(:,l).GT.min_sechiba)
1118       var_tmp(:,nslmd)=var_tmp(:,nslmd)+carbon(:,iactive,l,nslmd)*veget_cov_max(:,l)/depth_deepsoil(:,l)
1119    ENDWHERE
1120 ENDDO
1121 DO j=2,nslm
1122   DO l= 1,nvm
1123    var_tmp(:,j)=var_tmp(:,j)+carbon(:,iactive,l,j)*veget_cov_max(:,l)/(z_soil(j)-z_soil(j-1))
1124   ENDDO
1125 ENDDO
1126 CALL xios_orchidee_send_field("CARBON_ACTIVE",var_tmp(:,:))
1127 
1128 var_tmp(:,:)=zero
1129 DO l= 1,nvm
1130    var_tmp(:,1)=var_tmp(:,1)+carbon(:,islow,l,1)*veget_cov_max(:,l)/z_soil(1)
1131        WHERE (depth_deepsoil(:,l).GT.min_sechiba)
1132       var_tmp(:,nslmd)=var_tmp(:,nslmd)+carbon(:,islow,l,nslmd)*veget_cov_max(:,l)/depth_deepsoil(:,l)
1133    ENDWHERE
1134 ENDDO
1135 DO j= 2,nslm
1136   DO l=1,nvm
1137    var_tmp(:,j)=var_tmp(:,j)+carbon(:,islow,l,j)*veget_cov_max(:,l)/(z_soil(j)-z_soil(j-1))
1138   ENDDO
1139 ENDDO
1140 CALL xios_orchidee_send_field("CARBON_SLOW",var_tmp(:,:))
1141 
1142 var_tmp(:,:)=zero
1143 DO l= 1,nvm
1144    var_tmp(:,1)=var_tmp(:,1)+carbon(:,ipassive,l,1)*veget_cov_max(:,l)/z_soil(1)
1145        WHERE (depth_deepsoil(:,l).GT.min_sechiba)
1146        var_tmp(:,nslmd)=var_tmp(:,nslmd)+carbon(:,ipassive,l,nslmd)*veget_cov_max(:,l)/depth_deepsoil(:,l)
1147    ENDWHERE
1148 ENDDO
1149 DO j= 2,nslm
1150   DO l=1,nvm
1151    var_tmp(:,j)=var_tmp(:,j)+carbon(:,ipassive,l,j)*veget_cov_max(:,l)/(z_soil(j)-z_soil(j-1))
1152   ENDDO
1153 ENDDO
1154 CALL xios_orchidee_send_field("CARBON_PASSIVE",var_tmp(:,:))
1155 
1156 var_tmp(:,:)=zero
1157 DO l=1,nvm
1158   DO k= 1,npool
1159      var_tmp(:,1)=var_tmp(:,1)+DOC(:,l,1,ifree,k,icarbon)*veget_cov_max(:,l)/(z_soil(1))
1160      WHERE (depth_deepsoil(:,l).GT.min_sechiba)
1161          var_tmp(:,nslmd)=var_tmp(:,nslmd)+DOC(:,l,nslmd,ifree,k,icarbon)*veget_cov_max(:,l)/depth_deepsoil(:,l)
1162      ENDWHERE
1163   ENDDO
1164 ENDDO
1165
1166 DO j= 2,nslm
1167   DO l=1,nvm
1168     DO k= 1,npool
1169        var_tmp(:,j)=var_tmp(:,j)+DOC(:,l,j,ifree,k,icarbon)*veget_cov_max(:,l)/ &
1170                              (z_soil(j)-z_soil(j-1))
1171     ENDDO
1172   ENDDO
1173 ENDDO
1174 CALL xios_orchidee_send_field("FREE_DOC",var_tmp(:,:))
1175 
1176 var_tmp(:,:)=zero
1177 DO l=1,nvm
1178    DO k=1,npool
1179       var_tmp(:,1)=var_tmp(:,1)+DOC(:,l,1,iadsorbed,k,icarbon)*veget_cov_max(:,l)/(z_soil(1))
1180       WHERE (depth_deepsoil(:,l).GT.min_sechiba)
1181           var_tmp(:,nslmd)=var_tmp(:,nslmd)+DOC(:,l,nslmd,iadsorbed,k,icarbon)*veget_cov_max(:,l)/depth_deepsoil(:,l)
1182       ENDWHERE
1183    ENDDO
1184 ENDDO
1185 DO j= 2,nslm
1186   DO l=1,nvm
1187     DO k=1,npool
1188        var_tmp(:,j)=var_tmp(:,j)+DOC(:,l,j,iadsorbed,k,icarbon)*veget_cov_max(:,l)/ &
1189                              (z_soil(j)-z_soil(j-1))
1190     ENDDO
1191   ENDDO
1192 ENDDO
1193 CALL xios_orchidee_send_field("ADSORBED_DOC",var_tmp(:,:))
1194        !
1195    CALL xios_orchidee_send_field("LITTERHUM",litterhum_daily)
1196    CALL xios_orchidee_send_field("TURNOVER_TIME",turnover_time)
1197    CALL xios_orchidee_send_field("PROD10",prod10)
1198    CALL xios_orchidee_send_field("FLUX10",flux10)
1199    CALL xios_orchidee_send_field("PROD100",prod100)
1200    CALL xios_orchidee_send_field("FLUX100",flux100)
1201    CALL xios_orchidee_send_field("CONVFLUX",convflux)
1202    CALL xios_orchidee_send_field("CFLUX_PROD10",cflux_prod10)
1203    CALL xios_orchidee_send_field("CFLUX_PROD100",cflux_prod100)
1204    CALL xios_orchidee_send_field("PROD10_HARVEST",prod10_harvest)
1205    CALL xios_orchidee_send_field("FLUX10_HARVEST",flux10_harvest)
1206    CALL xios_orchidee_send_field("PROD100_HARVEST",prod100_harvest)
1207    CALL xios_orchidee_send_field("FLUX100_HARVEST",flux100_harvest)
1208    CALL xios_orchidee_send_field("CONVFLUX_HARVEST",convflux_harvest)
1209    CALL xios_orchidee_send_field("CFLUX_PROD10_HARVEST",cflux_prod10_harvest)
1210    CALL xios_orchidee_send_field("CFLUX_PROD100_HARVEST",cflux_prod100_harvest)
1211    CALL xios_orchidee_send_field("WOOD_HARVEST",woodharvest/one_year*dt_days)
1212    CALL xios_orchidee_send_field("WOOD_HARVEST_PFT",woodharvestpft)
1213    CALL xios_orchidee_send_field("HARVEST_ABOVE",harvest_above)
1214    CALL xios_orchidee_send_field("VCMAX",vcmax)
1215    CALL xios_orchidee_send_field("AGE",age)
1216    CALL xios_orchidee_send_field("HEIGHT",height)
1217    CALL xios_orchidee_send_field("FIREINDEX",fireindex(:,:))
1218
1219    ! ipcc history
1220    ! Carbon stock transformed from gC/m2 into kgC/m2
1221    CALL xios_orchidee_send_field("cVeg",SUM(tot_live_biomass(:,:,icarbon)*veget_cov_max,dim=2)/1e3)
1222    CALL xios_orchidee_send_field("cVegGrass",sum_cVegGrass/1e3)
1223    CALL xios_orchidee_send_field("cVegCrop",sum_cVegCrop/1e3)
1224    CALL xios_orchidee_send_field("cVegTree",sum_cVegTree/1e3)
1225    CALL xios_orchidee_send_field("cOther",SUM(cOther(:,:,icarbon)*veget_cov_max,dim=2)/1e3)
1226    ! CALL xios_orchidee_send_field("cLitter",SUM(tot_litter_carb*veget_cov_max,dim=2)/1e3)
1227    CALL xios_orchidee_send_field("cLitterGrass",sum_cLitterGrass/1e3)
1228    CALL xios_orchidee_send_field("cLitterCrop",sum_cLitterCrop/1e3)
1229    CALL xios_orchidee_send_field("cLitterTree",sum_cLitterTree/1e3)
1230    ! CALL xios_orchidee_send_field("cSoil",SUM(tot_soil_carb*veget_cov_max,dim=2)/1e3)
1231    CALL xios_orchidee_send_field("cSoilGrass",sum_cSoilGrass/1e3)
1232    CALL xios_orchidee_send_field("cSoilCrop",sum_cSoilCrop/1e3)
1233    CALL xios_orchidee_send_field("cSoilTree",sum_cSoilTree/1e3)
1234    CALL xios_orchidee_send_field("cProduct",(prod10_total + prod100_total + prod10_harvest_total + prod100_harvest_total)/1e3)
1235
1236    cproductlut(:,:)=xios_default_val
1237    WHERE(fraclut(:,id_psl) > min_sechiba)
1238       cproductlut(:,id_psl)= ( prod10_total + prod100_total + prod10_harvest_total + prod100_harvest_total)/1e3 &
1239            / fraclut(:,id_psl)
1240    ENDWHERE
1241    CALL xios_orchidee_send_field("cproductlut",cproductlut)
1242
1243    CALL xios_orchidee_send_field("cMassVariation",carb_mass_variation/1e3/one_day)
1244
1245    CALL xios_orchidee_send_field("lai_ipcc",SUM(lai*veget_cov_max,dim=2)) ! m2/m2
1246   
1247    ! Carbon fluxes transformed from gC/m2/d into kgC/m2/s
1248    CALL xios_orchidee_send_field("gpp_ipcc",SUM(gpp_daily*veget_cov_max,dim=2)/1e3/one_day)
1249    CALL xios_orchidee_send_field("ra",SUM((resp_maint+resp_growth)*veget_cov_max,dim=2)/1e3/one_day)
1250    vartmp(:)=zero
1251    DO j = 2, nvm
1252       IF ( .NOT. is_tree(j) .AND. natural(j) ) THEN
1253          vartmp(:) = vartmp(:) + (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)
1254       ENDIF
1255    ENDDO
1256    CALL xios_orchidee_send_field("raGrass",vartmp/1e3/one_day)
1257
1258    vartmp(:)=zero
1259    DO j = 2, nvm
1260       IF (( .NOT. is_tree(j)) .AND. (.NOT. natural(j)) ) THEN
1261          vartmp(:) = vartmp(:) + (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)
1262       ENDIF
1263    ENDDO
1264    CALL xios_orchidee_send_field("raCrop",vartmp/1e3/one_day)
1265
1266    vartmp(:)=zero
1267    DO j = 2, nvm
1268       IF ( is_tree(j) ) THEN
1269          vartmp(:) = vartmp(:) + (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)
1270       ENDIF
1271    ENDDO
1272    CALL xios_orchidee_send_field("raTree",vartmp/1e3/one_day)
1273
1274    CALL xios_orchidee_send_field("npp_ipcc",SUM(npp_daily*veget_cov_max,dim=2)/1e3/one_day)
1275    vartmp(:)=zero
1276    DO j = 2, nvm
1277       IF ( .NOT. is_tree(j) .AND. natural(j) ) THEN
1278          vartmp(:) = vartmp(:) + npp_daily(:,j)*veget_cov_max(:,j)
1279       ENDIF
1280    ENDDO
1281    CALL xios_orchidee_send_field("nppGrass",vartmp/1e3/one_day)
1282
1283    DO j = 2, nvm
1284       IF ( (.NOT. is_tree(j)) .AND. (.NOT. natural(j)) ) THEN
1285          vartmp(:) = vartmp(:) + npp_daily(:,j)*veget_cov_max(:,j)
1286       ENDIF
1287    ENDDO
1288    CALL xios_orchidee_send_field("nppCrop",vartmp/1e3/one_day)
1289
1290    vartmp(:)=zero
1291    DO j = 2, nvm
1292       IF ( is_tree(j) ) THEN
1293          vartmp(:) = vartmp(:) + npp_daily(:,j)*veget_cov_max(:,j)
1294       ENDIF
1295    ENDDO
1296    CALL xios_orchidee_send_field("nppTree",vartmp/1e3/one_day)
1297
1298    CALL xios_orchidee_send_field("rh",SUM(resp_hetero*veget_cov_max,dim=2)/1e3/one_day)
1299    CALL xios_orchidee_send_field("fFire",SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day)
1300    ctotfirelut(:,:)=xios_default_val
1301    WHERE(fraclut(:,id_psl) > min_sechiba) 
1302       ctotfirelut(:,id_psl)= SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day &
1303            / fraclut(:,id_psl)
1304    ENDWHERE
1305    CALL xios_orchidee_send_field("ctotfirelut",ctotfirelut)
1306
1307    CALL xios_orchidee_send_field("fHarvest",harvest_above/1e3/one_day)
1308    Call xios_orchidee_send_field("fLuc",cflux_prod_total/1e3/one_day)
1309    CALL xios_orchidee_send_field("fWoodharvest",cflux_prod_harvest_total/1e3/one_day)
1310    CALL xios_orchidee_send_field("fDeforestToProduct",SUM((fDeforestToProduct-convfluxpft),dim=2)/1e3/one_day)
1311   
1312    flulccproductlut(:,:)=xios_default_val
1313    flulccproductlut(:,id_psl)=0.
1314    flulccproductlut(:,id_crp)=0.
1315    DO j=1,nvm
1316       IF(natural(j)) THEN
1317          flulccproductlut(:,id_psl)=flulccproductlut(:,id_psl) + fDeforestToProduct(:,j)
1318       ELSE
1319          flulccproductlut(:,id_crp)=flulccproductlut(:,id_crp) + fDeforestToProduct(:,j)
1320       ENDIF
1321    ENDDO
1322
1323    WHERE(fraclut(:,id_psl) > min_sechiba)
1324       flulccproductlut(:,id_psl)= (flulccproductlut(:,id_psl) &
1325            +SUM(fHarvestToProduct,dim=2)) &
1326            /1e3/one_day / fraclut(:,id_psl)
1327    ELSEWHERE
1328       flulccproductlut(:,id_psl)= xios_default_val
1329    ENDWHERE
1330    WHERE(fraclut(:,id_crp) > min_sechiba)
1331       flulccproductlut(:,id_crp)= flulccproductlut(:,id_crp) &
1332            /1e3/one_day / fraclut(:,id_crp)
1333    ELSEWHERE
1334       flulccproductlut(:,id_crp)= xios_default_val
1335    ENDWHERE
1336    CALL xios_orchidee_send_field("flulccproductlut",flulccproductlut)
1337
1338    flulccresiduelut(:,:)=xios_default_val
1339    flulccresiduelut(:,id_psl)=0.
1340    flulccresiduelut(:,id_crp)=0.
1341    DO j=1,nvm
1342       IF(natural(j)) THEN
1343          flulccresiduelut(:,id_psl) = flulccresiduelut(:,id_psl) + fLulccResidue(:,j)
1344       ELSE
1345          flulccresiduelut(:,id_crp) = flulccresiduelut(:,id_crp) + fLulccResidue(:,j)
1346       ENDIF
1347    ENDDO
1348
1349    WHERE(fraclut(:,id_psl) > min_sechiba)
1350       flulccresiduelut(:,id_psl)= flulccresiduelut(:,id_psl)/1e3/one_day/fraclut(:,id_psl)
1351    ELSEWHERE
1352       flulccresiduelut(:,id_psl)= xios_default_val
1353    ENDWHERE
1354    WHERE(fraclut(:,id_crp) > min_sechiba)
1355       flulccresiduelut(:,id_crp)= flulccresiduelut(:,id_crp)/1e3/one_day /fraclut(:,id_crp)
1356    ELSEWHERE
1357       flulccresiduelut(:,id_crp)= xios_default_val
1358    ENDWHERE
1359    CALL xios_orchidee_send_field("flulccresiduelut",flulccresiduelut)
1360    CALL xios_orchidee_send_field("flulccresidue",flulccresidue)
1361
1362    flulccatmlut(:,:)=xios_default_val
1363    flulccatmlut(:,id_psl)=0.
1364    flulccatmlut(:,id_crp)=0.
1365    DO j=1,nvm
1366       IF(natural(j)) THEN
1367          flulccatmlut(:,id_psl)=flulccatmlut(:,id_psl)+convfluxpft(:,j)
1368       ELSE
1369          flulccatmlut(:,id_crp)=flulccatmlut(:,id_crp)+convfluxpft(:,j)
1370       ENDIF
1371    ENDDO
1372
1373    WHERE(fraclut(:,id_psl) > min_sechiba)
1374       flulccatmlut(:,id_psl)= (flulccatmlut(:,id_psl)+cflux_prod10+cflux_prod100+cflux_prod_harvest_total) &
1375            /1e3/one_day / fraclut(:,id_psl)
1376    ELSEWHERE
1377       flulccatmlut(:,id_psl)= xios_default_val
1378    ENDWHERE
1379    WHERE(fraclut(:,id_crp) > min_sechiba)
1380       flulccatmlut(:,id_crp)= (flulccatmlut(:,id_crp)+harvest_above)/1e3/one_day / fraclut(:,id_crp)
1381    ELSEWHERE
1382       flulccatmlut(:,id_crp)= xios_default_val
1383    ENDWHERE
1384    CALL xios_orchidee_send_field("flulccatmlut",flulccatmlut)
1385
1386   ! co2_to_bm is not added as it is already included in gpp
1387    CALL xios_orchidee_send_field("nbp",(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) * &
1388          veget_cov_max,dim=2)-cflux_prod_total-cflux_prod_harvest_total-harvest_above)/1e3/one_day)
1389    CALL xios_orchidee_send_field("fVegLitter",SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*&
1390         veget_cov_max,dim=2)/1e3/one_day)
1391!    CALL xios_orchidee_send_field("fLitterSoil",SUM(SUM(soilcarbon_input,dim=2)*veget_cov_max,dim=2)/1e3/one_day)
1392
1393    ! Carbon stock transformed from gC/m2 into kgC/m2
1394    CALL xios_orchidee_send_field("cLeaf",SUM(biomass(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3)
1395    CALL xios_orchidee_send_field("cStem",SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*&
1396          veget_cov_max,dim=2)/1e3)
1397    CALL xios_orchidee_send_field("cWood",SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon) + &
1398         biomass(:,:,isapbelow,icarbon) + biomass(:,:,iheartbelow,icarbon))* &
1399         veget_cov_max,dim=2)/1e3)
1400    CALL xios_orchidee_send_field("cRoot",SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + &
1401         biomass(:,:,iheartbelow,icarbon) )*veget_cov_max,dim=2)/1e3)
1402    CALL xios_orchidee_send_field("cMisc",SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*&
1403         veget_cov_max,dim=2)/1e3)
1404    ! CALL xios_orchidee_send_field("cLitterAbove",SUM((litter(:,istructural,:,iabove,icarbon)+&
1405         ! litter(:,imetabolic,:,iabove,icarbon))*veget_cov_max,dim=2)/1e3)
1406    ! CALL xios_orchidee_send_field("cLitterBelow",SUM((litter(:,istructural,:,ibelow,icarbon)+&
1407         ! litter(:,imetabolic,:,ibelow,icarbon))*veget_cov_max,dim=2)/1e3)
1408    ! CALL xios_orchidee_send_field("cSoilFast",SUM(carbon(:,iactive,:)*veget_cov_max,dim=2)/1e3)
1409    ! CALL xios_orchidee_send_field("cSoilMedium",SUM(carbon(:,islow,:)*veget_cov_max,dim=2)/1e3)
1410    ! CALL xios_orchidee_send_field("cSoilSlow",SUM(carbon(:,ipassive,:)*veget_cov_max,dim=2)/1e3)
1411
1412    ! DO k = 1, ncarb
1413       ! csoilpools(:,k) = SUM(carbon(:,k,:)*veget_cov_max,dim=2)/1e3
1414    ! END DO
1415    ! CALL xios_orchidee_send_field("cSoilPools",csoilpools)
1416
1417    ! Vegetation fractions [0,100]
1418    CALL xios_orchidee_send_field("landCoverFrac",veget_cov_max*100)
1419
1420    ! Carbon fluxes transformed from gC/m2/d into kgC/m2/s
1421    CALL xios_orchidee_send_field("rGrowth",SUM(resp_growth*veget_cov_max,dim=2)/1e3/one_day)
1422    CALL xios_orchidee_send_field("rMaint",SUM(resp_maint*veget_cov_max,dim=2)/1e3/one_day)
1423    CALL xios_orchidee_send_field("nppLeaf",SUM(bm_alloc(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3/one_day)
1424    ! nppStem : from wood above surface
1425    CALL xios_orchidee_send_field("nppStem",SUM((bm_alloc(:,:,isapabove,icarbon) + bm_alloc(:,:,iheartabove,icarbon)) &
1426         *veget_cov_max,dim=2)/1e3/one_day)
1427    ! nppWood : from wood above and below surface
1428    CALL xios_orchidee_send_field("nppWood",SUM((bm_alloc(:,:,isapabove,icarbon) + bm_alloc(:,:,iheartabove,icarbon) + &
1429         bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iheartbelow,icarbon)) &
1430         *veget_cov_max,dim=2)/1e3/one_day)
1431    ! nppRoot : from wood below surface and fine roots
1432    CALL xios_orchidee_send_field("nppRoot",SUM((bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iheartbelow,icarbon) + &
1433         bm_alloc(:,:,iroot,icarbon)) * veget_cov_max,dim=2)/1e3/one_day)
1434    CALL xios_orchidee_send_field("nppOther",SUM(( bm_alloc(:,:,ifruit,icarbon) + bm_alloc(:,:,icarbres,icarbon) ) * &
1435         veget_cov_max,dim=2)/1e3/one_day)
1436
1437    ! Calculate variables according to LUMIP specifications.
1438    clitterlut(:,:) = 0 
1439    csoillut(:,:) = 0 
1440    cveglut(:,:) = 0
1441    lailut(:,:) = 0
1442    ralut(:,:) = 0
1443    rhlut(:,:) = 0
1444    npplut(:,:) = 0
1445
1446    DO j=1,nvm
1447       IF (natural(j)) THEN
1448          ! clitterlut(:,id_psl) = clitterlut(:,id_psl) + tot_litter_carb(:,j)*veget_cov_max(:,j)/1e3
1449          ! csoillut(:,id_psl) = csoillut(:,id_psl) + tot_soil_carb(:,j)*veget_cov_max(:,j)/1e3
1450          cveglut(:,id_psl) = cveglut(:,id_psl) + tot_live_biomass(:,j,icarbon)*veget_cov_max(:,j)/1e3
1451          lailut(:,id_psl) = lailut(:,id_psl) + lai(:,j)*veget_cov_max(:,j)
1452          ralut(:,id_psl) = ralut(:,id_psl) + &
1453               (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)/1e3/one_day
1454          rhlut(:,id_psl) = rhlut(:,id_psl) + &
1455               resp_hetero(:,j)*veget_cov_max(:,j)/1e3/one_day
1456          npplut(:,id_psl) = npplut(:,id_psl) + &
1457               npp_daily(:,j)*veget_cov_max(:,j)/1e3/one_day
1458       ELSE
1459          ! clitterlut(:,id_crp) = clitterlut(:,id_crp) + tot_litter_carb(:,j)*veget_cov_max(:,j)/1e3
1460          ! csoillut(:,id_crp) = csoillut(:,id_crp) + tot_soil_carb(:,j)*veget_cov_max(:,j)/1e3
1461          cveglut(:,id_crp) = cveglut(:,id_crp) + tot_live_biomass(:,j,icarbon)*veget_cov_max(:,j)/1e3
1462          lailut(:,id_crp) = lailut(:,id_crp) + lai(:,j)*veget_cov_max(:,j)
1463          ralut(:,id_crp) = ralut(:,id_crp) + &
1464               (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)/1e3/one_day
1465          rhlut(:,id_crp) = rhlut(:,id_crp) + &
1466               resp_hetero(:,j)*veget_cov_max(:,j)/1e3/one_day
1467          npplut(:,id_crp) = npplut(:,id_crp) + &
1468               npp_daily(:,j)*veget_cov_max(:,j)/1e3/one_day
1469       END IF
1470    END DO
1471
1472    WHERE (fraclut(:,id_psl)>min_sechiba)
1473       clitterlut(:,id_psl) = clitterlut(:,id_psl)/fraclut(:,id_psl)
1474       csoillut(:,id_psl) = csoillut(:,id_psl)/fraclut(:,id_psl)
1475       cveglut(:,id_psl) = cveglut(:,id_psl)/fraclut(:,id_psl)
1476       lailut(:,id_psl) = lailut(:,id_psl)/fraclut(:,id_psl)
1477       ralut(:,id_psl) = ralut(:,id_psl)/fraclut(:,id_psl)
1478       rhlut(:,id_psl) = rhlut(:,id_psl)/fraclut(:,id_psl)
1479       npplut(:,id_psl) = npplut(:,id_psl)/fraclut(:,id_psl)
1480    ELSEWHERE
1481       clitterlut(:,id_psl) = xios_default_val
1482       csoillut(:,id_psl) = xios_default_val
1483       cveglut(:,id_psl) = xios_default_val
1484       lailut(:,id_psl) = xios_default_val
1485       ralut(:,id_psl) = xios_default_val
1486       rhlut(:,id_psl) = xios_default_val
1487       npplut(:,id_psl) = xios_default_val
1488    END WHERE
1489
1490    WHERE (fraclut(:,id_crp)>min_sechiba)
1491       clitterlut(:,id_crp) = clitterlut(:,id_crp)/fraclut(:,id_crp)
1492       csoillut(:,id_crp) = csoillut(:,id_crp)/fraclut(:,id_crp)
1493       cveglut(:,id_crp) = cveglut(:,id_crp)/fraclut(:,id_crp)
1494       lailut(:,id_crp) = lailut(:,id_crp)/fraclut(:,id_crp)
1495       ralut(:,id_crp) = ralut(:,id_crp)/fraclut(:,id_crp)
1496       rhlut(:,id_crp) = rhlut(:,id_crp)/fraclut(:,id_crp)
1497       npplut(:,id_crp) = npplut(:,id_crp)/fraclut(:,id_crp)
1498    ELSEWHERE
1499       clitterlut(:,id_crp) = xios_default_val
1500       csoillut(:,id_crp) = xios_default_val
1501       cveglut(:,id_crp) = xios_default_val
1502       lailut(:,id_crp) = xios_default_val
1503       ralut(:,id_crp) = xios_default_val
1504       rhlut(:,id_crp) = xios_default_val
1505       npplut(:,id_crp) = xios_default_val
1506    END WHERE
1507
1508    clitterlut(:,id_pst) = xios_default_val
1509    clitterlut(:,id_urb) = xios_default_val
1510    csoillut(:,id_pst)   = xios_default_val
1511    csoillut(:,id_urb)   = xios_default_val
1512    cveglut(:,id_pst)    = xios_default_val
1513    cveglut(:,id_urb)    = xios_default_val
1514    lailut(:,id_pst)     = xios_default_val
1515    lailut(:,id_urb)     = xios_default_val
1516    ralut(:,id_pst)      = xios_default_val
1517    ralut(:,id_urb)      = xios_default_val
1518    rhlut(:,id_pst)      = xios_default_val
1519    rhlut(:,id_urb)      = xios_default_val
1520    npplut(:,id_pst)     = xios_default_val
1521    npplut(:,id_urb)     = xios_default_val
1522
1523    CALL xios_orchidee_send_field("clitterlut",clitterlut)
1524    CALL xios_orchidee_send_field("csoillut",csoillut)
1525    CALL xios_orchidee_send_field("cveglut",cveglut)
1526    CALL xios_orchidee_send_field("lailut",lailut)
1527    CALL xios_orchidee_send_field("ralut",ralut)
1528    CALL xios_orchidee_send_field("rhlut",rhlut)
1529    CALL xios_orchidee_send_field("npplut",npplut)
1530
1531    CALL histwrite_p (hist_id_stomate, 'RESOLUTION_X', itime, &
1532         resolution(:,1), npts, hori_index)
1533    CALL histwrite_p (hist_id_stomate, 'RESOLUTION_Y', itime, &
1534         resolution(:,2), npts, hori_index)
1535    CALL histwrite_p (hist_id_stomate, 'CONTFRAC', itime, &
1536         contfrac(:), npts, hori_index)
1537
1538    ! CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AB', itime, &
1539         ! litter(:,istructural,:,iabove,icarbon), npts*nvm, horipft_index)
1540    ! CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AB', itime, &
1541         ! litter(:,imetabolic,:,iabove,icarbon), npts*nvm, horipft_index)
1542    ! CALL histwrite_p (hist_id_stomate, 'LITTER_STR_BE', itime, &
1543         ! litter(:,istructural,:,ibelow,icarbon), npts*nvm, horipft_index)
1544    ! CALL histwrite_p (hist_id_stomate, 'LITTER_MET_BE', itime, &
1545         ! litter(:,imetabolic,:,ibelow,icarbon), npts*nvm, horipft_index)
1546
1547    CALL histwrite_p (hist_id_stomate, 'DEADLEAF_COVER', itime, &
1548         deadleaf_cover, npts, hori_index)
1549
1550    ! CALL histwrite_p (hist_id_stomate, 'TOTAL_SOIL_CARB', itime, &
1551         ! tot_litter_soil_carb, npts*nvm, horipft_index)
1552    ! CALL histwrite_p (hist_id_stomate, 'CARBON_ACTIVE', itime, &
1553         ! carbon(:,iactive,:), npts*nvm, horipft_index)
1554    ! CALL histwrite_p (hist_id_stomate, 'CARBON_SLOW', itime, &
1555         ! carbon(:,islow,:), npts*nvm, horipft_index)
1556    ! CALL histwrite_p (hist_id_stomate, 'CARBON_PASSIVE', itime, &
1557         ! carbon(:,ipassive,:), npts*nvm, horipft_index)
1558
1559    CALL histwrite_p (hist_id_stomate, 'T2M_MONTH', itime, &
1560         t2m_month, npts, hori_index)
1561    CALL histwrite_p (hist_id_stomate, 'T2M_WEEK', itime, &
1562         t2m_week, npts, hori_index)
1563    CALL histwrite_p (hist_id_stomate, 'TSEASON', itime, &
1564         Tseason, npts, hori_index)
1565    CALL histwrite_p (hist_id_stomate, 'TMIN_SPRING_TIME', itime, &
1566         Tmin_spring_time, npts*nvm, horipft_index)
1567    CALL histwrite_p (hist_id_stomate, 'ONSET_DATE', itime, &
1568         onset_date(:,:), npts*nvm, horipft_index)
1569    CALL histwrite_p (hist_id_stomate, 'FPC_MAX', itime, &
1570         fpc_max, npts*nvm, horipft_index)
1571    CALL histwrite_p (hist_id_stomate, 'MAXFPC_LASTYEAR', itime, &
1572         maxfpc_lastyear, npts*nvm, horipft_index)
1573    CALL histwrite_p (hist_id_stomate, 'HET_RESP', itime, &
1574         resp_hetero(:,:), npts*nvm, horipft_index)
1575    CALL histwrite_p (hist_id_stomate, 'FIREINDEX', itime, &
1576         fireindex(:,:), npts*nvm, horipft_index)
1577    CALL histwrite_p (hist_id_stomate, 'LITTERHUM', itime, &
1578         litterhum_daily, npts, hori_index)
1579    CALL histwrite_p (hist_id_stomate, 'CO2_FIRE', itime, &
1580         co2_fire, npts*nvm, horipft_index)
1581    CALL histwrite_p (hist_id_stomate, 'CO2_TAKEN', itime, &
1582         co2_to_bm, npts*nvm, horipft_index)
1583    ! land cover change
1584    CALL histwrite_p (hist_id_stomate, 'CONVFLUX', itime, &
1585         convflux, npts, hori_index)
1586    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10', itime, &
1587         cflux_prod10, npts, hori_index)
1588    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100', itime, &
1589         cflux_prod100, npts, hori_index)
1590    CALL histwrite_p (hist_id_stomate, 'CONVFLUX_HARVEST', itime, &
1591         convflux_harvest, npts, hori_index)
1592    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10_HARVEST', itime, &
1593         cflux_prod10_harvest, npts, hori_index)
1594    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100_HARVES', itime, &
1595         cflux_prod100_harvest, npts, hori_index)
1596    CALL histwrite_p (hist_id_stomate, 'WOOD_HARVEST', itime, &
1597         woodharvest/one_year*dt_days, npts, hori_index)
1598    CALL histwrite_p (hist_id_stomate, 'WOOD_HARVEST_PFT', itime, &
1599         woodharvestpft, npts*nvm, horipft_index)
1600
1601    CALL histwrite_p (hist_id_stomate, 'HARVEST_ABOVE', itime, &
1602         harvest_above, npts, hori_index)
1603
1604    CALL histwrite_p (hist_id_stomate, 'LAI', itime, &
1605         lai, npts*nvm, horipft_index)
1606    CALL histwrite_p (hist_id_stomate, 'VEGET_COV_MAX', itime, &
1607         veget_cov_max, npts*nvm, horipft_index)
1608    CALL histwrite_p (hist_id_stomate, 'NPP', itime, &
1609         npp_daily, npts*nvm, horipft_index)
1610    CALL histwrite_p (hist_id_stomate, 'GPP', itime, &
1611         gpp_daily, npts*nvm, horipft_index)
1612    CALL histwrite_p (hist_id_stomate, 'IND', itime, &
1613         ind, npts*nvm, horipft_index)
1614    CALL histwrite_p (hist_id_stomate, 'CN_IND', itime, &
1615         cn_ind, npts*nvm, horipft_index)
1616    CALL histwrite_p (hist_id_stomate, 'WOODMASS_IND', itime, &
1617         woodmass_ind, npts*nvm, horipft_index)
1618    CALL histwrite_p (hist_id_stomate, 'TOTAL_M', itime, &
1619         tot_live_biomass(:,:,icarbon), npts*nvm, horipft_index)
1620    CALL histwrite_p (hist_id_stomate, 'LEAF_M', itime, &
1621         biomass(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1622    CALL histwrite_p (hist_id_stomate, 'SAP_M_AB', itime, &
1623         biomass(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1624    CALL histwrite_p (hist_id_stomate, 'SAP_M_BE', itime, &
1625         biomass(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
1626    CALL histwrite_p (hist_id_stomate, 'HEART_M_AB', itime, &
1627         biomass(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
1628    CALL histwrite_p (hist_id_stomate, 'HEART_M_BE', itime, &
1629         biomass(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
1630    CALL histwrite_p (hist_id_stomate, 'ROOT_M', itime, &
1631         biomass(:,:,iroot,icarbon), npts*nvm, horipft_index)
1632    CALL histwrite_p (hist_id_stomate, 'FRUIT_M', itime, &
1633         biomass(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1634    CALL histwrite_p (hist_id_stomate, 'RESERVE_M', itime, &
1635         biomass(:,:,icarbres,icarbon), npts*nvm, horipft_index)
1636    CALL histwrite_p (hist_id_stomate, 'TOTAL_TURN', itime, &
1637         tot_turnover(:,:,icarbon), npts*nvm, horipft_index)
1638    CALL histwrite_p (hist_id_stomate, 'LEAF_TURN', itime, &
1639         turnover_daily(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1640    CALL histwrite_p (hist_id_stomate, 'SAP_AB_TURN', itime, &
1641         turnover_daily(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1642    CALL histwrite_p (hist_id_stomate, 'ROOT_TURN', itime, &
1643         turnover_daily(:,:,iroot,icarbon), npts*nvm, horipft_index)
1644    CALL histwrite_p (hist_id_stomate, 'FRUIT_TURN', itime, &
1645         turnover_daily(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1646    CALL histwrite_p (hist_id_stomate, 'TOTAL_BM_LITTER', itime, &
1647         tot_bm_to_litter(:,:,icarbon), npts*nvm, horipft_index)
1648    CALL histwrite_p (hist_id_stomate, 'LEAF_BM_LITTER', itime, &
1649         bm_to_litter(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1650    CALL histwrite_p (hist_id_stomate, 'SAP_AB_BM_LITTER', itime, &
1651         bm_to_litter(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1652    CALL histwrite_p (hist_id_stomate, 'SAP_BE_BM_LITTER', itime, &
1653         bm_to_litter(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
1654    CALL histwrite_p (hist_id_stomate, 'HEART_AB_BM_LITTER', itime, &
1655         bm_to_litter(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
1656    CALL histwrite_p (hist_id_stomate, 'HEART_BE_BM_LITTER', itime, &
1657         bm_to_litter(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
1658    CALL histwrite_p (hist_id_stomate, 'ROOT_BM_LITTER', itime, &
1659         bm_to_litter(:,:,iroot,icarbon), npts*nvm, horipft_index)
1660    CALL histwrite_p (hist_id_stomate, 'FRUIT_BM_LITTER', itime, &
1661         bm_to_litter(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1662    CALL histwrite_p (hist_id_stomate, 'RESERVE_BM_LITTER', itime, &
1663         bm_to_litter(:,:,icarbres,icarbon), npts*nvm, horipft_index)
1664    CALL histwrite_p (hist_id_stomate, 'MAINT_RESP', itime, &
1665         resp_maint, npts*nvm, horipft_index)
1666    CALL histwrite_p (hist_id_stomate, 'GROWTH_RESP', itime, &
1667         resp_growth, npts*nvm, horipft_index)
1668    CALL histwrite_p (hist_id_stomate, 'AGE', itime, &
1669         age, npts*nvm, horipft_index)
1670    CALL histwrite_p (hist_id_stomate, 'HEIGHT', itime, &
1671         height, npts*nvm, horipft_index)
1672    CALL histwrite_p (hist_id_stomate, 'MOISTRESS', itime, &
1673         moiavail_week, npts*nvm, horipft_index)
1674    CALL histwrite_p (hist_id_stomate, 'VCMAX', itime, &
1675         vcmax, npts*nvm, horipft_index)
1676    CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME', itime, &
1677         turnover_time, npts*nvm, horipft_index)
1678    ! land cover change
1679    CALL histwrite_p (hist_id_stomate, 'PROD10', itime, &
1680         prod10, npts*11, horip11_index)
1681    CALL histwrite_p (hist_id_stomate, 'PROD100', itime, &
1682         prod100, npts*101, horip101_index)
1683    CALL histwrite_p (hist_id_stomate, 'FLUX10', itime, &
1684         flux10, npts*10, horip10_index)
1685    CALL histwrite_p (hist_id_stomate, 'FLUX100', itime, &
1686         flux100, npts*100, horip100_index)
1687    CALL histwrite_p (hist_id_stomate, 'PROD10_HARVEST', itime, &
1688         prod10_harvest, npts*11, horip11_index)
1689    CALL histwrite_p (hist_id_stomate, 'PROD100_HARVEST', itime, &
1690         prod100_harvest, npts*101, horip101_index)
1691    CALL histwrite_p (hist_id_stomate, 'FLUX10_HARVEST', itime, &
1692         flux10_harvest, npts*10, horip10_index)
1693    CALL histwrite_p (hist_id_stomate, 'FLUX100_HARVEST', itime, &
1694         flux100_harvest, npts*100, horip100_index)
1695
1696    IF ( hist_id_stomate_IPCC > 0 ) THEN
1697       vartmp(:)=SUM(tot_live_biomass(:,:,icarbon)*veget_cov_max,dim=2)/1e3
1698       CALL histwrite_p (hist_id_stomate_IPCC, "cVeg", itime, &
1699            vartmp, npts, hori_index)
1700       ! vartmp(:)=SUM(tot_litter_carb*veget_cov_max,dim=2)/1e3
1701       ! CALL histwrite_p (hist_id_stomate_IPCC, "cLitter", itime, &
1702            ! vartmp, npts, hori_index)
1703       ! vartmp(:)=SUM(tot_soil_carb*veget_cov_max,dim=2)/1e3
1704       ! CALL histwrite_p (hist_id_stomate_IPCC, "cSoil", itime, &
1705            ! vartmp, npts, hori_index)
1706       vartmp(:)=(prod10_total + prod100_total + prod10_harvest_total + prod100_harvest_total)/1e3
1707       CALL histwrite_p (hist_id_stomate_IPCC, "cProduct", itime, &
1708            vartmp, npts, hori_index)
1709       vartmp(:)=carb_mass_variation/1e3/one_day
1710       CALL histwrite_p (hist_id_stomate_IPCC, "cMassVariation", itime, &
1711            vartmp, npts, hori_index)
1712       vartmp(:)=SUM(lai*veget_cov_max,dim=2)
1713       CALL histwrite_p (hist_id_stomate_IPCC, "lai", itime, &
1714            vartmp, npts, hori_index)
1715       vartmp(:)=SUM(gpp_daily*veget_cov_max,dim=2)/1e3/one_day
1716       CALL histwrite_p (hist_id_stomate_IPCC, "gpp", itime, &
1717            vartmp, npts, hori_index)
1718       vartmp(:)=SUM((resp_maint+resp_growth)*veget_cov_max,dim=2)/1e3/one_day
1719       CALL histwrite_p (hist_id_stomate_IPCC, "ra", itime, &
1720            vartmp, npts, hori_index)
1721       vartmp(:)=SUM(npp_daily*veget_cov_max,dim=2)/1e3/one_day
1722       CALL histwrite_p (hist_id_stomate_IPCC, "npp", itime, &
1723            vartmp, npts, hori_index)
1724       vartmp(:)=SUM(resp_hetero*veget_cov_max,dim=2)/1e3/one_day
1725       CALL histwrite_p (hist_id_stomate_IPCC, "rh", itime, &
1726            vartmp, npts, hori_index)
1727       vartmp(:)=SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day
1728       CALL histwrite_p (hist_id_stomate_IPCC, "fFire", itime, &
1729            vartmp, npts, hori_index)
1730       vartmp(:)=harvest_above/1e3/one_day
1731       CALL histwrite_p (hist_id_stomate_IPCC, "fHarvest", itime, &
1732            vartmp, npts, hori_index)
1733       vartmp(:)=cflux_prod_total/1e3/one_day
1734       CALL histwrite_p (hist_id_stomate_IPCC, "fLuc", itime, &
1735            vartmp, npts, hori_index)
1736       vartmp(:)=cflux_prod_harvest_total/1e3/one_day
1737       CALL histwrite_p (hist_id_stomate_IPCC, "fWoodharvest", itime, &
1738            vartmp, npts, hori_index)
1739       ! co2_to_bm is not added as it is already included in gpp
1740       vartmp(:)=(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
1741            &        *veget_cov_max,dim=2)-cflux_prod_total-cflux_prod_harvest_total-harvest_above)/1e3/one_day
1742       CALL histwrite_p (hist_id_stomate_IPCC, "nbp", itime, &
1743            vartmp, npts, hori_index)
1744       vartmp(:)=SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*veget_cov_max,dim=2)/1e3/one_day
1745       CALL histwrite_p (hist_id_stomate_IPCC, "fVegLitter", itime, &
1746            vartmp, npts, hori_index)
1747       ! vartmp(:)=SUM(SUM(soilcarbon_input,dim=2)*veget_cov_max,dim=2)/1e3/one_day
1748       ! CALL histwrite_p (hist_id_stomate_IPCC, "fLitterSoil", itime, &
1749            ! vartmp, npts, hori_index)
1750       vartmp(:)=SUM(biomass(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3
1751       CALL histwrite_p (hist_id_stomate_IPCC, "cLeaf", itime, &
1752            vartmp, npts, hori_index)
1753       vartmp(:)=SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*veget_cov_max,dim=2)/1e3
1754       CALL histwrite_p (hist_id_stomate_IPCC, "cStem", itime, &
1755            vartmp, npts, hori_index)
1756       vartmp(:)=SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + biomass(:,:,iheartbelow,icarbon) ) &
1757            &        *veget_cov_max,dim=2)/1e3
1758       CALL histwrite_p (hist_id_stomate_IPCC, "cRoot", itime, &
1759            vartmp, npts, hori_index)
1760       vartmp(:)=SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*veget_cov_max,dim=2)/1e3
1761       CALL histwrite_p (hist_id_stomate_IPCC, "cMisc", itime, &
1762            vartmp, npts, hori_index)
1763       ! vartmp(:)=SUM((litter(:,istructural,:,iabove,icarbon)+litter(:,imetabolic,:,iabove,icarbon))*&
1764            ! veget_cov_max,dim=2)/1e3
1765       ! CALL histwrite_p (hist_id_stomate_IPCC, "cLitterAbove", itime, &
1766            ! vartmp, npts, hori_index)
1767       ! vartmp(:)=SUM((litter(:,istructural,:,ibelow,icarbon)+litter(:,imetabolic,:,ibelow,icarbon))*&
1768            ! veget_cov_max,dim=2)/1e3
1769       ! CALL histwrite_p (hist_id_stomate_IPCC, "cLitterBelow", itime, &
1770            ! vartmp, npts, hori_index)
1771       ! vartmp(:)=SUM(carbon(:,iactive,:)*veget_cov_max,dim=2)/1e3
1772       ! CALL histwrite_p (hist_id_stomate_IPCC, "cSoilFast", itime, &
1773            ! vartmp, npts, hori_index)
1774       ! vartmp(:)=SUM(carbon(:,islow,:)*veget_cov_max,dim=2)/1e3
1775       ! CALL histwrite_p (hist_id_stomate_IPCC, "cSoilMedium", itime, &
1776            ! vartmp, npts, hori_index)
1777       ! vartmp(:)=SUM(carbon(:,ipassive,:)*veget_cov_max,dim=2)/1e3
1778       ! CALL histwrite_p (hist_id_stomate_IPCC, "cSoilSlow", itime, &
1779            ! vartmp, npts, hori_index)
1780       DO j=1,nvm
1781          histvar(:,j)=veget_cov_max(:,j)*100
1782       ENDDO
1783       CALL histwrite_p (hist_id_stomate_IPCC, "landCoverFrac", itime, &
1784            histvar, npts*nvm, horipft_index)
1785       !-
1786       vartmp(:)=zero
1787       DO j = 2,nvm
1788          IF (is_deciduous(j)) THEN
1789             vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1790          ENDIF
1791       ENDDO
1792       CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimDec", itime, &
1793            vartmp, npts, hori_index)
1794       !-
1795       vartmp(:)=zero
1796       DO j = 2,nvm
1797          IF (is_evergreen(j)) THEN
1798             vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1799          ENDIF
1800       ENDDO
1801       CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimEver", itime, &
1802            vartmp, npts, hori_index)
1803       !-
1804       vartmp(:)=zero
1805       DO j = 2,nvm
1806          IF ( .NOT.(is_c4(j)) ) THEN
1807             vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1808          ENDIF
1809       ENDDO
1810       CALL histwrite_p (hist_id_stomate_IPCC, "c3PftFrac", itime, &
1811            vartmp, npts, hori_index)
1812       !-
1813       vartmp(:)=zero
1814       DO j = 2,nvm
1815          IF ( is_c4(j) ) THEN
1816             vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1817          ENDIF
1818       ENDDO
1819       CALL histwrite_p (hist_id_stomate_IPCC, "c4PftFrac", itime, &
1820            vartmp, npts, hori_index)
1821       !-
1822       vartmp(:)=SUM(resp_growth*veget_cov_max,dim=2)/1e3/one_day
1823       CALL histwrite_p (hist_id_stomate_IPCC, "rGrowth", itime, &
1824            vartmp, npts, hori_index)
1825       vartmp(:)=SUM(resp_maint*veget_cov_max,dim=2)/1e3/one_day
1826       CALL histwrite_p (hist_id_stomate_IPCC, "rMaint", itime, &
1827            vartmp, npts, hori_index)
1828       vartmp(:)=SUM(bm_alloc(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3/one_day
1829       CALL histwrite_p (hist_id_stomate_IPCC, "nppLeaf", itime, &
1830            vartmp, npts, hori_index)
1831       vartmp(:)=SUM((bm_alloc(:,:,isapabove,icarbon) + bm_alloc(:,:,iheartabove,icarbon))*veget_cov_max,dim=2)/1e3/one_day
1832       CALL histwrite_p (hist_id_stomate_IPCC, "nppStem", itime, &
1833            vartmp, npts, hori_index)
1834       vartmp(:)=SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )*veget_cov_max,dim=2)/1e3/one_day
1835       CALL histwrite_p (hist_id_stomate_IPCC, "nppRoot", itime, &
1836            vartmp, npts, hori_index)
1837
1838       CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_X', itime, &
1839            resolution(:,1), npts, hori_index)
1840       CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_Y', itime, &
1841            resolution(:,2), npts, hori_index)
1842       CALL histwrite_p (hist_id_stomate_IPCC, 'CONTFRAC', itime, &
1843            contfrac(:), npts, hori_index)
1844
1845    ENDIF
1846
1847    IF (printlev>=4) WRITE(numout,*) 'Leaving stomate_lpj'
1848
1849  END SUBROUTINE StomateLpj
1850
1851
1852!! ================================================================================================================================
1853!! SUBROUTINE   : harvest
1854!!
1855!>\BRIEF        Harvest of croplands
1856!!
1857!! DESCRIPTION  : To take into account biomass harvest from crop (mainly to take
1858!! into account for the reduced litter input and then decreased soil carbon. it is a
1859!! constant (40\%) fraction of above ground biomass.
1860!!
1861!! RECENT CHANGE(S) : None
1862!!
1863!! MAIN OUTPUT VARIABLE(S): ::harvest_above the harvested biomass
1864!!
1865!! REFERENCE(S) :
1866!! - Piao, S., P. Ciais, P. Friedlingstein, N. de Noblet-Ducoudre, P. Cadule, N. Viovy, and T. Wang. 2009.
1867!!   Spatiotemporal patterns of terrestrial carbon cycle during the 20th century. Global Biogeochemical
1868!!   Cycles 23:doi:10.1029/2008GB003339.
1869!!
1870!! FLOWCHART    : None
1871!! \n
1872!_ ================================================================================================================================
1873
1874  SUBROUTINE harvest(npts, dt_days, veget_cov_max, &
1875       bm_to_litter, turnover_daily, &
1876       harvest_above)
1877
1878  !! 0. Variable and parameter declaration
1879
1880    !! 0.1 Input variables
1881
1882    INTEGER, INTENT(in)                                    :: npts            !! Domain size (unitless)
1883    REAL(r_std), INTENT(in)                                :: dt_days         !! Time step (days)                               
1884    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: veget_cov_max       !! new "maximal" coverage fraction of a PFT (LAI ->
1885                                                                              !! infinity) on ground @tex $(m^2 m^{-2})$ @endtex
1886   
1887   !! 0.2 Output variables
1888   
1889   !! 0.3 Modified variables
1890
1891    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! [DISPENSABLE] conversion of biomass to litter
1892                                                                                     !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1893    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily   !! Turnover rates
1894                                                                               !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1895    REAL(r_std), DIMENSION(npts), INTENT(inout)            :: harvest_above    !! harvest above ground biomass for agriculture
1896                                                                               !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1897    !! 0.4 Local variables
1898
1899    INTEGER(i_std)                                         :: i, j, k, l, m    !! indices                       
1900    REAL(r_std)                                            :: above_old        !! biomass of previous time step
1901                                                                               !! @tex $(gC m^{-2})$ @endtex
1902!_ ================================================================================================================================
1903
1904  !! 1. Yearly initialisation
1905
1906    above_old             = zero
1907    harvest_above         = zero
1908
1909    DO i = 1, npts
1910       DO j = 1,nvm
1911          IF (.NOT. natural(j)) THEN
1912             above_old = turnover_daily(i,j,ileaf,icarbon) + turnover_daily(i,j,isapabove,icarbon) + &
1913                  &       turnover_daily(i,j,iheartabove,icarbon) + turnover_daily(i,j,ifruit,icarbon) + &
1914                  &       turnover_daily(i,j,icarbres,icarbon) + turnover_daily(i,j,isapbelow,icarbon) + &
1915                  &       turnover_daily(i,j,iheartbelow,icarbon) + turnover_daily(i,j,iroot,icarbon)
1916
1917             turnover_daily(i,j,ileaf,icarbon) = turnover_daily(i,j,ileaf,icarbon)*frac_turnover_daily
1918             turnover_daily(i,j,isapabove,icarbon) = turnover_daily(i,j,isapabove,icarbon)*frac_turnover_daily
1919             turnover_daily(i,j,isapbelow,icarbon) = turnover_daily(i,j,isapbelow,icarbon)*frac_turnover_daily
1920             turnover_daily(i,j,iheartabove,icarbon) = turnover_daily(i,j,iheartabove,icarbon)*frac_turnover_daily
1921             turnover_daily(i,j,iheartbelow,icarbon) = turnover_daily(i,j,iheartbelow,icarbon)*frac_turnover_daily
1922             turnover_daily(i,j,iroot,icarbon) = turnover_daily(i,j,iroot,icarbon)*frac_turnover_daily
1923             turnover_daily(i,j,ifruit,icarbon) = turnover_daily(i,j,ifruit,icarbon)*frac_turnover_daily
1924             turnover_daily(i,j,icarbres,icarbon) = turnover_daily(i,j,icarbres,icarbon)*frac_turnover_daily
1925             harvest_above(i)  = harvest_above(i) + veget_cov_max(i,j) * above_old *(un - frac_turnover_daily)
1926          ENDIF
1927       ENDDO
1928    ENDDO
1929
1930!!$    harvest_above = harvest_above
1931  END SUBROUTINE harvest
1932END MODULE stomate_lpj
Note: See TracBrowser for help on using the repository browser.