source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_stomate/stomate_lpj.f90 @ 8398

Last change on this file since 8398 was 4189, checked in by nicolas.vuichard, 7 years ago

Fixes the bug reported in ticket #3272 as proposed by Albert & Fabienne
  • Property svn:keywords set to HeadURL Date Author Revision
File size: 77.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!! SVN          :
19!! $HeadURL$
20!! $Date$
21!! $Revision$
22!! \n
23!_ ================================================================================================================================
24
25MODULE stomate_lpj
26
27  ! modules used:
28
29  USE ioipsl_para
30  USE xios_orchidee
31  USE grid
32  USE stomate_data
33  USE constantes
34  USE constantes_soil
35  USE pft_parameters
36  USE lpj_constraints
37  USE lpj_pftinout
38  USE lpj_kill
39  USE lpj_crown
40  USE lpj_fire
41  USE lpj_gap
42  USE lpj_light
43  USE lpj_establish
44  USE lpj_cover
45  USE stomate_prescribe
46  USE stomate_phenology
47  USE stomate_alloc
48  USE stomate_npp
49  USE stomate_turnover
50  USE stomate_litter
51  USE stomate_soilcarbon
52  USE stomate_vmax
53  USE stomate_lcchange
54!  USE 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, dead_leaves, carbon, lignin_struc, &
151       veget_max, veget_max_new, npp_longterm, lm_lastyearmax, veget_lastlight, &
152       everywhere, need_adjacent, RIP_time, &
153       lai, rprof,npp_daily, turnover_daily, turnover_time,&
154       control_moist, control_temp, soilcarbon_input, &
155       co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, &
156       height, deadleaf_cover, vcmax, &
157       bm_to_litter, &
158       prod10,prod100,flux10, flux100, &
159       convflux,cflux_prod10,cflux_prod100, harvest_above, carb_mass_total, &
160       fpc_max, MatrixA, &
161       Tseason, Tmin_spring_time, begin_leaves, onset_date)
162   
163  !! 0. Variable and parameter declaration
164
165    !! 0.1 input
166
167    INTEGER(i_std), INTENT(in)                                 :: npts                 !! Domain size (unitless)
168    REAL(r_std), INTENT(in)                                    :: dt_days              !! Time step of Stomate (days)
169    INTEGER(i_std), DIMENSION(npts,NbNeighb), INTENT(in)       :: neighbours           !! Indices of the 8 neighbours of each grid
170                                                                                       !! point [1=North and then clockwise]
171    REAL(r_std), DIMENSION(npts,2), INTENT(in)                 :: resolution           !! Resolution at each grid point (m) 
172                                                                                       !! [1=E-W, 2=N-S]
173    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: clay                 !! Clay fraction (0 to 1, unitless)
174    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: herbivores           !! Time constant of probability of a leaf to
175                                                                                       !! be eaten by a herbivore (days)
176    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: tsurf_daily          !! Daily surface temperatures (K)
177    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: tsoil_daily          !! Daily soil temperatures (K)
178    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_daily            !! Daily 2 meter temperatures (K)
179    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_min_daily        !! Daily minimum 2 meter temperatures (K)
180    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: litterhum_daily      !! Daily litter humidity (0 to 1, unitless)
181    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: soilhum_daily        !! Daily soil humidity (0 to 1, unitless)
182    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: maxmoiavail_lastyear !! Last year's maximum moisture availability
183                                                                                       !! (0 to 1, unitless)
184    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: minmoiavail_lastyear !! Last year's minimum moisture availability
185                                                                                       !! (0 to 1, unitless)
186    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: gdd0_lastyear        !! Last year's GDD0 (K)
187    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: precip_lastyear      !! Lastyear's precipitation
188                                                                                       !! @tex $(mm year^{-1})$ @endtex
189                                                                                       !! to determine if establishment possible
190    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: moiavail_month       !! "Monthly" moisture availability (0 to 1,
191                                                                                       !! unitless)
192    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: moiavail_week        !! "Weekly" moisture availability
193                                                                                       !! (0 to 1, unitless)
194    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_longterm         !! "Long term" 2 meter reference
195                                                                                       !! temperatures (K)
196    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_month            !! "Monthly" 2-meter temperatures (K)
197    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: t2m_week             !! "Weekly" 2-meter temperatures (K)
198    REAL(r_std), DIMENSION(npts), INTENT(in)                   :: Tseason              !! "seasonal" 2-meter temperatures (K)
199    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: Tmin_spring_time     !! Number of days after begin_leaves (leaf onset)
200    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: onset_date           !! Date in the year at when the leaves started to grow(begin_leaves)
201    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: tsoil_month          !! "Monthly" soil temperatures (K)
202    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: soilhum_month        !! "Monthly" soil humidity
203                                                                                       !! (0 to 1, unitless)
204    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gdd_m5_dormance      !! Growing degree days (K), threshold -5 deg
205                                                                                       !! C (for phenology)
206    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gdd_from_growthinit  !! growing degree days, since growthinit for crops
207    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: gdd_midwinter        !! Growing degree days (K), since midwinter
208                                                                                       !! (for phenology) - this is written to the history files
209    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: ncd_dormance         !! Number of chilling days (days), since
210                                                                                       !! leaves were lost (for phenology)
211    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: ngd_minus5           !! Number of growing days (days), threshold
212                                                                                       !! -5 deg C (for phenology)
213    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)        :: turnover_longterm    !! "Long term" turnover rate 
214                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
215    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: gpp_daily            !! Daily gross primary productivity 
216                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
217    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: time_hum_min         !! Time elapsed since strongest moisture
218                                                                                       !! availability (days)
219    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: maxfpc_lastyear      !! Last year's maximum foliage projected
220                                                                                       !! coverage for each natural PFT,
221                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
222    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)        :: resp_maint_part      !! Maintenance respiration of different
223                                                                                       !! plant parts 
224                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
225    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: fpc_max              !! "Maximal" coverage fraction of a PFT (LAI
226                                                                                       !! -> infinity) on ground 
227                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
228    REAL(r_std), DIMENSION(npts,nvm),INTENT(in)                :: veget_max_new        !! Old timestep "maximal" coverage fraction of a PFT
229
230  !! 0.2 Output variables
231   
232    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: npp_daily            !! Net primary productivity
233                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
234    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out) :: turnover_daily       !! Turnover rates
235                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
236    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: co2_to_bm            !! CO2 taken up from atmosphere when
237                                                                                       !! introducing a new PFT (introduced for
238                                                                                       !! carbon balance closure)
239                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
240    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: co2_fire             !! Carbon emitted into the atmosphere by
241                                                                                       !! fire (living and dead biomass) 
242                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
243    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: resp_hetero          !! Heterotrophic respiration
244                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
245    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: resp_maint           !! Maintenance respiration 
246                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
247    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: resp_growth          !! Growth respiration 
248                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
249   
250    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: deadleaf_cover       !! Fraction of soil covered by dead leaves
251                                                                                       !! (0 to 1, unitless)
252    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: vcmax                !! Maximum rate of carboxylation
253    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out):: bm_to_litter      !! Conversion of biomass to litter
254                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
255    LOGICAL, DIMENSION(npts,nvm), INTENT(out)                  :: begin_leaves         !! signal to start putting leaves on (true/false)
256
257    !! 0.3 Modified variables
258   
259    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: height               !! Height of vegetation (m)
260    REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout)          :: control_moist        !! Moisture control of heterotrophic
261                                                                                       !! respiration (0 to 1, unitless)
262    REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout)          :: control_temp         !! Temperature control of heterotrophic
263                                                                                       !! respiration, above and below
264                                                                                       !! (0 to 1, unitless)
265    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout)      :: soilcarbon_input     !! Quantity of carbon going into carbon
266                                                                                       !! pools from litter decomposition 
267                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
268    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: lai                  !! Leaf area index OF AN INDIVIDUAL PLANT,
269                                                                                       !! where a PFT contains n indentical plants
270                                                                                       !! i.e., using the mean individual approach
271                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
272    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: rprof                !! Prescribed root depth (m)
273    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: PFTpresent           !! Tab indicating which PFTs are present in
274                                                                                       !! each pixel
275    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: age                  !! Age (years)   
276    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: fireindex            !! Probability of fire (0 to 1, unitless)
277    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: firelitter           !! Longer term litter above the ground that
278                                                                                       !! can be burned, @tex $(gC m^{-2})$ @endtex
279    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_age             !! Leaf age (days)
280    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_frac            !! Fraction of leaves in leaf age class,
281                                                                                       !! (0 to 1, unitless)
282    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass        !! Biomass @tex $(gC m^{-2})$ @endtex
283    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: ind                  !! Density of individuals
284                                                                                       !! @tex $(m^{-2})$ @endtex
285    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: adapted              !! Adaptation of PFT (killed if too cold)
286                                                                                       !! (0 to 1, unitless)
287    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: regenerate           !! "Fitness": Winter sufficiently cold for
288                                                                                       !! PFT regeneration ? (0 to 1, unitless)
289    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: senescence           !! Flag for setting senescence stage (only
290                                                                                       !! for deciduous trees)
291    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: when_growthinit      !! How many days ago was the beginning of
292                                                                                       !! the growing season (days)
293    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)      :: litterpart           !! Fraction of litter above the ground
294                                                                                       !! belonging to different PFTs
295                                                                                       !! (0 to 1, unitless)
296    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout):: litter     !! Metabolic and structural litter, above
297                                                                                       !! and below ground
298                                                                                       !! @tex $(gC m^{-2})$ @endtex
299    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)      :: dead_leaves          !! Dead leaves on ground, per PFT, metabolic
300                                                                                       !! and structural, 
301                                                                                       !! @tex $(gC m^{-2})$ @endtex
302    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout)      :: carbon               !! Carbon pool: active, slow, or passive,
303                                                                                       !! @tex $(gC m^{-2})$ @endtex 
304    REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout)      :: lignin_struc         !! Ratio of Lignin/Carbon in structural
305                                                                                       !! litter, above and below ground, 
306                                                                                       !! @tex $(gC m^{-2})$ @endtex
307    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: veget_max            !! "Maximal" coverage fraction of a PFT (LAI
308                                                                                       !! -> infinity) on ground
309    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: npp_longterm         !! "Long term" mean yearly primary
310                                                                                       !! productivity
311                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
312    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: lm_lastyearmax       !! Last year's maximum leaf mass, for each
313                                                                                       !! PFT @tex $(gC m^{-2})$ @endtex
314    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: veget_lastlight      !! Vegetation fractions (on ground) after
315                                                                                       !! last light competition 
316                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
317    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: everywhere           !! Is the PFT everywhere in the grid box or
318                                                                                       !! very localized (after its introduction)
319                                                                                       !! (unitless)
320    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                :: need_adjacent        !! In order for this PFT to be introduced,
321                                                                                       !! does it have to be present in an
322                                                                                       !! adjacent grid box?
323    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: RIP_time             !! How much time ago was the PFT eliminated
324                                                                                       !! for the last time (y)
325    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: turnover_time        !! Turnover_time of leaves for grasses
326                                                                                       !! (days)
327    REAL(r_std),DIMENSION(npts,0:10), INTENT(inout)            :: prod10               !! Products remaining in the 10
328                                                                                       !! year-turnover pool after the annual
329                                                                                       !! release for each compartment (10
330                                                                                       !! + 1 : input from year of land cover
331                                                                                       !! change) @tex $(gC m^{-2})$ @endtex
332    REAL(r_std),DIMENSION(npts,0:100), INTENT(inout)           :: prod100              !! Products remaining in the 100
333                                                                                       !! year-turnover pool after the annual
334                                                                                       !! release for each compartment (100
335                                                                                       !! + 1 : input from year of land cover
336                                                                                       !! change) @tex $(gC m^{-2})$ @endtex
337    REAL(r_std),DIMENSION(npts,10), INTENT(inout)              :: flux10               !! Annual release from the 10
338                                                                                       !! year-turnover pool compartments 
339                                                                                       !! @tex $(gC m^{-2})$ @endtex
340    REAL(r_std),DIMENSION(npts,100), INTENT(inout)             :: flux100              !! Annual release from the 100
341                                                                                       !! year-turnover pool compartments 
342                                                                                       !! @tex $(gC m^{-2})$ @endtex
343    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: convflux             !! Release during first year following land
344                                                                                       !! cover change @tex $(gC m^{-2})$ @endtex
345    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: cflux_prod10         !! Total annual release from the 10
346                                                                                       !! year-turnover pool
347                                                                                       !! @tex $(gC m^{-2})$ @endtex
348    REAL(r_std),DIMENSION(npts), INTENT(inout)                 :: cflux_prod100        !! Total annual release from the 100
349                                                                                       !! year-turnover pool
350                                                                                       !! @tex $(gC m^{-2})$ @endtex
351    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: harvest_above        !! Harvest above ground biomass for
352                                                                                       !! agriculture @tex $(gC m^{-2})$ @endtex
353    REAL(r_std), DIMENSION(npts), INTENT(inout)                :: carb_mass_total      !! Carbon Mass total (soil, litter, veg)
354                                                                                       !! @tex $(gC m^{-2})$ @endtex 
355    REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA         !! Matrix containing the fluxes 
356                                                                                       !! between the carbon pools
357                                                                                       !! per sechiba time step
358                                                                                       !! @tex $(gC.m^2.day^{-1})$ @endtex
359
360    !! 0.4 Local variables
361
362    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_bm_to_litter    !! Total conversion of biomass to litter
363                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
364    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_live_biomass    !! Total living biomass 
365                                                                                       !! @tex $(gC m{-2})$ @endtex
366    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements)           :: bm_alloc            !! Biomass increase, i.e. NPP per plant part
367                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
368    REAL(r_std), DIMENSION(npts,nvm,nelements)                  :: tot_turnover        !! Total turnover rate 
369                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
370    REAL(r_std), DIMENSION(npts,nvm)                            :: tot_litter_soil_carb!! Total soil and litter carbon 
371                                                                                       !! @tex $(gC m^{-2})$ @endtex
372    REAL(r_std), DIMENSION(npts,nvm)                            :: tot_litter_carb     !! Total litter carbon
373                                                                                       !! @tex $(gC m^{-2})$ @endtex
374    REAL(r_std), DIMENSION(npts,nvm)                            :: tot_soil_carb       !! Total soil carbon 
375                                                                                       !! @tex $(gC m^{-2})$ @endtex
376    REAL(r_std), DIMENSION(npts)                                :: carb_mass_variation !! Carbon Mass variation 
377                                                                                       !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
378    REAL(r_std), DIMENSION(npts,nvm)                            :: cn_ind              !! Crown area of individuals
379                                                                                       !! @tex $(m^{2})$ @endtex
380    REAL(r_std), DIMENSION(npts,nvm)                            :: woodmass_ind        !! Woodmass of individuals (gC)
381    REAL(r_std), DIMENSION(npts,nvm,nparts)                     :: f_alloc             !! Fraction that goes into plant part
382                                                                                       !! (0 to 1, unitless)
383    REAL(r_std), DIMENSION(npts)                                :: avail_tree          !! Space availability for trees
384                                                                                       !! (0 to 1, unitless)
385    REAL(r_std), DIMENSION(npts)                                :: avail_grass         !! Space availability for grasses
386                                                                                       !! (0 to 1, unitless)
387    INTEGER                                                     :: j,k
388    REAL(r_std),DIMENSION(npts)                                 :: prod10_total        !! Total products remaining in the pool
389                                                                                       !! after the annual release
390                                                                                       !! @tex $(gC m^{-2})$ @endtex
391    REAL(r_std),DIMENSION(npts)                                 :: prod100_total       !! Total products remaining in the pool
392                                                                                       !! after the annual release
393                                                                                       !! @tex $(gC m^{-2})$ @endtex
394    REAL(r_std),DIMENSION(npts)                                 :: cflux_prod_total    !! Total flux from conflux and the 10/100
395                                                                                       !! year-turnover pool
396                                                                                       !! @tex $(gC m^{-2} year^{-1})$ @endtex
397    REAL(r_std),DIMENSION(npts,nvm)                             :: veget_max_tmp       !! "Maximal" coverage fraction of a PFT 
398                                                                                       !! (LAI-> infinity) on ground (unitless)
399    REAL(r_std), DIMENSION(npts,nvm)                            :: mortality           !! Fraction of individual dying this time
400                                                                                       !! step (0 to 1, unitless)
401    REAL(r_std), DIMENSION(npts)                                :: vartmp              !! Temporary variable used to add history
402    REAL(r_std), DIMENSION(npts,nvm)                            :: histvar             !! History variables
403!_ ================================================================================================================================
404
405    IF (printlev>=3) WRITE(numout,*) 'Entering stomate_lpj'
406
407 
408  !! 1. Initializations
409   
410    !! 1.1 Initialize variables to zero
411    co2_to_bm(:,:) = zero
412    co2_fire(:,:) = zero
413    npp_daily(:,:) = zero
414    resp_maint(:,:) = zero
415    resp_growth(:,:) = zero
416    harvest_above(:) = zero
417    bm_to_litter(:,:,:,:) = zero
418    cn_ind(:,:) = zero
419    woodmass_ind(:,:) = zero
420    turnover_daily(:,:,:,:) = zero
421   
422    !! 1.2  Initialize variables to veget_max
423    veget_max_tmp(:,:) = veget_max(:,:)
424
425    !! 1.3 Calculate some vegetation characteristics
426   
427    !! 1.3.1 Calculate some vegetation characteristics
428    !        Calculate cn_ind (individual crown mass) and individual height from
429    !        state variables if running DGVM or dynamic mortality in static cover mode
430    !??        Explain (maybe in the header once) why you mulitply with veget_max in the DGVM
431    !??        and why you don't multiply with veget_max in stomate.
432    IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
433       IF(ok_dgvm) THEN
434          WHERE (ind(:,:).GT.min_stomate)
435             woodmass_ind(:,:) = &
436                  ((biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
437                  +biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)) & 
438                  *veget_max(:,:))/ind(:,:)
439          ENDWHERE
440       ELSE
441          WHERE (ind(:,:).GT.min_stomate)
442             woodmass_ind(:,:) = &
443                  (biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
444                  +biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon))/ind(:,:)
445          ENDWHERE
446       ENDIF
447
448       CALL crown (npts,  PFTpresent, &
449            ind, biomass, woodmass_ind, &
450            veget_max, cn_ind, height)
451    ENDIF
452
453    !! 1.3.2 Prescribe characteristics if the vegetation is not dynamic
454    !        IF the DGVM is not activated, the density of individuals and their crown
455    !        areas don't matter, but they should be defined for the case we switch on
456    !        the DGVM afterwards. At the first call, if the DGVM is not activated,
457    !        impose a minimum biomass for prescribed PFTs and declare them present.
458    CALL prescribe (npts, &
459         veget_max, dt_days, PFTpresent, everywhere, when_growthinit, &
460         biomass, leaf_frac, ind, cn_ind, co2_to_bm)
461
462
463  !! 2. Climatic constraints for PFT presence and regenerativeness
464
465    !   Call this even when DGVM is not activated so that "adapted" and "regenerate"
466    !   are kept up to date for the moment when the DGVM is activated.
467    CALL constraints (npts, dt_days, &
468         t2m_month, t2m_min_daily,when_growthinit, Tseason, &
469         adapted, regenerate)
470
471   
472  !! 3. Determine introduction and elimination of PTS based on climate criteria
473 
474    IF ( ok_dgvm ) THEN
475     
476       !! 3.1 Calculate introduction and elimination
477       CALL pftinout (npts, dt_days, adapted, regenerate, &
478            neighbours, veget_max, &
479            biomass, ind, cn_ind, age, leaf_frac, npp_longterm, lm_lastyearmax, senescence, &
480            PFTpresent, everywhere, when_growthinit, need_adjacent, RIP_time, &
481            co2_to_bm, &
482            avail_tree, avail_grass)
483
484       !! 3.2 Reset attributes for eliminated PFTs.
485       !     This also kills PFTs that had 0 leafmass during the last year. The message
486       !     "... after pftinout" is misleading in this case.
487       CALL kill (npts, 'pftinout  ', lm_lastyearmax, &
488            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
489            lai, age, leaf_age, leaf_frac, npp_longterm, &
490            when_growthinit, everywhere, veget_max, bm_to_litter)
491
492       
493       !! 3.3 Calculate new crown area and diameter
494       !      Calculate new crown area, diameter and maximum vegetation cover**[No longer used in the subroutine]
495       !      unsure whether this is really required
496       !      - in theory this could ONLY be done at the END of stomate_lpj
497       !      calculate woodmass of individual tree
498       WHERE ((ind(:,:).GT.min_stomate))
499          WHERE  ( veget_max(:,:) .GT. min_stomate)
500             woodmass_ind(:,:) = &
501                  ((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
502                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))*veget_max(:,:))/ind(:,:)
503          ELSEWHERE
504             woodmass_ind(:,:) =(biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
505                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
506          ENDWHERE
507
508       ENDWHERE
509       
510       ! Calculate crown area and diameter for all PFTs (including the newly established)
511       CALL crown (npts, PFTpresent, &
512            ind, biomass, woodmass_ind, &
513            veget_max, cn_ind, height)
514
515    ENDIF
516   
517  !! 4. Phenology
518
519    !! 4.1 Write values to history file
520    !      Current values for ::when_growthinit
521    CALL xios_orchidee_send_field("WHEN_GROWTHINIT",when_growthinit)
522
523    CALL histwrite_p (hist_id_stomate, 'WHEN_GROWTHINIT', itime, when_growthinit, npts*nvm, horipft_index)
524
525    ! Set and write values for ::PFTpresent
526    WHERE(PFTpresent)
527       histvar=un
528    ELSEWHERE
529       histvar=zero
530    ENDWHERE
531
532    CALL xios_orchidee_send_field("PFTPRESENT",histvar)
533
534    CALL histwrite_p (hist_id_stomate, 'PFTPRESENT', itime, histvar, npts*nvm, horipft_index)
535
536    ! Set and write values for gdd_midwinter
537    WHERE(gdd_midwinter.EQ.undef)
538       histvar=val_exp
539    ELSEWHERE
540       histvar=gdd_midwinter
541    ENDWHERE
542
543    CALL xios_orchidee_send_field("GDD_MIDWINTER",histvar)
544
545    CALL histwrite_p (hist_id_stomate, 'GDD_MIDWINTER', itime, histvar, npts*nvm, horipft_index)
546
547    ! Set and write values for gdd_m5_dormance
548    WHERE(gdd_m5_dormance.EQ.undef)
549       histvar=val_exp
550    ELSEWHERE
551       histvar=gdd_m5_dormance
552    ENDWHERE
553   
554    CALL xios_orchidee_send_field('GDD_M5_DORMANCE',histvar)
555    CALL histwrite_p (hist_id_stomate, 'GDD_M5_DORMANCE', itime, histvar, npts*nvm, horipft_index)
556
557    ! Set and write values for ncd_dormance
558    WHERE(ncd_dormance.EQ.undef)
559       histvar=val_exp
560    ELSEWHERE
561       histvar=ncd_dormance
562    ENDWHERE
563
564    CALL xios_orchidee_send_field("NCD_DORMANCE",histvar)
565
566    CALL histwrite_p (hist_id_stomate, 'NCD_DORMANCE', itime, histvar, npts*nvm, horipft_index)
567
568    !! 4.2 Calculate phenology
569    CALL phenology (npts, dt_days, PFTpresent, &
570         veget_max, &
571         t2m_longterm, t2m_month, t2m_week, gpp_daily, &
572         maxmoiavail_lastyear, minmoiavail_lastyear, &
573         moiavail_month, moiavail_week, &
574         gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
575         senescence, time_hum_min, &
576         biomass, leaf_frac, leaf_age, &
577         when_growthinit, co2_to_bm, &
578         begin_leaves)
579   
580  !! 5. Allocate C to different plant parts
581   
582    CALL alloc (npts, dt_days, &
583         lai, veget_max, senescence, when_growthinit, &
584         moiavail_week, tsoil_month, soilhum_month, &
585         biomass, age, leaf_age, leaf_frac, rprof, f_alloc)
586
587  !! 6. NPP, maintenance and growth respiration
588
589    !! 6.1 Calculate NPP and respiration terms
590    CALL npp_calc (npts, dt_days, &
591         PFTpresent, &
592         t2m_daily, tsoil_daily, lai, rprof, &
593         gpp_daily, f_alloc, bm_alloc, resp_maint_part,&
594         biomass, leaf_age, leaf_frac, age, &
595         resp_maint, resp_growth, npp_daily)
596
597    !! 6.2 Kill slow growing PFTs in DGVM or STOMATE with constant mortality
598    IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
599       CALL kill (npts, 'npp       ', lm_lastyearmax,  &
600            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
601            lai, age, leaf_age, leaf_frac, npp_longterm, &
602            when_growthinit, everywhere, veget_max, bm_to_litter)
603
604       !! 6.2.1 Update wood biomass     
605       !        For the DGVM
606       IF(ok_dgvm) THEN
607          WHERE (ind(:,:).GT.min_stomate)
608             woodmass_ind(:,:) = &
609                  ((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
610                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon)) & 
611                  *veget_max(:,:))/ind(:,:)
612          ENDWHERE
613
614       ! For all pixels with individuals
615       ELSE
616          WHERE (ind(:,:).GT.min_stomate)
617             woodmass_ind(:,:) = &
618                  (biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
619                  + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
620          ENDWHERE
621       ENDIF ! ok_dgvm
622
623       !! 6.2.2 New crown area and maximum vegetation cover after growth
624       CALL crown (npts, PFTpresent, &
625            ind, biomass, woodmass_ind,&
626            veget_max, cn_ind, height)
627
628    ENDIF ! ok_dgvm
629   
630  !! 7. fire
631
632    !! 7.1. Burn PFTs
633    CALL fire (npts, dt_days, litterpart, &
634         litterhum_daily, t2m_daily, lignin_struc, veget_max, &
635         fireindex, firelitter, biomass, ind, &
636         litter, dead_leaves, bm_to_litter, &
637         co2_fire, MatrixA)
638
639    !! 7.2 Kill PFTs in DGVM
640    IF ( ok_dgvm ) THEN
641
642       ! reset attributes for eliminated PFTs
643       CALL kill (npts, 'fire      ', lm_lastyearmax, &
644            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
645            lai, age, leaf_age, leaf_frac, npp_longterm, &
646            when_growthinit, everywhere, veget_max, bm_to_litter)
647
648    ENDIF ! ok_dgvm
649 
650  !! 8. Tree mortality
651
652    ! Does not depend on age, therefore does not change crown area.
653    CALL gap (npts, dt_days, &
654         npp_longterm, turnover_longterm, lm_lastyearmax, &
655         PFTpresent, t2m_min_daily, Tmin_spring_time, &
656         biomass, ind, bm_to_litter, mortality)
657
658
659    IF ( ok_dgvm ) THEN
660
661       ! reset attributes for eliminated PFTs
662       CALL kill (npts, 'gap       ', lm_lastyearmax, &
663            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
664            lai, age, leaf_age, leaf_frac, npp_longterm, &
665            when_growthinit, everywhere, veget_max, bm_to_litter)
666
667    ENDIF
668
669  !! 9. Leaf senescence, new lai and other turnover processes
670
671    CALL turn (npts, dt_days, PFTpresent, &
672         herbivores, &
673         maxmoiavail_lastyear, minmoiavail_lastyear, &
674         moiavail_week,  moiavail_month,t2m_longterm, t2m_month, t2m_week, veget_max, &
675         gdd_from_growthinit, leaf_age, leaf_frac, age, lai, biomass, &
676         turnover_daily, senescence,turnover_time)
677
678    !! 10. Light competition
679   
680    !! If not using constant mortality then kill with light competition
681!    IF ( ok_dgvm .OR. .NOT.(lpj_gap_const_mort) ) THEN
682    IF ( ok_dgvm ) THEN
683 
684       !! 10.1 Light competition
685       CALL light (npts, dt_days, &
686            veget_max, fpc_max, PFTpresent, cn_ind, lai, maxfpc_lastyear, &
687            lm_lastyearmax, ind, biomass, veget_lastlight, bm_to_litter, mortality)
688       
689       !! 10.2 Reset attributes for eliminated PFTs
690       CALL kill (npts, 'light     ', lm_lastyearmax, &
691            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
692            lai, age, leaf_age, leaf_frac, npp_longterm, &
693            when_growthinit, everywhere, veget_max, bm_to_litter)
694
695    ENDIF
696
697   
698  !! 11. Establishment of saplings
699   
700    IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort ) THEN
701
702       !! 11.1 Establish new plants
703       CALL establish (npts, dt_days, PFTpresent, regenerate, &
704            neighbours, resolution, need_adjacent, herbivores, &
705            precip_lastyear, gdd0_lastyear, lm_lastyearmax, &
706            cn_ind, lai, avail_tree, avail_grass, npp_longterm, &
707            leaf_age, leaf_frac, &
708            ind, biomass, age, everywhere, co2_to_bm, veget_max, woodmass_ind, &
709            mortality, bm_to_litter)
710
711       !! 11.2 Calculate new crown area (and maximum vegetation cover)
712       CALL crown (npts, PFTpresent, &
713            ind, biomass, woodmass_ind, &
714            veget_max, cn_ind, height)
715
716    ENDIF
717
718  !! 12. Calculate final LAI and vegetation cover
719   
720    CALL cover (npts, cn_ind, ind, biomass, &
721         veget_max, veget_max_tmp, lai, &
722         litter, carbon, turnover_daily, bm_to_litter, &
723         co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily)
724
725  !! 13. Update litter pools to account for harvest
726 
727    ! the whole litter stuff:
728    !    litter update, lignin content, PFT parts, litter decay,
729    !    litter heterotrophic respiration, dead leaf soil cover.
730    !    No vertical discretisation in the soil for litter decay.\n
731    ! added by shilong for harvest
732    IF(harvest_agri) THEN
733       CALL harvest(npts, dt_days, veget_max, &
734            bm_to_litter, turnover_daily, &
735            harvest_above)
736    ENDIF
737
738    !! 14. Land cover change if it is time to do so
739    !! The flag do_now_lcchange is set in slowproc_main at the same time as the vegetation is read from file.
740    !! The vegetation fractions are not updated yet and will be updated in the end of sechiba_main.
741    IF (do_now_stomate_lcchange) THEN
742       CALL lcchange_main (npts, dt_days, veget_max, veget_max_new, &
743            biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
744            co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,flux10,flux100, &
745            prod10,prod100,convflux,cflux_prod10,cflux_prod100,leaf_frac,&
746            npp_longterm, lm_lastyearmax, litter, carbon)
747       do_now_stomate_lcchange=.FALSE.
748
749       ! Set the flag done_stomate_lcchange to be used in the end of sechiba_main to update the fractions.
750       done_stomate_lcchange=.TRUE.
751    ENDIF
752
753
754    !! 15. Calculate vcmax
755
756    CALL vmax (npts, dt_days, &
757         leaf_age, leaf_frac, &
758         vcmax)
759
760    !MM déplacement pour initialisation correcte des grandeurs cumulées :
761    cflux_prod_total(:) = convflux(:) + cflux_prod10(:) + cflux_prod100(:)
762    prod10_total(:)=SUM(prod10,dim=2)
763    prod100_total(:)=SUM(prod100,dim=2)
764   
765  !! 16. Total heterotrophic respiration
766
767    tot_soil_carb(:,:) = zero
768    tot_litter_carb(:,:) = zero
769    DO j=2,nvm
770
771       tot_litter_carb(:,j) = tot_litter_carb(:,j) + (litter(:,istructural,j,iabove,icarbon) + &
772            &          litter(:,imetabolic,j,iabove,icarbon) + &
773            &          litter(:,istructural,j,ibelow,icarbon) + litter(:,imetabolic,j,ibelow,icarbon))
774
775       tot_soil_carb(:,j) = tot_soil_carb(:,j) + (carbon(:,iactive,j) + &
776            &          carbon(:,islow,j)+  carbon(:,ipassive,j))
777
778    ENDDO
779    tot_litter_soil_carb(:,:) = tot_litter_carb(:,:) + tot_soil_carb(:,:)
780
781!!$     DO k = 1, nelements ! Loop over # elements
782!!$        tot_live_biomass(:,:,k) = biomass(:,:,ileaf,k) + biomass(:,:,isapabove,k) + biomass(:,:,isapbelow,k) +&
783!!$             &                    biomass(:,:,iheartabove,k) + biomass(:,:,iheartbelow,k) + &
784!!$             &                    biomass(:,:,iroot,k) + biomass(:,:,ifruit,k) + biomass(:,:,icarbres,k)
785!!$    END DO ! Loop over # elements
786
787    tot_live_biomass(:,:,:) = biomass(:,:,ileaf,:) + biomass(:,:,isapabove,:) + biomass(:,:,isapbelow,:) +&
788             &                    biomass(:,:,iheartabove,:) + biomass(:,:,iheartbelow,:) + &
789             &                    biomass(:,:,iroot,:) + biomass(:,:,ifruit,:) + biomass(:,:,icarbres,:)
790
791
792    tot_turnover(:,:,:) = turnover_daily(:,:,ileaf,:) + turnover_daily(:,:,isapabove,:) + &
793         &         turnover_daily(:,:,isapbelow,:) + turnover_daily(:,:,iheartabove,:) + &
794         &         turnover_daily(:,:,iheartbelow,:) + turnover_daily(:,:,iroot,:) + &
795         &         turnover_daily(:,:,ifruit,:) + turnover_daily(:,:,icarbres,:)
796
797    tot_bm_to_litter(:,:,:) = bm_to_litter(:,:,ileaf,:) + bm_to_litter(:,:,isapabove,:) +&
798         &             bm_to_litter(:,:,isapbelow,:) + bm_to_litter(:,:,iheartbelow,:) +&
799         &             bm_to_litter(:,:,iheartabove,:) + bm_to_litter(:,:,iroot,:) + &
800         &             bm_to_litter(:,:,ifruit,:) + bm_to_litter(:,:,icarbres,:)
801
802    carb_mass_variation(:)=-carb_mass_total(:)
803    carb_mass_total(:)=SUM((tot_live_biomass(:,:,icarbon)+tot_litter_carb+tot_soil_carb)*veget_max,dim=2) + &
804         &                 (prod10_total + prod100_total)
805    carb_mass_variation(:)=carb_mass_total(:)+carb_mass_variation(:)
806   
807  !! 17. Write history
808
809    CALL xios_orchidee_send_field("RESOLUTION_X",resolution(:,1))
810    CALL xios_orchidee_send_field("RESOLUTION_Y",resolution(:,2))
811    CALL xios_orchidee_send_field("CONTFRAC_STOMATE",contfrac(:))
812    CALL xios_orchidee_send_field("T2M_MONTH",t2m_month)
813    CALL xios_orchidee_send_field("T2M_WEEK",t2m_week)
814    CALL xios_orchidee_send_field("TSEASON",Tseason)
815    CALL xios_orchidee_send_field("TMIN_SPRING_TIME", Tmin_spring_time)
816    CALL xios_orchidee_send_field("ONSET_DATE",onset_date)
817    CALL xios_orchidee_send_field("FPC_MAX",fpc_max)
818    CALL xios_orchidee_send_field("MAXFPC_LASTYEAR",maxfpc_lastyear)
819    CALL xios_orchidee_send_field("HET_RESP",resp_hetero(:,:))
820    CALL xios_orchidee_send_field("CO2_FIRE",co2_fire)
821    CALL xios_orchidee_send_field("CO2_TAKEN",co2_to_bm)
822    CALL xios_orchidee_send_field("LAI",lai)
823    CALL xios_orchidee_send_field("VEGET_MAX",veget_max)
824    CALL xios_orchidee_send_field("NPP_STOMATE",npp_daily)
825    CALL xios_orchidee_send_field("GPP",gpp_daily)
826    CALL xios_orchidee_send_field("IND",ind)
827    CALL xios_orchidee_send_field("CN_IND",cn_ind)
828    CALL xios_orchidee_send_field("WOODMASS_IND",woodmass_ind)
829    CALL xios_orchidee_send_field("TOTAL_M",tot_live_biomass)
830    CALL xios_orchidee_send_field("MOISTRESS",moiavail_week)
831    CALL xios_orchidee_send_field("LEAF_M",biomass(:,:,ileaf,icarbon))
832    CALL xios_orchidee_send_field("SAP_M_AB",biomass(:,:,isapabove,icarbon))
833    CALL xios_orchidee_send_field("SAP_M_BE",biomass(:,:,isapbelow,icarbon))
834    CALL xios_orchidee_send_field("HEART_M_AB",biomass(:,:,iheartabove,icarbon))
835    CALL xios_orchidee_send_field("HEART_M_BE",biomass(:,:,iheartbelow,icarbon))
836    CALL xios_orchidee_send_field("ROOT_M",biomass(:,:,iroot,icarbon))
837    CALL xios_orchidee_send_field("FRUIT_M",biomass(:,:,ifruit,icarbon))
838    CALL xios_orchidee_send_field("RESERVE_M",biomass(:,:,icarbres,icarbon))
839    CALL xios_orchidee_send_field("TOTAL_TURN",tot_turnover)
840    CALL xios_orchidee_send_field("LEAF_TURN",turnover_daily(:,:,ileaf,icarbon))
841    CALL xios_orchidee_send_field("MAINT_RESP",resp_maint)
842    CALL xios_orchidee_send_field("GROWTH_RESP",resp_growth)
843    CALL xios_orchidee_send_field("SAP_AB_TURN",turnover_daily(:,:,isapabove,icarbon))
844    CALL xios_orchidee_send_field("ROOT_TURN",turnover_daily(:,:,iroot,icarbon))
845    CALL xios_orchidee_send_field("FRUIT_TURN",turnover_daily(:,:,ifruit,icarbon))
846    CALL xios_orchidee_send_field("TOTAL_BM_LITTER",tot_bm_to_litter(:,:,icarbon))
847    CALL xios_orchidee_send_field("LEAF_BM_LITTER",bm_to_litter(:,:,ileaf,icarbon))
848    CALL xios_orchidee_send_field("SAP_AB_BM_LITTER",bm_to_litter(:,:,isapabove,icarbon))
849    CALL xios_orchidee_send_field("SAP_BE_BM_LITTER",bm_to_litter(:,:,isapbelow,icarbon))
850    CALL xios_orchidee_send_field("HEART_AB_BM_LITTER",bm_to_litter(:,:,iheartabove,icarbon))
851    CALL xios_orchidee_send_field("HEART_BE_BM_LITTER",bm_to_litter(:,:,iheartbelow,icarbon))
852    CALL xios_orchidee_send_field("ROOT_BM_LITTER",bm_to_litter(:,:,iroot,icarbon))
853    CALL xios_orchidee_send_field("FRUIT_BM_LITTER",bm_to_litter(:,:,ifruit,icarbon))
854    CALL xios_orchidee_send_field("RESERVE_BM_LITTER",bm_to_litter(:,:,icarbres,icarbon))
855    CALL xios_orchidee_send_field("LITTER_STR_AB",litter(:,istructural,:,iabove,icarbon))
856    CALL xios_orchidee_send_field("LITTER_MET_AB",litter(:,imetabolic,:,iabove,icarbon))
857    CALL xios_orchidee_send_field("LITTER_STR_BE",litter(:,istructural,:,ibelow,icarbon))
858    CALL xios_orchidee_send_field("LITTER_MET_BE",litter(:,imetabolic,:,ibelow,icarbon))
859    CALL xios_orchidee_send_field("DEADLEAF_COVER",deadleaf_cover)
860    CALL xios_orchidee_send_field("TOTAL_SOIL_CARB",tot_litter_soil_carb)
861    CALL xios_orchidee_send_field("CARBON_ACTIVE",carbon(:,iactive,:))
862    CALL xios_orchidee_send_field("CARBON_SLOW",carbon(:,islow,:))
863    CALL xios_orchidee_send_field("CARBON_PASSIVE",carbon(:,ipassive,:))
864    CALL xios_orchidee_send_field("LITTERHUM",litterhum_daily)
865    CALL xios_orchidee_send_field("TURNOVER_TIME",turnover_time)
866    CALL xios_orchidee_send_field("PROD10",prod10)
867    CALL xios_orchidee_send_field("FLUX10",flux10)
868    CALL xios_orchidee_send_field("PROD100",prod100)
869    CALL xios_orchidee_send_field("FLUX100",flux100)
870    CALL xios_orchidee_send_field("CONVFLUX",convflux)
871    CALL xios_orchidee_send_field("CFLUX_PROD10",cflux_prod10)
872    CALL xios_orchidee_send_field("CFLUX_PROD100",cflux_prod100)
873    CALL xios_orchidee_send_field("HARVEST_ABOVE",harvest_above)
874    CALL xios_orchidee_send_field("VCMAX",vcmax)
875    CALL xios_orchidee_send_field("AGE",age)
876    CALL xios_orchidee_send_field("HEIGHT",height)
877    CALL xios_orchidee_send_field("FIREINDEX",fireindex(:,:))
878! ipcc history
879     CALL xios_orchidee_send_field("cVeg",SUM(tot_live_biomass(:,:,icarbon)*veget_max,dim=2)/1e3*contfrac)
880     CALL xios_orchidee_send_field("cLitter",SUM(tot_litter_carb*veget_max,dim=2)/1e3*contfrac)
881     CALL xios_orchidee_send_field("cSoil",SUM(tot_soil_carb*veget_max,dim=2)/1e3*contfrac)
882     CALL xios_orchidee_send_field("cProduct",(prod10_total + prod100_total)/1e3)
883     CALL xios_orchidee_send_field("cMassVariation",carb_mass_variation/1e3/one_day*contfrac)
884     CALL xios_orchidee_send_field("lai_ipcc",SUM(lai*veget_max,dim=2)*contfrac)
885     CALL xios_orchidee_send_field("gpp_ipcc",SUM(gpp_daily*veget_max,dim=2)/1e3/one_day*contfrac)
886     CALL xios_orchidee_send_field("ra",SUM((resp_maint+resp_growth)*veget_max,dim=2)/1e3/one_day*contfrac)
887     CALL xios_orchidee_send_field("npp_ipcc",SUM(npp_daily*veget_max,dim=2)/1e3/one_day*contfrac)
888     CALL xios_orchidee_send_field("rh",SUM(resp_hetero*veget_max,dim=2)/1e3/one_day*contfrac)
889     CALL xios_orchidee_send_field("fFire",SUM(co2_fire*veget_max,dim=2)/1e3/one_day*contfrac)
890     CALL xios_orchidee_send_field("fHarvest",harvest_above/1e3/one_day*contfrac)
891     CALL xios_orchidee_send_field("fLuc",cflux_prod_total/1e3/one_day*contfrac)
892     CALL xios_orchidee_send_field("nbp",(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
893            &        *veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_day*contfrac)
894     CALL xios_orchidee_send_field("fVegLitter",SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*&
895          veget_max,dim=2)/1e3/one_day*contfrac)
896     CALL xios_orchidee_send_field("fLitterSoil",SUM(SUM(soilcarbon_input,dim=2)*veget_max,dim=2)/1e3/one_day*contfrac)
897     CALL xios_orchidee_send_field("cLeaf",SUM(biomass(:,:,ileaf,icarbon)*veget_max,dim=2)/1e3*contfrac)
898     CALL xios_orchidee_send_field("cWood",SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*&
899          veget_max,dim=2)/1e3*contfrac)
900     CALL xios_orchidee_send_field("cRoot",SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + &
901          biomass(:,:,iheartbelow,icarbon) )*veget_max,dim=2)/1e3*contfrac)
902     CALL xios_orchidee_send_field("cMisc",SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*&
903          veget_max,dim=2)/1e3*contfrac)
904     CALL xios_orchidee_send_field("cLitterAbove",SUM((litter(:,istructural,:,iabove,icarbon)+&
905          litter(:,imetabolic,:,iabove,icarbon))*veget_max,dim=2)/1e3*contfrac)
906     CALL xios_orchidee_send_field("cLitterBelow",SUM((litter(:,istructural,:,ibelow,icarbon)+&
907          litter(:,imetabolic,:,ibelow,icarbon))*veget_max,dim=2)/1e3*contfrac)
908     CALL xios_orchidee_send_field("cSoilFast",SUM(carbon(:,iactive,:)*veget_max,dim=2)/1e3*contfrac)
909     CALL xios_orchidee_send_field("cSoilMedium",SUM(carbon(:,islow,:)*veget_max,dim=2)/1e3*contfrac)
910     CALL xios_orchidee_send_field("cSoilSlow",SUM(carbon(:,ipassive,:)*veget_max,dim=2)/1e3*contfrac)
911       DO j=1,nvm
912          histvar(:,j)=veget_max(:,j)*contfrac(:)*100
913       ENDDO
914     CALL xios_orchidee_send_field("landCoverFrac",histvar)
915       vartmp(:)=zero
916       DO j = 2,nvm
917          IF (is_deciduous(j)) THEN
918             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
919          ENDIF
920       ENDDO
921     CALL xios_orchidee_send_field("treeFracPrimDec",vartmp)
922       vartmp(:)=zero
923       DO j = 2,nvm
924          IF (is_evergreen(j)) THEN
925             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
926          ENDIF
927       ENDDO
928     CALL xios_orchidee_send_field("treeFracPrimEver",vartmp)
929       vartmp(:)=zero
930       DO j = 2,nvm
931          IF ( .NOT.(is_c4(j)) ) THEN
932             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
933          ENDIF
934       ENDDO
935     CALL xios_orchidee_send_field("c3PftFrac",vartmp)
936       vartmp(:)=zero
937       DO j = 2,nvm
938          IF ( is_c4(j) ) THEN
939             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
940          ENDIF
941       ENDDO
942     CALL xios_orchidee_send_field("c4PftFrac",vartmp)
943     CALL xios_orchidee_send_field("rGrowth",SUM(resp_growth*veget_max,dim=2)/1e3/one_day*contfrac)
944     CALL xios_orchidee_send_field("rMaint",SUM(resp_maint*veget_max,dim=2)/1e3/one_day*contfrac)
945     CALL xios_orchidee_send_field("nppLeaf",SUM(bm_alloc(:,:,ileaf,icarbon)*veget_max,dim=2)/1e3/one_day*contfrac)
946     CALL xios_orchidee_send_field("nppWood",SUM(bm_alloc(:,:,isapabove,icarbon)*veget_max,dim=2)/1e3/one_day*contfrac)
947     CALL xios_orchidee_send_field("nppRoot",SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )*&
948          veget_max,dim=2)/1e3/one_day*contfrac)     
949
950
951    CALL histwrite_p (hist_id_stomate, 'RESOLUTION_X', itime, &
952         resolution(:,1), npts, hori_index)
953    CALL histwrite_p (hist_id_stomate, 'RESOLUTION_Y', itime, &
954         resolution(:,2), npts, hori_index)
955    CALL histwrite_p (hist_id_stomate, 'CONTFRAC', itime, &
956         contfrac(:), npts, hori_index)
957
958    CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AB', itime, &
959         litter(:,istructural,:,iabove,icarbon), npts*nvm, horipft_index)
960    CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AB', itime, &
961         litter(:,imetabolic,:,iabove,icarbon), npts*nvm, horipft_index)
962    CALL histwrite_p (hist_id_stomate, 'LITTER_STR_BE', itime, &
963         litter(:,istructural,:,ibelow,icarbon), npts*nvm, horipft_index)
964    CALL histwrite_p (hist_id_stomate, 'LITTER_MET_BE', itime, &
965         litter(:,imetabolic,:,ibelow,icarbon), npts*nvm, horipft_index)
966
967    CALL histwrite_p (hist_id_stomate, 'DEADLEAF_COVER', itime, &
968         deadleaf_cover, npts, hori_index)
969
970    CALL histwrite_p (hist_id_stomate, 'TOTAL_SOIL_CARB', itime, &
971         tot_litter_soil_carb, npts*nvm, horipft_index)
972    CALL histwrite_p (hist_id_stomate, 'CARBON_ACTIVE', itime, &
973         carbon(:,iactive,:), npts*nvm, horipft_index)
974    CALL histwrite_p (hist_id_stomate, 'CARBON_SLOW', itime, &
975         carbon(:,islow,:), npts*nvm, horipft_index)
976    CALL histwrite_p (hist_id_stomate, 'CARBON_PASSIVE', itime, &
977         carbon(:,ipassive,:), npts*nvm, horipft_index)
978
979    CALL histwrite_p (hist_id_stomate, 'T2M_MONTH', itime, &
980         t2m_month, npts, hori_index)
981    CALL histwrite_p (hist_id_stomate, 'T2M_WEEK', itime, &
982         t2m_week, npts, hori_index)
983    CALL histwrite_p (hist_id_stomate, 'TSEASON', itime, &
984         Tseason, npts, hori_index)
985    CALL histwrite_p (hist_id_stomate, 'TMIN_SPRING_TIME', itime, &
986         Tmin_spring_time, npts*nvm, horipft_index)
987    CALL histwrite_p (hist_id_stomate, 'ONSET_DATE', itime, &
988         onset_date(:,:), npts*nvm, horipft_index)
989    CALL histwrite_p (hist_id_stomate, 'FPC_MAX', itime, &
990         fpc_max, npts*nvm, horipft_index)
991    CALL histwrite_p (hist_id_stomate, 'MAXFPC_LASTYEAR', itime, &
992         maxfpc_lastyear, npts*nvm, horipft_index)
993    CALL histwrite_p (hist_id_stomate, 'HET_RESP', itime, &
994         resp_hetero(:,:), npts*nvm, horipft_index)
995    CALL histwrite_p (hist_id_stomate, 'FIREINDEX', itime, &
996         fireindex(:,:), npts*nvm, horipft_index)
997    CALL histwrite_p (hist_id_stomate, 'LITTERHUM', itime, &
998         litterhum_daily, npts, hori_index)
999    CALL histwrite_p (hist_id_stomate, 'CO2_FIRE', itime, &
1000         co2_fire, npts*nvm, horipft_index)
1001    CALL histwrite_p (hist_id_stomate, 'CO2_TAKEN', itime, &
1002         co2_to_bm, npts*nvm, horipft_index)
1003    ! land cover change
1004    CALL histwrite_p (hist_id_stomate, 'CONVFLUX', itime, &
1005         convflux, npts, hori_index)
1006    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10', itime, &
1007         cflux_prod10, npts, hori_index)
1008    CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100', itime, &
1009         cflux_prod100, npts, hori_index)
1010    CALL histwrite_p (hist_id_stomate, 'HARVEST_ABOVE', itime, &
1011         harvest_above, npts, hori_index)
1012
1013    CALL histwrite_p (hist_id_stomate, 'LAI', itime, &
1014         lai, npts*nvm, horipft_index)
1015    CALL histwrite_p (hist_id_stomate, 'VEGET_MAX', itime, &
1016         veget_max, npts*nvm, horipft_index)
1017    CALL histwrite_p (hist_id_stomate, 'NPP', itime, &
1018         npp_daily, npts*nvm, horipft_index)
1019    CALL histwrite_p (hist_id_stomate, 'GPP', itime, &
1020         gpp_daily, npts*nvm, horipft_index)
1021    CALL histwrite_p (hist_id_stomate, 'IND', itime, &
1022         ind, npts*nvm, horipft_index)
1023    CALL histwrite_p (hist_id_stomate, 'CN_IND', itime, &
1024         cn_ind, npts*nvm, horipft_index)
1025    CALL histwrite_p (hist_id_stomate, 'WOODMASS_IND', itime, &
1026         woodmass_ind, npts*nvm, horipft_index)
1027    CALL histwrite_p (hist_id_stomate, 'TOTAL_M', itime, &
1028         tot_live_biomass(:,:,icarbon), npts*nvm, horipft_index)
1029    CALL histwrite_p (hist_id_stomate, 'LEAF_M', itime, &
1030         biomass(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1031    CALL histwrite_p (hist_id_stomate, 'SAP_M_AB', itime, &
1032         biomass(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1033    CALL histwrite_p (hist_id_stomate, 'SAP_M_BE', itime, &
1034         biomass(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
1035    CALL histwrite_p (hist_id_stomate, 'HEART_M_AB', itime, &
1036         biomass(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
1037    CALL histwrite_p (hist_id_stomate, 'HEART_M_BE', itime, &
1038         biomass(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
1039    CALL histwrite_p (hist_id_stomate, 'ROOT_M', itime, &
1040         biomass(:,:,iroot,icarbon), npts*nvm, horipft_index)
1041    CALL histwrite_p (hist_id_stomate, 'FRUIT_M', itime, &
1042         biomass(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1043    CALL histwrite_p (hist_id_stomate, 'RESERVE_M', itime, &
1044         biomass(:,:,icarbres,icarbon), npts*nvm, horipft_index)
1045    CALL histwrite_p (hist_id_stomate, 'TOTAL_TURN', itime, &
1046         tot_turnover(:,:,icarbon), npts*nvm, horipft_index)
1047    CALL histwrite_p (hist_id_stomate, 'LEAF_TURN', itime, &
1048         turnover_daily(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1049    CALL histwrite_p (hist_id_stomate, 'SAP_AB_TURN', itime, &
1050         turnover_daily(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1051    CALL histwrite_p (hist_id_stomate, 'ROOT_TURN', itime, &
1052         turnover_daily(:,:,iroot,icarbon), npts*nvm, horipft_index)
1053    CALL histwrite_p (hist_id_stomate, 'FRUIT_TURN', itime, &
1054         turnover_daily(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1055    CALL histwrite_p (hist_id_stomate, 'TOTAL_BM_LITTER', itime, &
1056         tot_bm_to_litter(:,:,icarbon), npts*nvm, horipft_index)
1057    CALL histwrite_p (hist_id_stomate, 'LEAF_BM_LITTER', itime, &
1058         bm_to_litter(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1059    CALL histwrite_p (hist_id_stomate, 'SAP_AB_BM_LITTER', itime, &
1060         bm_to_litter(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1061    CALL histwrite_p (hist_id_stomate, 'SAP_BE_BM_LITTER', itime, &
1062         bm_to_litter(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
1063    CALL histwrite_p (hist_id_stomate, 'HEART_AB_BM_LITTER', itime, &
1064         bm_to_litter(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
1065    CALL histwrite_p (hist_id_stomate, 'HEART_BE_BM_LITTER', itime, &
1066         bm_to_litter(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
1067    CALL histwrite_p (hist_id_stomate, 'ROOT_BM_LITTER', itime, &
1068         bm_to_litter(:,:,iroot,icarbon), npts*nvm, horipft_index)
1069    CALL histwrite_p (hist_id_stomate, 'FRUIT_BM_LITTER', itime, &
1070         bm_to_litter(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1071    CALL histwrite_p (hist_id_stomate, 'RESERVE_BM_LITTER', itime, &
1072         bm_to_litter(:,:,icarbres,icarbon), npts*nvm, horipft_index)
1073    CALL histwrite_p (hist_id_stomate, 'MAINT_RESP', itime, &
1074         resp_maint, npts*nvm, horipft_index)
1075    CALL histwrite_p (hist_id_stomate, 'GROWTH_RESP', itime, &
1076         resp_growth, npts*nvm, horipft_index)
1077    CALL histwrite_p (hist_id_stomate, 'AGE', itime, &
1078         age, npts*nvm, horipft_index)
1079    CALL histwrite_p (hist_id_stomate, 'HEIGHT', itime, &
1080         height, npts*nvm, horipft_index)
1081    CALL histwrite_p (hist_id_stomate, 'MOISTRESS', itime, &
1082         moiavail_week, npts*nvm, horipft_index)
1083    CALL histwrite_p (hist_id_stomate, 'VCMAX', itime, &
1084         vcmax, npts*nvm, horipft_index)
1085    CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME', itime, &
1086         turnover_time, npts*nvm, horipft_index)
1087    ! land cover change
1088    CALL histwrite_p (hist_id_stomate, 'PROD10', itime, &
1089         prod10, npts*11, horip11_index)
1090    CALL histwrite_p (hist_id_stomate, 'PROD100', itime, &
1091         prod100, npts*101, horip101_index)
1092    CALL histwrite_p (hist_id_stomate, 'FLUX10', itime, &
1093         flux10, npts*10, horip10_index)
1094    CALL histwrite_p (hist_id_stomate, 'FLUX100', itime, &
1095         flux100, npts*100, horip100_index)
1096
1097    IF ( hist_id_stomate_IPCC > 0 ) THEN
1098       vartmp(:)=SUM(tot_live_biomass(:,:,icarbon)*veget_max,dim=2)/1e3*contfrac
1099       CALL histwrite_p (hist_id_stomate_IPCC, "cVeg", itime, &
1100            vartmp, npts, hori_index)
1101       vartmp(:)=SUM(tot_litter_carb*veget_max,dim=2)/1e3*contfrac
1102       CALL histwrite_p (hist_id_stomate_IPCC, "cLitter", itime, &
1103            vartmp, npts, hori_index)
1104       vartmp(:)=SUM(tot_soil_carb*veget_max,dim=2)/1e3*contfrac
1105       CALL histwrite_p (hist_id_stomate_IPCC, "cSoil", itime, &
1106            vartmp, npts, hori_index)
1107       vartmp(:)=(prod10_total + prod100_total)/1e3
1108       CALL histwrite_p (hist_id_stomate_IPCC, "cProduct", itime, &
1109            vartmp, npts, hori_index)
1110       vartmp(:)=carb_mass_variation/1e3/one_day*contfrac
1111       CALL histwrite_p (hist_id_stomate_IPCC, "cMassVariation", itime, &
1112            vartmp, npts, hori_index)
1113       vartmp(:)=SUM(lai*veget_max,dim=2)*contfrac
1114       CALL histwrite_p (hist_id_stomate_IPCC, "lai", itime, &
1115            vartmp, npts, hori_index)
1116       vartmp(:)=SUM(gpp_daily*veget_max,dim=2)/1e3/one_day*contfrac
1117       CALL histwrite_p (hist_id_stomate_IPCC, "gpp", itime, &
1118            vartmp, npts, hori_index)
1119       vartmp(:)=SUM((resp_maint+resp_growth)*veget_max,dim=2)/1e3/one_day*contfrac
1120       CALL histwrite_p (hist_id_stomate_IPCC, "ra", itime, &
1121            vartmp, npts, hori_index)
1122       vartmp(:)=SUM(npp_daily*veget_max,dim=2)/1e3/one_day*contfrac
1123       CALL histwrite_p (hist_id_stomate_IPCC, "npp", itime, &
1124            vartmp, npts, hori_index)
1125       vartmp(:)=SUM(resp_hetero*veget_max,dim=2)/1e3/one_day*contfrac
1126       CALL histwrite_p (hist_id_stomate_IPCC, "rh", itime, &
1127            vartmp, npts, hori_index)
1128       vartmp(:)=SUM(co2_fire*veget_max,dim=2)/1e3/one_day*contfrac
1129       CALL histwrite_p (hist_id_stomate_IPCC, "fFire", itime, &
1130            vartmp, npts, hori_index)
1131       vartmp(:)=harvest_above/1e3/one_day*contfrac
1132       CALL histwrite_p (hist_id_stomate_IPCC, "fHarvest", itime, &
1133            vartmp, npts, hori_index)
1134       vartmp(:)=cflux_prod_total/1e3/one_day*contfrac
1135       CALL histwrite_p (hist_id_stomate_IPCC, "fLuc", itime, &
1136            vartmp, npts, hori_index)
1137       vartmp(:)=(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
1138            &        *veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_day*contfrac
1139       CALL histwrite_p (hist_id_stomate_IPCC, "nbp", itime, &
1140            vartmp, npts, hori_index)
1141       vartmp(:)=SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*veget_max,dim=2)/1e3/one_day*contfrac
1142       CALL histwrite_p (hist_id_stomate_IPCC, "fVegLitter", itime, &
1143            vartmp, npts, hori_index)
1144       vartmp(:)=SUM(SUM(soilcarbon_input,dim=2)*veget_max,dim=2)/1e3/one_day*contfrac
1145       CALL histwrite_p (hist_id_stomate_IPCC, "fLitterSoil", itime, &
1146            vartmp, npts, hori_index)
1147       vartmp(:)=SUM(biomass(:,:,ileaf,icarbon)*veget_max,dim=2)/1e3*contfrac
1148       CALL histwrite_p (hist_id_stomate_IPCC, "cLeaf", itime, &
1149            vartmp, npts, hori_index)
1150       vartmp(:)=SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*veget_max,dim=2)/1e3*contfrac
1151       CALL histwrite_p (hist_id_stomate_IPCC, "cWood", itime, &
1152            vartmp, npts, hori_index)
1153       vartmp(:)=SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + biomass(:,:,iheartbelow,icarbon) ) &
1154            &        *veget_max,dim=2)/1e3*contfrac
1155       CALL histwrite_p (hist_id_stomate_IPCC, "cRoot", itime, &
1156            vartmp, npts, hori_index)
1157       vartmp(:)=SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*veget_max,dim=2)/1e3*contfrac
1158       CALL histwrite_p (hist_id_stomate_IPCC, "cMisc", itime, &
1159            vartmp, npts, hori_index)
1160       vartmp(:)=SUM((litter(:,istructural,:,iabove,icarbon)+litter(:,imetabolic,:,iabove,icarbon))*&
1161            veget_max,dim=2)/1e3*contfrac
1162       CALL histwrite_p (hist_id_stomate_IPCC, "cLitterAbove", itime, &
1163            vartmp, npts, hori_index)
1164       vartmp(:)=SUM((litter(:,istructural,:,ibelow,icarbon)+litter(:,imetabolic,:,ibelow,icarbon))*&
1165            veget_max,dim=2)/1e3*contfrac
1166       CALL histwrite_p (hist_id_stomate_IPCC, "cLitterBelow", itime, &
1167            vartmp, npts, hori_index)
1168       vartmp(:)=SUM(carbon(:,iactive,:)*veget_max,dim=2)/1e3*contfrac
1169       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilFast", itime, &
1170            vartmp, npts, hori_index)
1171       vartmp(:)=SUM(carbon(:,islow,:)*veget_max,dim=2)/1e3*contfrac
1172       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilMedium", itime, &
1173            vartmp, npts, hori_index)
1174       vartmp(:)=SUM(carbon(:,ipassive,:)*veget_max,dim=2)/1e3*contfrac
1175       CALL histwrite_p (hist_id_stomate_IPCC, "cSoilSlow", itime, &
1176            vartmp, npts, hori_index)
1177       DO j=1,nvm
1178          histvar(:,j)=veget_max(:,j)*contfrac(:)*100
1179       ENDDO
1180       CALL histwrite_p (hist_id_stomate_IPCC, "landCoverFrac", itime, &
1181            histvar, npts*nvm, horipft_index)
1182       !-
1183       vartmp(:)=zero
1184       DO j = 2,nvm
1185          IF (is_deciduous(j)) THEN
1186             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
1187          ENDIF
1188       ENDDO
1189       CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimDec", itime, &
1190            vartmp, npts, hori_index)
1191       !-
1192       vartmp(:)=zero
1193       DO j = 2,nvm
1194          IF (is_evergreen(j)) THEN
1195             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
1196          ENDIF
1197       ENDDO
1198       CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimEver", itime, &
1199            vartmp, npts, hori_index)
1200       !-
1201       vartmp(:)=zero
1202       DO j = 2,nvm
1203          IF ( .NOT.(is_c4(j)) ) THEN
1204             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
1205          ENDIF
1206       ENDDO
1207       CALL histwrite_p (hist_id_stomate_IPCC, "c3PftFrac", itime, &
1208            vartmp, npts, hori_index)
1209       !-
1210       vartmp(:)=zero
1211       DO j = 2,nvm
1212          IF ( is_c4(j) ) THEN
1213             vartmp(:) = vartmp(:) + veget_max(:,j)*contfrac*100
1214          ENDIF
1215       ENDDO
1216       CALL histwrite_p (hist_id_stomate_IPCC, "c4PftFrac", itime, &
1217            vartmp, npts, hori_index)
1218       !-
1219       vartmp(:)=SUM(resp_growth*veget_max,dim=2)/1e3/one_day*contfrac
1220       CALL histwrite_p (hist_id_stomate_IPCC, "rGrowth", itime, &
1221            vartmp, npts, hori_index)
1222       vartmp(:)=SUM(resp_maint*veget_max,dim=2)/1e3/one_day*contfrac
1223       CALL histwrite_p (hist_id_stomate_IPCC, "rMaint", itime, &
1224            vartmp, npts, hori_index)
1225       vartmp(:)=SUM(bm_alloc(:,:,ileaf,icarbon)*veget_max,dim=2)/1e3/one_day*contfrac
1226       CALL histwrite_p (hist_id_stomate_IPCC, "nppLeaf", itime, &
1227            vartmp, npts, hori_index)
1228       vartmp(:)=SUM(bm_alloc(:,:,isapabove,icarbon)*veget_max,dim=2)/1e3/one_day*contfrac
1229       CALL histwrite_p (hist_id_stomate_IPCC, "nppWood", itime, &
1230            vartmp, npts, hori_index)
1231       vartmp(:)=SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )*veget_max,dim=2)/1e3/one_day*contfrac
1232       CALL histwrite_p (hist_id_stomate_IPCC, "nppRoot", itime, &
1233            vartmp, npts, hori_index)
1234
1235       CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_X', itime, &
1236            resolution(:,1), npts, hori_index)
1237       CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_Y', itime, &
1238            resolution(:,2), npts, hori_index)
1239       CALL histwrite_p (hist_id_stomate_IPCC, 'CONTFRAC', itime, &
1240            contfrac(:), npts, hori_index)
1241
1242    ENDIF
1243
1244    IF (printlev>=4) WRITE(numout,*) 'Leaving stomate_lpj'
1245
1246  END SUBROUTINE StomateLpj
1247
1248
1249!! ================================================================================================================================
1250!! SUBROUTINE   : harvest
1251!!
1252!>\BRIEF        Harvest of croplands
1253!!
1254!! DESCRIPTION  : To take into account biomass harvest from crop (mainly to take
1255!! into account for the reduced litter input and then decreased soil carbon. it is a
1256!! constant (40\%) fraction of above ground biomass.
1257!!
1258!! RECENT CHANGE(S) : None
1259!!
1260!! MAIN OUTPUT VARIABLE(S): ::harvest_above the harvested biomass
1261!!
1262!! REFERENCE(S) :
1263!! - Piao, S., P. Ciais, P. Friedlingstein, N. de Noblet-Ducoudre, P. Cadule, N. Viovy, and T. Wang. 2009.
1264!!   Spatiotemporal patterns of terrestrial carbon cycle during the 20th century. Global Biogeochemical
1265!!   Cycles 23:doi:10.1029/2008GB003339.
1266!!
1267!! FLOWCHART    : None
1268!! \n
1269!_ ================================================================================================================================
1270
1271  SUBROUTINE harvest(npts, dt_days, veget_max, &
1272       bm_to_litter, turnover_daily, &
1273       harvest_above)
1274
1275  !! 0. Variable and parameter declaration
1276
1277    !! 0.1 Input variables
1278
1279    INTEGER, INTENT(in)                                    :: npts            !! Domain size (unitless)
1280    REAL(r_std), INTENT(in)                                :: dt_days         !! Time step (days)                               
1281    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: veget_max       !! new "maximal" coverage fraction of a PFT (LAI ->
1282                                                                              !! infinity) on ground @tex $(m^2 m^{-2})$ @endtex
1283   
1284   !! 0.2 Output variables
1285   
1286   !! 0.3 Modified variables
1287
1288    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! [DISPENSABLE] conversion of biomass to litter
1289                                                                                     !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1290    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily   !! Turnover rates
1291                                                                               !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1292    REAL(r_std), DIMENSION(npts), INTENT(inout)            :: harvest_above    !! harvest above ground biomass for agriculture
1293                                                                               !! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
1294    !! 0.4 Local variables
1295
1296    INTEGER(i_std)                                         :: i, j, k, l, m    !! indices                       
1297    REAL(r_std)                                            :: above_old        !! biomass of previous time step
1298                                                                               !! @tex $(gC m^{-2})$ @endtex
1299!_ ================================================================================================================================
1300
1301  !! 1. Yearly initialisation
1302
1303    above_old             = zero
1304    harvest_above         = zero
1305
1306    DO i = 1, npts
1307       DO j = 1,nvm
1308          IF (.NOT. natural(j)) THEN
1309             above_old = turnover_daily(i,j,ileaf,icarbon) + turnover_daily(i,j,isapabove,icarbon) + &
1310                  &       turnover_daily(i,j,iheartabove,icarbon) + turnover_daily(i,j,ifruit,icarbon) + &
1311                  &       turnover_daily(i,j,icarbres,icarbon) + turnover_daily(i,j,isapbelow,icarbon) + &
1312                  &       turnover_daily(i,j,iheartbelow,icarbon) + turnover_daily(i,j,iroot,icarbon)
1313
1314             turnover_daily(i,j,ileaf,icarbon) = turnover_daily(i,j,ileaf,icarbon)*frac_turnover_daily
1315             turnover_daily(i,j,isapabove,icarbon) = turnover_daily(i,j,isapabove,icarbon)*frac_turnover_daily
1316             turnover_daily(i,j,isapbelow,icarbon) = turnover_daily(i,j,isapbelow,icarbon)*frac_turnover_daily
1317             turnover_daily(i,j,iheartabove,icarbon) = turnover_daily(i,j,iheartabove,icarbon)*frac_turnover_daily
1318             turnover_daily(i,j,iheartbelow,icarbon) = turnover_daily(i,j,iheartbelow,icarbon)*frac_turnover_daily
1319             turnover_daily(i,j,iroot,icarbon) = turnover_daily(i,j,iroot,icarbon)*frac_turnover_daily
1320             turnover_daily(i,j,ifruit,icarbon) = turnover_daily(i,j,ifruit,icarbon)*frac_turnover_daily
1321             turnover_daily(i,j,icarbres,icarbon) = turnover_daily(i,j,icarbres,icarbon)*frac_turnover_daily
1322             harvest_above(i)  = harvest_above(i) + veget_max(i,j) * above_old *(un - frac_turnover_daily)
1323          ENDIF
1324       ENDDO
1325    ENDDO
1326
1327!!$    harvest_above = harvest_above
1328  END SUBROUTINE harvest
1329END MODULE stomate_lpj
Note: See TracBrowser for help on using the repository browser.