source: branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_lpj.f90 @ 7326

Last change on this file since 7326 was 7326, checked in by josefine.ghattas, 3 years ago

Corrected bug on carbon balance closure. See ticket #785
Integration in branch 2_2 done by P. Cadule

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