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

Last change on this file since 6151 was 6151, checked in by fabienne.maignan, 5 years ago

Adding outputs for respiration of biomass compartments

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