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

Last change on this file since 8543 was 8543, checked in by josefine.ghattas, 4 weeks ago

Integrated changset [8542] : Removed output variable CONTFRAC_STOMATE which was redundant with contfrac. Now in file_def, in stomate files use the same variable output from sechiba. Add a comment for RESOLUTION_X and RESOLUTION_Y to remember that these variables are deactivated automatically when running on an unstructured grid.
  • Property svn:keywords set to HeadURL Date Author Revision
File size: 102.1 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("T2M_MONTH",t2m_month)
948    CALL xios_orchidee_send_field("T2M_WEEK",t2m_week)
949    CALL xios_orchidee_send_field("TSEASON",Tseason)
950    CALL xios_orchidee_send_field("TMIN_SPRING_TIME", Tmin_spring_time)
951    CALL xios_orchidee_send_field("ONSET_DATE",onset_date)
952    CALL xios_orchidee_send_field("FPC_MAX",fpc_max)
953    CALL xios_orchidee_send_field("MAXFPC_LASTYEAR",maxfpc_lastyear)
954    CALL xios_orchidee_send_field("HET_RESP",resp_hetero(:,:))
955    CALL xios_orchidee_send_field("CO2_FIRE",co2_fire)
956    CALL xios_orchidee_send_field("CO2_TAKEN",co2_to_bm)
957    CALL xios_orchidee_send_field("LAI",lai)
958    CALL xios_orchidee_send_field("VEGET_COV_MAX",veget_cov_max)
959    CALL xios_orchidee_send_field("NPP_STOMATE",npp_daily+co2_to_bm)
960    CALL xios_orchidee_send_field("GPP",gpp_daily+co2_to_bm)
961    CALL xios_orchidee_send_field("IND",ind)
962    CALL xios_orchidee_send_field("CN_IND",cn_ind)
963    CALL xios_orchidee_send_field("WOODMASS_IND",woodmass_ind)
964    CALL xios_orchidee_send_field("TOTAL_M",tot_live_biomass)
965    CALL xios_orchidee_send_field("MOISTRESS",moiavail_week)
966    CALL xios_orchidee_send_field("LEAF_M",biomass(:,:,ileaf,icarbon))
967    CALL xios_orchidee_send_field("SAP_M_AB",biomass(:,:,isapabove,icarbon))
968    CALL xios_orchidee_send_field("SAP_M_BE",biomass(:,:,isapbelow,icarbon))
969    CALL xios_orchidee_send_field("HEART_M_AB",biomass(:,:,iheartabove,icarbon))
970    CALL xios_orchidee_send_field("HEART_M_BE",biomass(:,:,iheartbelow,icarbon))
971    CALL xios_orchidee_send_field("ROOT_M",biomass(:,:,iroot,icarbon))
972    CALL xios_orchidee_send_field("FRUIT_M",biomass(:,:,ifruit,icarbon))
973    CALL xios_orchidee_send_field("RESERVE_M",biomass(:,:,icarbres,icarbon))
974    CALL xios_orchidee_send_field("TOTAL_TURN",tot_turnover)
975    CALL xios_orchidee_send_field("LEAF_TURN",turnover_daily(:,:,ileaf,icarbon))
976    CALL xios_orchidee_send_field("MAINT_RESP",resp_maint)
977    CALL xios_orchidee_send_field("GROWTH_RESP",resp_growth)
978    CALL xios_orchidee_send_field("SAP_AB_TURN",turnover_daily(:,:,isapabove,icarbon))
979    CALL xios_orchidee_send_field("ROOT_TURN",turnover_daily(:,:,iroot,icarbon))
980    CALL xios_orchidee_send_field("FRUIT_TURN",turnover_daily(:,:,ifruit,icarbon))
981    CALL xios_orchidee_send_field("TOTAL_BM_LITTER",tot_bm_to_litter(:,:,icarbon))
982    CALL xios_orchidee_send_field("LEAF_BM_LITTER",bm_to_litter(:,:,ileaf,icarbon))
983    CALL xios_orchidee_send_field("SAP_AB_BM_LITTER",bm_to_litter(:,:,isapabove,icarbon))
984    CALL xios_orchidee_send_field("SAP_BE_BM_LITTER",bm_to_litter(:,:,isapbelow,icarbon))
985    CALL xios_orchidee_send_field("HEART_AB_BM_LITTER",bm_to_litter(:,:,iheartabove,icarbon))
986    CALL xios_orchidee_send_field("HEART_BE_BM_LITTER",bm_to_litter(:,:,iheartbelow,icarbon))
987    CALL xios_orchidee_send_field("ROOT_BM_LITTER",bm_to_litter(:,:,iroot,icarbon))
988    CALL xios_orchidee_send_field("FRUIT_BM_LITTER",bm_to_litter(:,:,ifruit,icarbon))
989    CALL xios_orchidee_send_field("RESERVE_BM_LITTER",bm_to_litter(:,:,icarbres,icarbon))
990    CALL xios_orchidee_send_field("LITTER_STR_AB",litter(:,istructural,:,iabove,icarbon))
991    CALL xios_orchidee_send_field("LITTER_MET_AB",litter(:,imetabolic,:,iabove,icarbon))
992    CALL xios_orchidee_send_field("LITTER_STR_BE",litter(:,istructural,:,ibelow,icarbon))
993    CALL xios_orchidee_send_field("LITTER_MET_BE",litter(:,imetabolic,:,ibelow,icarbon))
994    CALL xios_orchidee_send_field("DEADLEAF_COVER",deadleaf_cover)
995    CALL xios_orchidee_send_field("TOTAL_SOIL_CARB",tot_litter_soil_carb)
996    CALL xios_orchidee_send_field("CARBON_ACTIVE",carbon(:,iactive,:))
997    CALL xios_orchidee_send_field("CARBON_SLOW",carbon(:,islow,:))
998    CALL xios_orchidee_send_field("CARBON_PASSIVE",carbon(:,ipassive,:))
999    CALL xios_orchidee_send_field("LITTERHUM",litterhum_daily)
1000    CALL xios_orchidee_send_field("TURNOVER_TIME",turnover_time)
1001    CALL xios_orchidee_send_field("PROD10",prod10)
1002    CALL xios_orchidee_send_field("FLUX10",flux10)
1003    CALL xios_orchidee_send_field("PROD100",prod100)
1004    CALL xios_orchidee_send_field("FLUX100",flux100)
1005    CALL xios_orchidee_send_field("CONVFLUX",convflux)
1006    CALL xios_orchidee_send_field("CFLUX_PROD10",cflux_prod10)
1007    CALL xios_orchidee_send_field("CFLUX_PROD100",cflux_prod100)
1008    CALL xios_orchidee_send_field("PROD10_HARVEST",prod10_harvest)
1009    CALL xios_orchidee_send_field("FLUX10_HARVEST",flux10_harvest)
1010    CALL xios_orchidee_send_field("PROD100_HARVEST",prod100_harvest)
1011    CALL xios_orchidee_send_field("FLUX100_HARVEST",flux100_harvest)
1012    CALL xios_orchidee_send_field("CONVFLUX_HARVEST",convflux_harvest)
1013    CALL xios_orchidee_send_field("CFLUX_PROD10_HARVEST",cflux_prod10_harvest)
1014    CALL xios_orchidee_send_field("CFLUX_PROD100_HARVEST",cflux_prod100_harvest)
1015    IF (do_wood_harvest) THEN
1016       CALL xios_orchidee_send_field("WOOD_HARVEST",woodharvest/one_year*dt_days)
1017       CALL xios_orchidee_send_field("WOOD_HARVEST_PFT",woodharvestpft)
1018    END IF
1019    CALL xios_orchidee_send_field("HARVEST_ABOVE",harvest_above)
1020    CALL xios_orchidee_send_field("VCMAX",vcmax)
1021    CALL xios_orchidee_send_field("AGE",age)
1022    CALL xios_orchidee_send_field("HEIGHT",height)
1023    CALL xios_orchidee_send_field("FIREINDEX",fireindex(:,:))
1024
1025    ! ipcc history
1026    ! Carbon stock transformed from gC/m2 into kgC/m2
1027    CALL xios_orchidee_send_field("cVeg",SUM(tot_live_biomass(:,:,icarbon)*veget_cov_max,dim=2)/1e3)
1028    CALL xios_orchidee_send_field("cVegGrass",sum_cVegGrass/1e3)
1029    CALL xios_orchidee_send_field("cVegCrop",sum_cVegCrop/1e3)
1030    CALL xios_orchidee_send_field("cVegTree",sum_cVegTree/1e3)
1031    CALL xios_orchidee_send_field("cOther",SUM(cOther(:,:,icarbon)*veget_cov_max,dim=2)/1e3)
1032    CALL xios_orchidee_send_field("cLitter",SUM(tot_litter_carb*veget_cov_max,dim=2)/1e3)
1033    CALL xios_orchidee_send_field("cLitterGrass",sum_cLitterGrass/1e3)
1034    CALL xios_orchidee_send_field("cLitterCrop",sum_cLitterCrop/1e3)
1035    CALL xios_orchidee_send_field("cLitterTree",sum_cLitterTree/1e3)
1036    CALL xios_orchidee_send_field("cSoil",SUM(tot_soil_carb*veget_cov_max,dim=2)/1e3)
1037    CALL xios_orchidee_send_field("cSoilGrass",sum_cSoilGrass/1e3)
1038    CALL xios_orchidee_send_field("cSoilCrop",sum_cSoilCrop/1e3)
1039    CALL xios_orchidee_send_field("cSoilTree",sum_cSoilTree/1e3)
1040    CALL xios_orchidee_send_field("cProduct",(prod10_total + prod100_total + prod10_harvest_total + prod100_harvest_total)/1e3)
1041
1042    cproductlut(:,:)=xios_default_val
1043    WHERE(fraclut(:,id_psl) > min_sechiba)
1044       cproductlut(:,id_psl)= ( prod10_total + prod100_total + prod10_harvest_total + prod100_harvest_total)/1e3 &
1045            / fraclut(:,id_psl)
1046    ENDWHERE
1047    CALL xios_orchidee_send_field("cproductlut",cproductlut)
1048
1049    CALL xios_orchidee_send_field("cMassVariation",carb_mass_variation/1e3/one_day)
1050
1051    CALL xios_orchidee_send_field("lai_ipcc",SUM(lai*veget_cov_max,dim=2)) ! m2/m2
1052   
1053    ! Carbon fluxes transformed from gC/m2/d into kgC/m2/s
1054    CALL xios_orchidee_send_field("gpp_ipcc",SUM((gpp_daily+co2_to_bm)*veget_cov_max,dim=2)/1e3/one_day)
1055    CALL xios_orchidee_send_field("ra",SUM((resp_maint+resp_growth)*veget_cov_max,dim=2)/1e3/one_day)
1056    vartmp(:)=zero
1057    DO j = 2, nvm
1058       IF ( .NOT. is_tree(j) .AND. natural(j) ) THEN
1059          vartmp(:) = vartmp(:) + (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)
1060       ENDIF
1061    ENDDO
1062    CALL xios_orchidee_send_field("raGrass",vartmp/1e3/one_day)
1063
1064    vartmp(:)=zero
1065    DO j = 2, nvm
1066       IF (( .NOT. is_tree(j)) .AND. (.NOT. natural(j)) ) THEN
1067          vartmp(:) = vartmp(:) + (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)
1068       ENDIF
1069    ENDDO
1070    CALL xios_orchidee_send_field("raCrop",vartmp/1e3/one_day)
1071
1072    vartmp(:)=zero
1073    DO j = 2, nvm
1074       IF ( is_tree(j) ) THEN
1075          vartmp(:) = vartmp(:) + (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)
1076       ENDIF
1077    ENDDO
1078    CALL xios_orchidee_send_field("raTree",vartmp/1e3/one_day)
1079
1080    CALL xios_orchidee_send_field("npp_ipcc",SUM((npp_daily+co2_to_bm)*veget_cov_max,dim=2)/1e3/one_day)
1081
1082    vartmp(:)=zero
1083    DO j = 2, nvm
1084       IF ( .NOT. is_tree(j) .AND. natural(j) ) THEN
1085          vartmp(:) = vartmp(:) +( npp_daily(:,j)+co2_to_bm(:,j))*veget_cov_max(:,j)
1086       ENDIF
1087    ENDDO
1088    CALL xios_orchidee_send_field("nppGrass",vartmp/1e3/one_day)
1089
1090    DO j = 2, nvm
1091       IF ( (.NOT. is_tree(j)) .AND. (.NOT. natural(j)) ) THEN
1092          vartmp(:) = vartmp(:) + (npp_daily(:,j)+co2_to_bm(:,j))*veget_cov_max(:,j)
1093       ENDIF
1094    ENDDO
1095    CALL xios_orchidee_send_field("nppCrop",vartmp/1e3/one_day)
1096
1097    vartmp(:)=zero
1098    DO j = 2, nvm
1099       IF ( is_tree(j) ) THEN
1100          vartmp(:) = vartmp(:) + (npp_daily(:,j)+co2_to_bm(:,j))*veget_cov_max(:,j)
1101       ENDIF
1102    ENDDO
1103    CALL xios_orchidee_send_field("nppTree",vartmp/1e3/one_day)
1104
1105    CALL xios_orchidee_send_field("rh",SUM(resp_hetero*veget_cov_max,dim=2)/1e3/one_day)
1106    CALL xios_orchidee_send_field("rhLitter",SUM(resp_hetero_litter*veget_cov_max,dim=2)/1e3/one_day)
1107    CALL xios_orchidee_send_field("rhSoil",SUM(resp_hetero_soil*veget_cov_max,dim=2)/1e3/one_day)
1108    CALL xios_orchidee_send_field("HET_RESP_SOIL",resp_hetero_soil(:,:))
1109    CALL xios_orchidee_send_field("HET_RESP_LITTER",resp_hetero_litter(:,:))
1110
1111    CALL xios_orchidee_send_field("fFire",SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day)
1112    ctotfirelut(:,:)=xios_default_val
1113    WHERE(fraclut(:,id_psl) > min_sechiba) 
1114       ctotfirelut(:,id_psl)= SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day &
1115            / fraclut(:,id_psl)
1116    ENDWHERE
1117    CALL xios_orchidee_send_field("ctotfirelut",ctotfirelut)
1118
1119    CALL xios_orchidee_send_field("fHarvest",harvest_above/1e3/one_day)
1120    Call xios_orchidee_send_field("fLuc",cflux_prod_total/1e3/one_day)
1121    CALL xios_orchidee_send_field("fWoodharvest",cflux_prod_harvest_total/1e3/one_day)
1122    CALL xios_orchidee_send_field("fDeforestToProduct",SUM((fDeforestToProduct-convfluxpft),dim=2)/1e3/one_day)
1123   
1124    flulccproductlut(:,:)=xios_default_val
1125    flulccproductlut(:,id_psl)=0.
1126    flulccproductlut(:,id_crp)=0.
1127    DO j=1,nvm
1128       IF(natural(j)) THEN
1129          flulccproductlut(:,id_psl)=flulccproductlut(:,id_psl) + fDeforestToProduct(:,j)
1130       ELSE
1131          flulccproductlut(:,id_crp)=flulccproductlut(:,id_crp) + fDeforestToProduct(:,j)
1132       ENDIF
1133    ENDDO
1134
1135    WHERE(fraclut(:,id_psl) > min_sechiba)
1136       flulccproductlut(:,id_psl)= (flulccproductlut(:,id_psl) &
1137            +SUM(fHarvestToProduct,dim=2)) &
1138            /1e3/one_day / fraclut(:,id_psl)
1139    ELSEWHERE
1140       flulccproductlut(:,id_psl)= xios_default_val
1141    ENDWHERE
1142    WHERE(fraclut(:,id_crp) > min_sechiba)
1143       flulccproductlut(:,id_crp)= flulccproductlut(:,id_crp) &
1144            /1e3/one_day / fraclut(:,id_crp)
1145    ELSEWHERE
1146       flulccproductlut(:,id_crp)= xios_default_val
1147    ENDWHERE
1148    CALL xios_orchidee_send_field("flulccproductlut",flulccproductlut)
1149
1150    flulccresiduelut(:,:)=xios_default_val
1151    flulccresiduelut(:,id_psl)=0.
1152    flulccresiduelut(:,id_crp)=0.
1153    DO j=1,nvm
1154       IF(natural(j)) THEN
1155          flulccresiduelut(:,id_psl) = flulccresiduelut(:,id_psl) + fLulccResidue(:,j)
1156       ELSE
1157          flulccresiduelut(:,id_crp) = flulccresiduelut(:,id_crp) + fLulccResidue(:,j)
1158       ENDIF
1159    ENDDO
1160
1161    WHERE(fraclut(:,id_psl) > min_sechiba)
1162       flulccresiduelut(:,id_psl)= flulccresiduelut(:,id_psl)/1e3/one_day/fraclut(:,id_psl)
1163    ELSEWHERE
1164       flulccresiduelut(:,id_psl)= xios_default_val
1165    ENDWHERE
1166    WHERE(fraclut(:,id_crp) > min_sechiba)
1167       flulccresiduelut(:,id_crp)= flulccresiduelut(:,id_crp)/1e3/one_day /fraclut(:,id_crp)
1168    ELSEWHERE
1169       flulccresiduelut(:,id_crp)= xios_default_val
1170    ENDWHERE
1171    CALL xios_orchidee_send_field("flulccresiduelut",flulccresiduelut)
1172    CALL xios_orchidee_send_field("flulccresidue",flulccresidue)
1173
1174    flulccatmlut(:,:)=xios_default_val
1175    flulccatmlut(:,id_psl)=0.
1176    flulccatmlut(:,id_crp)=0.
1177    DO j=1,nvm
1178       IF(natural(j)) THEN
1179          flulccatmlut(:,id_psl)=flulccatmlut(:,id_psl)+convfluxpft(:,j)
1180       ELSE
1181          flulccatmlut(:,id_crp)=flulccatmlut(:,id_crp)+convfluxpft(:,j)
1182       ENDIF
1183    ENDDO
1184
1185    WHERE(fraclut(:,id_psl) > min_sechiba)
1186       flulccatmlut(:,id_psl)= (flulccatmlut(:,id_psl)+cflux_prod10+cflux_prod100+cflux_prod_harvest_total) &
1187            /1e3/one_day / fraclut(:,id_psl)
1188    ELSEWHERE
1189       flulccatmlut(:,id_psl)= xios_default_val
1190    ENDWHERE
1191    WHERE(fraclut(:,id_crp) > min_sechiba)
1192       flulccatmlut(:,id_crp)= (flulccatmlut(:,id_crp)+harvest_above)/1e3/one_day / fraclut(:,id_crp)
1193    ELSEWHERE
1194       flulccatmlut(:,id_crp)= xios_default_val
1195    ENDWHERE
1196    CALL xios_orchidee_send_field("flulccatmlut",flulccatmlut)
1197
1198    CALL xios_orchidee_send_field("nbp",(SUM((gpp_daily+co2_to_bm-(resp_maint+resp_growth+resp_hetero)-co2_fire) * &
1199          veget_cov_max,dim=2)-cflux_prod_total-cflux_prod_harvest_total-harvest_above)/1e3/one_day)
1200    CALL xios_orchidee_send_field("fVegLitter",SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*&
1201         veget_cov_max,dim=2)/1e3/one_day)
1202    CALL xios_orchidee_send_field("fLitterSoil",SUM(SUM(soilcarbon_input,dim=2)*veget_cov_max,dim=2)/1e3/one_day)
1203
1204    ! Carbon stock transformed from gC/m2 into kgC/m2
1205    CALL xios_orchidee_send_field("cLeaf",SUM(biomass(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3)
1206    CALL xios_orchidee_send_field("cStem",SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*&
1207          veget_cov_max,dim=2)/1e3)
1208    CALL xios_orchidee_send_field("cWood",SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon) + &
1209         biomass(:,:,isapbelow,icarbon) + biomass(:,:,iheartbelow,icarbon))* &
1210         veget_cov_max,dim=2)/1e3)
1211    CALL xios_orchidee_send_field("cRoot",SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + &
1212         biomass(:,:,iheartbelow,icarbon) )*veget_cov_max,dim=2)/1e3)
1213    CALL xios_orchidee_send_field("cMisc",SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*&
1214         veget_cov_max,dim=2)/1e3)
1215    CALL xios_orchidee_send_field("cLitterAbove",SUM((litter(:,istructural,:,iabove,icarbon)+&
1216         litter(:,imetabolic,:,iabove,icarbon))*veget_cov_max,dim=2)/1e3)
1217    CALL xios_orchidee_send_field("cLitterBelow",SUM((litter(:,istructural,:,ibelow,icarbon)+&
1218         litter(:,imetabolic,:,ibelow,icarbon))*veget_cov_max,dim=2)/1e3)
1219    CALL xios_orchidee_send_field("cSoilFast",SUM(carbon(:,iactive,:)*veget_cov_max,dim=2)/1e3)
1220    CALL xios_orchidee_send_field("cSoilMedium",SUM(carbon(:,islow,:)*veget_cov_max,dim=2)/1e3)
1221    CALL xios_orchidee_send_field("cSoilSlow",SUM(carbon(:,ipassive,:)*veget_cov_max,dim=2)/1e3)
1222
1223    DO k = 1, ncarb
1224       csoilpools(:,k) = SUM(carbon(:,k,:)*veget_cov_max,dim=2)/1e3
1225    END DO
1226    CALL xios_orchidee_send_field("cSoilPools",csoilpools)
1227
1228    ! Vegetation fractions [0,100]
1229    CALL xios_orchidee_send_field("landCoverFrac",veget_cov_max*100)
1230
1231    ! Carbon fluxes transformed from gC/m2/d into kgC/m2/s
1232    CALL xios_orchidee_send_field("rGrowth",SUM(resp_growth*veget_cov_max,dim=2)/1e3/one_day)
1233    CALL xios_orchidee_send_field("rMaint",SUM(resp_maint*veget_cov_max,dim=2)/1e3/one_day)
1234    CALL xios_orchidee_send_field("nppLeaf",SUM(bm_alloc(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3/one_day)
1235    ! nppStem : from wood above surface
1236    CALL xios_orchidee_send_field("nppStem",SUM((bm_alloc(:,:,isapabove,icarbon) + bm_alloc(:,:,iheartabove,icarbon)) &
1237         *veget_cov_max,dim=2)/1e3/one_day)
1238    ! nppWood : from wood above and below surface
1239    CALL xios_orchidee_send_field("nppWood",SUM((bm_alloc(:,:,isapabove,icarbon) + bm_alloc(:,:,iheartabove,icarbon) + &
1240         bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iheartbelow,icarbon)) &
1241         *veget_cov_max,dim=2)/1e3/one_day)
1242    ! nppRoot : from wood below surface and fine roots
1243    CALL xios_orchidee_send_field("nppRoot",SUM((bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iheartbelow,icarbon) + &
1244         bm_alloc(:,:,iroot,icarbon)) * veget_cov_max,dim=2)/1e3/one_day)
1245    CALL xios_orchidee_send_field("nppOther",SUM(( bm_alloc(:,:,ifruit,icarbon) + bm_alloc(:,:,icarbres,icarbon) ) * &
1246         veget_cov_max,dim=2)/1e3/one_day)
1247
1248    ! Calculate variables according to LUMIP specifications.
1249    clitterlut(:,:) = 0 
1250    csoillut(:,:) = 0 
1251    cveglut(:,:) = 0
1252    lailut(:,:) = 0
1253    ralut(:,:) = 0
1254    rhlut(:,:) = 0
1255    npplut(:,:) = 0
1256
1257    DO j=1,nvm
1258       IF (natural(j)) THEN
1259          clitterlut(:,id_psl) = clitterlut(:,id_psl) + tot_litter_carb(:,j)*veget_cov_max(:,j)/1e3
1260          csoillut(:,id_psl) = csoillut(:,id_psl) + tot_soil_carb(:,j)*veget_cov_max(:,j)/1e3
1261          cveglut(:,id_psl) = cveglut(:,id_psl) + tot_live_biomass(:,j,icarbon)*veget_cov_max(:,j)/1e3
1262          lailut(:,id_psl) = lailut(:,id_psl) + lai(:,j)*veget_cov_max(:,j)
1263          ralut(:,id_psl) = ralut(:,id_psl) + &
1264               (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)/1e3/one_day
1265          rhlut(:,id_psl) = rhlut(:,id_psl) + &
1266               resp_hetero(:,j)*veget_cov_max(:,j)/1e3/one_day
1267          npplut(:,id_psl) = npplut(:,id_psl) + &
1268               (npp_daily(:,j)+co2_to_bm(:,j))*veget_cov_max(:,j)/1e3/one_day
1269       ELSE
1270          clitterlut(:,id_crp) = clitterlut(:,id_crp) + tot_litter_carb(:,j)*veget_cov_max(:,j)/1e3
1271          csoillut(:,id_crp) = csoillut(:,id_crp) + tot_soil_carb(:,j)*veget_cov_max(:,j)/1e3
1272          cveglut(:,id_crp) = cveglut(:,id_crp) + tot_live_biomass(:,j,icarbon)*veget_cov_max(:,j)/1e3
1273          lailut(:,id_crp) = lailut(:,id_crp) + lai(:,j)*veget_cov_max(:,j)
1274          ralut(:,id_crp) = ralut(:,id_crp) + &
1275               (resp_maint(:,j)+resp_growth(:,j))*veget_cov_max(:,j)/1e3/one_day
1276          rhlut(:,id_crp) = rhlut(:,id_crp) + &
1277               resp_hetero(:,j)*veget_cov_max(:,j)/1e3/one_day
1278          npplut(:,id_crp) = npplut(:,id_crp) + &
1279               (npp_daily(:,j)+co2_to_bm(:,j))*veget_cov_max(:,j)/1e3/one_day
1280       END IF
1281    END DO
1282
1283    WHERE (fraclut(:,id_psl)>min_sechiba)
1284       clitterlut(:,id_psl) = clitterlut(:,id_psl)/fraclut(:,id_psl)
1285       csoillut(:,id_psl) = csoillut(:,id_psl)/fraclut(:,id_psl)
1286       cveglut(:,id_psl) = cveglut(:,id_psl)/fraclut(:,id_psl)
1287       lailut(:,id_psl) = lailut(:,id_psl)/fraclut(:,id_psl)
1288       ralut(:,id_psl) = ralut(:,id_psl)/fraclut(:,id_psl)
1289       rhlut(:,id_psl) = rhlut(:,id_psl)/fraclut(:,id_psl)
1290       npplut(:,id_psl) = npplut(:,id_psl)/fraclut(:,id_psl)
1291    ELSEWHERE
1292       clitterlut(:,id_psl) = xios_default_val
1293       csoillut(:,id_psl) = xios_default_val
1294       cveglut(:,id_psl) = xios_default_val
1295       lailut(:,id_psl) = xios_default_val
1296       ralut(:,id_psl) = xios_default_val
1297       rhlut(:,id_psl) = xios_default_val
1298       npplut(:,id_psl) = xios_default_val
1299    END WHERE
1300
1301    WHERE (fraclut(:,id_crp)>min_sechiba)
1302       clitterlut(:,id_crp) = clitterlut(:,id_crp)/fraclut(:,id_crp)
1303       csoillut(:,id_crp) = csoillut(:,id_crp)/fraclut(:,id_crp)
1304       cveglut(:,id_crp) = cveglut(:,id_crp)/fraclut(:,id_crp)
1305       lailut(:,id_crp) = lailut(:,id_crp)/fraclut(:,id_crp)
1306       ralut(:,id_crp) = ralut(:,id_crp)/fraclut(:,id_crp)
1307       rhlut(:,id_crp) = rhlut(:,id_crp)/fraclut(:,id_crp)
1308       npplut(:,id_crp) = npplut(:,id_crp)/fraclut(:,id_crp)
1309    ELSEWHERE
1310       clitterlut(:,id_crp) = xios_default_val
1311       csoillut(:,id_crp) = xios_default_val
1312       cveglut(:,id_crp) = xios_default_val
1313       lailut(:,id_crp) = xios_default_val
1314       ralut(:,id_crp) = xios_default_val
1315       rhlut(:,id_crp) = xios_default_val
1316       npplut(:,id_crp) = xios_default_val
1317    END WHERE
1318
1319    clitterlut(:,id_pst) = xios_default_val
1320    clitterlut(:,id_urb) = xios_default_val
1321    csoillut(:,id_pst)   = xios_default_val
1322    csoillut(:,id_urb)   = xios_default_val
1323    cveglut(:,id_pst)    = xios_default_val
1324    cveglut(:,id_urb)    = xios_default_val
1325    lailut(:,id_pst)     = xios_default_val
1326    lailut(:,id_urb)     = xios_default_val
1327    ralut(:,id_pst)      = xios_default_val
1328    ralut(:,id_urb)      = xios_default_val
1329    rhlut(:,id_pst)      = xios_default_val
1330    rhlut(:,id_urb)      = xios_default_val
1331    npplut(:,id_pst)     = xios_default_val
1332    npplut(:,id_urb)     = xios_default_val
1333
1334    CALL xios_orchidee_send_field("clitterlut",clitterlut)
1335    CALL xios_orchidee_send_field("csoillut",csoillut)
1336    CALL xios_orchidee_send_field("cveglut",cveglut)
1337    CALL xios_orchidee_send_field("lailut",lailut)
1338    CALL xios_orchidee_send_field("ralut",ralut)
1339    CALL xios_orchidee_send_field("rhlut",rhlut)
1340    CALL xios_orchidee_send_field("npplut",npplut)
1341
1342    CALL histwrite_p (hist_id_stomate, 'RESOLUTION_X', itime, &
1343         resolution(:,1), npts, hori_index)
1344    CALL histwrite_p (hist_id_stomate, 'RESOLUTION_Y', itime, &
1345         resolution(:,2), npts, hori_index)
1346    CALL histwrite_p (hist_id_stomate, 'CONTFRAC', itime, &
1347         contfrac(:), npts, hori_index)
1348
1349    CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AB', itime, &
1350         litter(:,istructural,:,iabove,icarbon), npts*nvm, horipft_index)
1351    CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AB', itime, &
1352         litter(:,imetabolic,:,iabove,icarbon), npts*nvm, horipft_index)
1353    CALL histwrite_p (hist_id_stomate, 'LITTER_STR_BE', itime, &
1354         litter(:,istructural,:,ibelow,icarbon), npts*nvm, horipft_index)
1355    CALL histwrite_p (hist_id_stomate, 'LITTER_MET_BE', itime, &
1356         litter(:,imetabolic,:,ibelow,icarbon), npts*nvm, horipft_index)
1357
1358    CALL histwrite_p (hist_id_stomate, 'DEADLEAF_COVER', itime, &
1359         deadleaf_cover, npts, hori_index)
1360
1361    CALL histwrite_p (hist_id_stomate, 'TOTAL_SOIL_CARB', itime, &
1362         tot_litter_soil_carb, npts*nvm, horipft_index)
1363    CALL histwrite_p (hist_id_stomate, 'CARBON_ACTIVE', itime, &
1364         carbon(:,iactive,:), npts*nvm, horipft_index)
1365    CALL histwrite_p (hist_id_stomate, 'CARBON_SLOW', itime, &
1366         carbon(:,islow,:), npts*nvm, horipft_index)
1367    CALL histwrite_p (hist_id_stomate, 'CARBON_PASSIVE', itime, &
1368         carbon(:,ipassive,:), npts*nvm, horipft_index)
1369
1370    CALL histwrite_p (hist_id_stomate, 'T2M_MONTH', itime, &
1371         t2m_month, npts, hori_index)
1372    CALL histwrite_p (hist_id_stomate, 'T2M_WEEK', itime, &
1373         t2m_week, npts, hori_index)
1374    CALL histwrite_p (hist_id_stomate, 'TSEASON', itime, &
1375         Tseason, npts, hori_index)
1376    CALL histwrite_p (hist_id_stomate, 'TMIN_SPRING_TIME', itime, &
1377         Tmin_spring_time, npts*nvm, horipft_index)
1378    CALL histwrite_p (hist_id_stomate, 'ONSET_DATE', itime, &
1379         onset_date(:,:), npts*nvm, horipft_index)
1380    CALL histwrite_p (hist_id_stomate, 'FPC_MAX', itime, &
1381         fpc_max, npts*nvm, horipft_index)
1382    CALL histwrite_p (hist_id_stomate, 'MAXFPC_LASTYEAR', itime, &
1383         maxfpc_lastyear, npts*nvm, horipft_index)
1384    CALL histwrite_p (hist_id_stomate, 'HET_RESP', itime, &
1385         resp_hetero(:,:), npts*nvm, horipft_index)
1386    CALL histwrite_p (hist_id_stomate, 'FIREINDEX', itime, &
1387         fireindex(:,:), npts*nvm, horipft_index)
1388    CALL histwrite_p (hist_id_stomate, 'LITTERHUM', itime, &
1389         litterhum_daily, npts, hori_index)
1390    CALL histwrite_p (hist_id_stomate, 'CO2_FIRE', itime, &
1391         co2_fire, npts*nvm, horipft_index)
1392    CALL histwrite_p (hist_id_stomate, 'CO2_TAKEN', itime, &
1393         co2_to_bm, npts*nvm, horipft_index)
1394    ! land cover change
1395    CALL histwrite_p (hist_id_stomate, 'CONVFLUX', itime, &
1396         convflux, npts, hori_index)
1397    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10', itime, &
1398         cflux_prod10, npts, hori_index)
1399    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100', itime, &
1400         cflux_prod100, npts, hori_index)
1401    CALL histwrite_p (hist_id_stomate, 'CONVFLUX_HARVEST', itime, &
1402         convflux_harvest, npts, hori_index)
1403    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10_HARVEST', itime, &
1404         cflux_prod10_harvest, npts, hori_index)
1405    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100_HARVES', itime, &
1406         cflux_prod100_harvest, npts, hori_index)
1407
1408    IF (do_wood_harvest) THEN
1409       CALL histwrite_p (hist_id_stomate, 'WOOD_HARVEST', itime, &
1410            woodharvest/one_year*dt_days, npts, hori_index)
1411       CALL histwrite_p (hist_id_stomate, 'WOOD_HARVEST_PFT', itime, &
1412            woodharvestpft, npts*nvm, horipft_index)
1413    END IF
1414
1415    CALL histwrite_p (hist_id_stomate, 'HARVEST_ABOVE', itime, &
1416         harvest_above, npts, hori_index)
1417
1418    CALL histwrite_p (hist_id_stomate, 'LAI', itime, &
1419         lai, npts*nvm, horipft_index)
1420    CALL histwrite_p (hist_id_stomate, 'VEGET_COV_MAX', itime, &
1421         veget_cov_max, npts*nvm, horipft_index)
1422    CALL histwrite_p (hist_id_stomate, 'NPP', itime, &
1423         npp_daily+co2_to_bm, npts*nvm, horipft_index)
1424    CALL histwrite_p (hist_id_stomate, 'GPP', itime, &
1425         gpp_daily+co2_to_bm, npts*nvm, horipft_index)
1426    CALL histwrite_p (hist_id_stomate, 'IND', itime, &
1427         ind, npts*nvm, horipft_index)
1428    CALL histwrite_p (hist_id_stomate, 'CN_IND', itime, &
1429         cn_ind, npts*nvm, horipft_index)
1430    CALL histwrite_p (hist_id_stomate, 'WOODMASS_IND', itime, &
1431         woodmass_ind, npts*nvm, horipft_index)
1432    CALL histwrite_p (hist_id_stomate, 'TOTAL_M', itime, &
1433         tot_live_biomass(:,:,icarbon), npts*nvm, horipft_index)
1434    CALL histwrite_p (hist_id_stomate, 'LEAF_M', itime, &
1435         biomass(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1436    CALL histwrite_p (hist_id_stomate, 'SAP_M_AB', itime, &
1437         biomass(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1438    CALL histwrite_p (hist_id_stomate, 'SAP_M_BE', itime, &
1439         biomass(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
1440    CALL histwrite_p (hist_id_stomate, 'HEART_M_AB', itime, &
1441         biomass(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
1442    CALL histwrite_p (hist_id_stomate, 'HEART_M_BE', itime, &
1443         biomass(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
1444    CALL histwrite_p (hist_id_stomate, 'ROOT_M', itime, &
1445         biomass(:,:,iroot,icarbon), npts*nvm, horipft_index)
1446    CALL histwrite_p (hist_id_stomate, 'FRUIT_M', itime, &
1447         biomass(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1448    CALL histwrite_p (hist_id_stomate, 'RESERVE_M', itime, &
1449         biomass(:,:,icarbres,icarbon), npts*nvm, horipft_index)
1450    CALL histwrite_p (hist_id_stomate, 'TOTAL_TURN', itime, &
1451         tot_turnover(:,:,icarbon), npts*nvm, horipft_index)
1452    CALL histwrite_p (hist_id_stomate, 'LEAF_TURN', itime, &
1453         turnover_daily(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1454    CALL histwrite_p (hist_id_stomate, 'SAP_AB_TURN', itime, &
1455         turnover_daily(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1456    CALL histwrite_p (hist_id_stomate, 'ROOT_TURN', itime, &
1457         turnover_daily(:,:,iroot,icarbon), npts*nvm, horipft_index)
1458    CALL histwrite_p (hist_id_stomate, 'FRUIT_TURN', itime, &
1459         turnover_daily(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1460    CALL histwrite_p (hist_id_stomate, 'TOTAL_BM_LITTER', itime, &
1461         tot_bm_to_litter(:,:,icarbon), npts*nvm, horipft_index)
1462    CALL histwrite_p (hist_id_stomate, 'LEAF_BM_LITTER', itime, &
1463         bm_to_litter(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1464    CALL histwrite_p (hist_id_stomate, 'SAP_AB_BM_LITTER', itime, &
1465         bm_to_litter(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1466    CALL histwrite_p (hist_id_stomate, 'SAP_BE_BM_LITTER', itime, &
1467         bm_to_litter(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
1468    CALL histwrite_p (hist_id_stomate, 'HEART_AB_BM_LITTER', itime, &
1469         bm_to_litter(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
1470    CALL histwrite_p (hist_id_stomate, 'HEART_BE_BM_LITTER', itime, &
1471         bm_to_litter(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
1472    CALL histwrite_p (hist_id_stomate, 'ROOT_BM_LITTER', itime, &
1473         bm_to_litter(:,:,iroot,icarbon), npts*nvm, horipft_index)
1474    CALL histwrite_p (hist_id_stomate, 'FRUIT_BM_LITTER', itime, &
1475         bm_to_litter(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1476    CALL histwrite_p (hist_id_stomate, 'RESERVE_BM_LITTER', itime, &
1477         bm_to_litter(:,:,icarbres,icarbon), npts*nvm, horipft_index)
1478    CALL histwrite_p (hist_id_stomate, 'MAINT_RESP', itime, &
1479         resp_maint, npts*nvm, horipft_index)
1480    CALL histwrite_p (hist_id_stomate, 'GROWTH_RESP', itime, &
1481         resp_growth, npts*nvm, horipft_index)
1482    CALL histwrite_p (hist_id_stomate, 'AGE', itime, &
1483         age, npts*nvm, horipft_index)
1484    CALL histwrite_p (hist_id_stomate, 'HEIGHT', itime, &
1485         height, npts*nvm, horipft_index)
1486    CALL histwrite_p (hist_id_stomate, 'MOISTRESS', itime, &
1487         moiavail_week, npts*nvm, horipft_index)
1488    CALL histwrite_p (hist_id_stomate, 'VCMAX', itime, &
1489         vcmax, npts*nvm, horipft_index)
1490    CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME', itime, &
1491         turnover_time, npts*nvm, horipft_index)
1492    ! land cover change
1493    CALL histwrite_p (hist_id_stomate, 'PROD10', itime, &
1494         prod10, npts*11, horip11_index)
1495    CALL histwrite_p (hist_id_stomate, 'PROD100', itime, &
1496         prod100, npts*101, horip101_index)
1497    CALL histwrite_p (hist_id_stomate, 'FLUX10', itime, &
1498         flux10, npts*10, horip10_index)
1499    CALL histwrite_p (hist_id_stomate, 'FLUX100', itime, &
1500         flux100, npts*100, horip100_index)
1501    CALL histwrite_p (hist_id_stomate, 'PROD10_HARVEST', itime, &
1502         prod10_harvest, npts*11, horip11_index)
1503    CALL histwrite_p (hist_id_stomate, 'PROD100_HARVEST', itime, &
1504         prod100_harvest, npts*101, horip101_index)
1505    CALL histwrite_p (hist_id_stomate, 'FLUX10_HARVEST', itime, &
1506         flux10_harvest, npts*10, horip10_index)
1507    CALL histwrite_p (hist_id_stomate, 'FLUX100_HARVEST', itime, &
1508         flux100_harvest, npts*100, horip100_index)
1509
1510    IF ( hist_id_stomate_IPCC > 0 ) THEN
1511       vartmp(:)=SUM(tot_live_biomass(:,:,icarbon)*veget_cov_max,dim=2)/1e3
1512       CALL histwrite_p (hist_id_stomate_IPCC, "cVeg", itime, &
1513            vartmp, npts, hori_index)
1514       vartmp(:)=SUM(tot_litter_carb*veget_cov_max,dim=2)/1e3
1515       CALL histwrite_p (hist_id_stomate_IPCC, "cLitter", itime, &
1516            vartmp, npts, hori_index)
1517       vartmp(:)=SUM(tot_soil_carb*veget_cov_max,dim=2)/1e3
1518       CALL histwrite_p (hist_id_stomate_IPCC, "cSoil", itime, &
1519            vartmp, npts, hori_index)
1520       vartmp(:)=(prod10_total + prod100_total + prod10_harvest_total + prod100_harvest_total)/1e3
1521       CALL histwrite_p (hist_id_stomate_IPCC, "cProduct", itime, &
1522            vartmp, npts, hori_index)
1523       vartmp(:)=carb_mass_variation/1e3/one_day
1524       CALL histwrite_p (hist_id_stomate_IPCC, "cMassVariation", itime, &
1525            vartmp, npts, hori_index)
1526       vartmp(:)=SUM(lai*veget_cov_max,dim=2)
1527       CALL histwrite_p (hist_id_stomate_IPCC, "lai", itime, &
1528            vartmp, npts, hori_index)
1529       vartmp(:)=SUM((gpp_daily+co2_to_bm)*veget_cov_max,dim=2)/1e3/one_day
1530       CALL histwrite_p (hist_id_stomate_IPCC, "gpp", itime, &
1531            vartmp, npts, hori_index)
1532       vartmp(:)=SUM((resp_maint+resp_growth)*veget_cov_max,dim=2)/1e3/one_day
1533       CALL histwrite_p (hist_id_stomate_IPCC, "ra", itime, &
1534            vartmp, npts, hori_index)
1535       vartmp(:)=SUM((npp_daily+co2_to_bm)*veget_cov_max,dim=2)/1e3/one_day
1536       CALL histwrite_p (hist_id_stomate_IPCC, "npp", itime, &
1537            vartmp, npts, hori_index)
1538       vartmp(:)=SUM(resp_hetero*veget_cov_max,dim=2)/1e3/one_day
1539       CALL histwrite_p (hist_id_stomate_IPCC, "rh", itime, &
1540            vartmp, npts, hori_index)
1541       vartmp(:)=SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day
1542       CALL histwrite_p (hist_id_stomate_IPCC, "fFire", itime, &
1543            vartmp, npts, hori_index)
1544       vartmp(:)=harvest_above/1e3/one_day
1545       CALL histwrite_p (hist_id_stomate_IPCC, "fHarvest", itime, &
1546            vartmp, npts, hori_index)
1547       vartmp(:)=cflux_prod_total/1e3/one_day
1548       CALL histwrite_p (hist_id_stomate_IPCC, "fLuc", itime, &
1549            vartmp, npts, hori_index)
1550       vartmp(:)=cflux_prod_harvest_total/1e3/one_day
1551       CALL histwrite_p (hist_id_stomate_IPCC, "fWoodharvest", itime, &
1552            vartmp, npts, hori_index)
1553       ! co2_to_bm is not added as it is already included in gpp
1554       vartmp(:)=(SUM((gpp_daily+co2_to_bm-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
1555            &        *veget_cov_max,dim=2)-cflux_prod_total-cflux_prod_harvest_total-harvest_above)/1e3/one_day
1556       CALL histwrite_p (hist_id_stomate_IPCC, "nbp", itime, &
1557            vartmp, npts, hori_index)
1558       vartmp(:)=SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*veget_cov_max,dim=2)/1e3/one_day
1559       CALL histwrite_p (hist_id_stomate_IPCC, "fVegLitter", itime, &
1560            vartmp, npts, hori_index)
1561       vartmp(:)=SUM(SUM(soilcarbon_input,dim=2)*veget_cov_max,dim=2)/1e3/one_day
1562       CALL histwrite_p (hist_id_stomate_IPCC, "fLitterSoil", itime, &
1563            vartmp, npts, hori_index)
1564       vartmp(:)=SUM(biomass(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3
1565       CALL histwrite_p (hist_id_stomate_IPCC, "cLeaf", itime, &
1566            vartmp, npts, hori_index)
1567       vartmp(:)=SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*veget_cov_max,dim=2)/1e3
1568       CALL histwrite_p (hist_id_stomate_IPCC, "cStem", itime, &
1569            vartmp, npts, hori_index)
1570       vartmp(:)=SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + biomass(:,:,iheartbelow,icarbon) ) &
1571            &        *veget_cov_max,dim=2)/1e3
1572       CALL histwrite_p (hist_id_stomate_IPCC, "cRoot", itime, &
1573            vartmp, npts, hori_index)
1574       vartmp(:)=SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*veget_cov_max,dim=2)/1e3
1575       CALL histwrite_p (hist_id_stomate_IPCC, "cMisc", itime, &
1576            vartmp, npts, hori_index)
1577       vartmp(:)=SUM((litter(:,istructural,:,iabove,icarbon)+litter(:,imetabolic,:,iabove,icarbon))*&
1578            veget_cov_max,dim=2)/1e3
1579       CALL histwrite_p (hist_id_stomate_IPCC, "cLitterAbove", itime, &
1580            vartmp, npts, hori_index)
1581       vartmp(:)=SUM((litter(:,istructural,:,ibelow,icarbon)+litter(:,imetabolic,:,ibelow,icarbon))*&
1582            veget_cov_max,dim=2)/1e3
1583       CALL histwrite_p (hist_id_stomate_IPCC, "cLitterBelow", itime, &
1584            vartmp, npts, hori_index)
1585       vartmp(:)=SUM(carbon(:,iactive,:)*veget_cov_max,dim=2)/1e3
1586       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilFast", itime, &
1587            vartmp, npts, hori_index)
1588       vartmp(:)=SUM(carbon(:,islow,:)*veget_cov_max,dim=2)/1e3
1589       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilMedium", itime, &
1590            vartmp, npts, hori_index)
1591       vartmp(:)=SUM(carbon(:,ipassive,:)*veget_cov_max,dim=2)/1e3
1592       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilSlow", itime, &
1593            vartmp, npts, hori_index)
1594       DO j=1,nvm
1595          histvar(:,j)=veget_cov_max(:,j)*100
1596       ENDDO
1597       CALL histwrite_p (hist_id_stomate_IPCC, "landCoverFrac", itime, &
1598            histvar, npts*nvm, horipft_index)
1599       !-
1600       vartmp(:)=zero
1601       DO j = 2,nvm
1602          IF (is_deciduous(j)) THEN
1603             vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1604          ENDIF
1605       ENDDO
1606       CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimDec", itime, &
1607            vartmp, npts, hori_index)
1608       !-
1609       vartmp(:)=zero
1610       DO j = 2,nvm
1611          IF (is_evergreen(j)) THEN
1612             vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1613          ENDIF
1614       ENDDO
1615       CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimEver", itime, &
1616            vartmp, npts, hori_index)
1617       !-
1618       vartmp(:)=zero
1619       DO j = 2,nvm
1620          IF ( .NOT.(is_c4(j)) ) THEN
1621             vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1622          ENDIF
1623       ENDDO
1624       CALL histwrite_p (hist_id_stomate_IPCC, "c3PftFrac", itime, &
1625            vartmp, npts, hori_index)
1626       !-
1627       vartmp(:)=zero
1628       DO j = 2,nvm
1629          IF ( is_c4(j) ) THEN
1630             vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1631          ENDIF
1632       ENDDO
1633       CALL histwrite_p (hist_id_stomate_IPCC, "c4PftFrac", itime, &
1634            vartmp, npts, hori_index)
1635       !-
1636       vartmp(:)=SUM(resp_growth*veget_cov_max,dim=2)/1e3/one_day
1637       CALL histwrite_p (hist_id_stomate_IPCC, "rGrowth", itime, &
1638            vartmp, npts, hori_index)
1639       vartmp(:)=SUM(resp_maint*veget_cov_max,dim=2)/1e3/one_day
1640       CALL histwrite_p (hist_id_stomate_IPCC, "rMaint", itime, &
1641            vartmp, npts, hori_index)
1642       vartmp(:)=SUM(bm_alloc(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3/one_day
1643       CALL histwrite_p (hist_id_stomate_IPCC, "nppLeaf", itime, &
1644            vartmp, npts, hori_index)
1645       vartmp(:)=SUM((bm_alloc(:,:,isapabove,icarbon) + bm_alloc(:,:,iheartabove,icarbon))*veget_cov_max,dim=2)/1e3/one_day
1646       CALL histwrite_p (hist_id_stomate_IPCC, "nppStem", itime, &
1647            vartmp, npts, hori_index)
1648       vartmp(:)=SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )*veget_cov_max,dim=2)/1e3/one_day
1649       CALL histwrite_p (hist_id_stomate_IPCC, "nppRoot", itime, &
1650            vartmp, npts, hori_index)
1651
1652       CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_X', itime, &
1653            resolution(:,1), npts, hori_index)
1654       CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_Y', itime, &
1655            resolution(:,2), npts, hori_index)
1656       CALL histwrite_p (hist_id_stomate_IPCC, 'CONTFRAC', itime, &
1657            contfrac(:), npts, hori_index)
1658
1659    ENDIF
1660
1661    IF (printlev>=4) WRITE(numout,*) 'Leaving stomate_lpj'
1662
1663  END SUBROUTINE StomateLpj
1664
1665
1666!! ================================================================================================================================
1667!! SUBROUTINE   : harvest
1668!!
1669!>\BRIEF        Harvest of croplands
1670!!
1671!! DESCRIPTION  : To take into account biomass harvest from crop (mainly to take
1672!! into account for the reduced litter input and then decreased soil carbon. it is a
1673!! constant (40\%) fraction of above ground biomass.
1674!!
1675!! RECENT CHANGE(S) : None
1676!!
1677!! MAIN OUTPUT VARIABLE(S): ::harvest_above the harvested biomass
1678!!
1679!! REFERENCE(S) :
1680!! - Piao, S., P. Ciais, P. Friedlingstein, N. de Noblet-Ducoudre, P. Cadule, N. Viovy, and T. Wang. 2009.
1681!!   Spatiotemporal patterns of terrestrial carbon cycle during the 20th century. Global Biogeochemical
1682!!   Cycles 23:doi:10.1029/2008GB003339.
1683!!
1684!! FLOWCHART    : None
1685!! \n
1686!_ ================================================================================================================================
1687
1688  SUBROUTINE harvest(npts, dt_days, veget_cov_max, &
1689       bm_to_litter, turnover_daily, &
1690       harvest_above)
1691
1692  !! 0. Variable and parameter declaration
1693
1694    !! 0.1 Input variables
1695
1696    INTEGER, INTENT(in)                                    :: npts            !! Domain size (unitless)
1697    REAL(r_std), INTENT(in)                                :: dt_days         !! Time step (days)                               
1698    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: veget_cov_max       !! new "maximal" coverage fraction of a PFT (LAI ->
1699                                                                              !! infinity) on ground @tex $(m^2 m^{-2})$ @endtex
1700   
1701   !! 0.2 Output variables
1702   
1703   !! 0.3 Modified variables
1704
1705    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! [DISPENSABLE] conversion of biomass to litter
1706                                                                                     !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1707    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily   !! Turnover rates
1708                                                                               !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1709    REAL(r_std), DIMENSION(npts), INTENT(inout)            :: harvest_above    !! harvest above ground biomass for agriculture
1710                                                                               !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1711    !! 0.4 Local variables
1712
1713    INTEGER(i_std)                                         :: i, j, k, l, m    !! indices                       
1714    REAL(r_std)                                            :: above_old        !! biomass of previous time step
1715                                                                               !! @tex $(gC m^{-2})$ @endtex
1716!_ ================================================================================================================================
1717
1718  !! 1. Yearly initialisation
1719
1720    above_old             = zero
1721    harvest_above         = zero
1722
1723    DO i = 1, npts
1724       DO j = 1,nvm
1725          IF (.NOT. natural(j)) THEN
1726             above_old = turnover_daily(i,j,ileaf,icarbon) + turnover_daily(i,j,isapabove,icarbon) + &
1727                  &       turnover_daily(i,j,iheartabove,icarbon) + turnover_daily(i,j,ifruit,icarbon) + &
1728                  &       turnover_daily(i,j,icarbres,icarbon) + turnover_daily(i,j,isapbelow,icarbon) + &
1729                  &       turnover_daily(i,j,iheartbelow,icarbon) + turnover_daily(i,j,iroot,icarbon)
1730
1731             turnover_daily(i,j,ileaf,icarbon) = turnover_daily(i,j,ileaf,icarbon)*frac_turnover_daily
1732             turnover_daily(i,j,isapabove,icarbon) = turnover_daily(i,j,isapabove,icarbon)*frac_turnover_daily
1733             turnover_daily(i,j,isapbelow,icarbon) = turnover_daily(i,j,isapbelow,icarbon)*frac_turnover_daily
1734             turnover_daily(i,j,iheartabove,icarbon) = turnover_daily(i,j,iheartabove,icarbon)*frac_turnover_daily
1735             turnover_daily(i,j,iheartbelow,icarbon) = turnover_daily(i,j,iheartbelow,icarbon)*frac_turnover_daily
1736             turnover_daily(i,j,iroot,icarbon) = turnover_daily(i,j,iroot,icarbon)*frac_turnover_daily
1737             turnover_daily(i,j,ifruit,icarbon) = turnover_daily(i,j,ifruit,icarbon)*frac_turnover_daily
1738             turnover_daily(i,j,icarbres,icarbon) = turnover_daily(i,j,icarbres,icarbon)*frac_turnover_daily
1739             harvest_above(i)  = harvest_above(i) + veget_cov_max(i,j) * above_old *(un - frac_turnover_daily)
1740          ENDIF
1741       ENDDO
1742    ENDDO
1743
1744!!$    harvest_above = harvest_above
1745  END SUBROUTINE harvest
1746END MODULE stomate_lpj
Note: See TracBrowser for help on using the repository browser.