source: branches/publications/ORCHIDEE-LEAK-r5919/src_stomate/stomate_lpj.f90 @ 5925

Last change on this file since 5925 was 5919, checked in by ronny.lauerwald, 5 years ago

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