source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/stomate.f90 @ 5816

Last change on this file since 5816 was 5816, checked in by jinfeng.chang, 5 years ago

copy ORCHIDEE-GMv3.2 for publication

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 242.1 KB
Line 
1! =================================================================================================================================
2! MODULE       : stomate
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       Groups the subroutines that: (1) initialize all variables in
10!! stomate, (2) read and write forcing files of stomate and the soil component,
11!! (3) aggregates and convert variables to handle the different time steps
12!! between sechiba and stomate, (4) call subroutines that govern major stomate
13!! processes (litter, soil, and vegetation dynamics) and (5) structures these tasks
14!! in stomate_main
15!!
16!!\n DESCRIPTION : None
17!!
18!! RECENT CHANGE(S) : None
19!!
20!! REFERENCE(S) : None
21!!
22!! SVN :
23!! $HeadURL$
24!! $Date$
25!! $Revision$
26!! \n
27!_ ================================================================================================================================
28
29MODULE stomate
30
31  ! Modules used:
32  USE netcdf
33  USE defprec
34  USE grid
35  USE constantes
36  USE constantes_soil
37  USE pft_parameters
38  USE stomate_io
39  USE stomate_data
40  USE stomate_season
41  USE stomate_lpj
42  USE stomate_litter
43  USE stomate_vmax
44  USE stomate_soilcarbon
45  USE stomate_resp
46  USE mod_orchidee_para
47  USE ioipsl_para 
48  USE xios_orchidee
49
50  USE matrix_resolution
51 
52  IMPLICIT NONE
53
54  ! Private & public routines
55
56  PRIVATE
57  PUBLIC stomate_main,stomate_clear,init_forcing, stomate_forcing_read, stomate_initialize, stomate_finalize
58
59
60  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: biomass              !! Biomass per ground area @tex $(gC m^{-2})$ @endtex
61!$OMP THREADPRIVATE(biomass)
62  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: veget_cov_max        !! Maximal fractional coverage: maximum share of a pixel
63                                                                         !! taken by a PFT
64!$OMP THREADPRIVATE(veget_cov_max)
65  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: ind                  !! Vegetation density, number of individuals per unit
66                                                                         !! ground area @tex $(m^{-2})$ @endtex
67!$OMP THREADPRIVATE(ind)
68  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: age                  !! Age of PFT it normalized by biomass - can increase and
69                                                                         !! decrease - (years)
70!$OMP THREADPRIVATE(age)
71  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: adapted              !! Winter too cold for PFT to survive (0-1, unitless)
72!$OMP THREADPRIVATE(adapted)
73  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: regenerate           !! Winter sufficiently cold to produce viable seeds
74                                                                         !! (0-1, unitless)
75!$OMP THREADPRIVATE(regenerate)
76  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: everywhere           !! Is the PFT everywhere in the grid box or very localized
77                                                                         !! (after its intoduction)
78!$OMP THREADPRIVATE(everywhere)
79  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: fireindex            !! Probability of fire (unitless)
80!$OMP THREADPRIVATE(fireindex)
81  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: veget_lastlight      !! Vegetation fractions (on ground) after last light
82                                                                         !! competition (unitless)
83!$OMP THREADPRIVATE(veget_lastlight)
84  REAL(r_std), ALLOCATABLE,SAVE,DIMENSION(:,:)   :: fpc_max              !! "maximal" coverage fraction of a grid box (LAI ->
85                                                                         !! infinity) on ground. [??CHECK??] It's set to zero here,
86                                                                         !! and then is used once in lpj_light.f90 to test if
87                                                                         !! fpc_nat is greater than it. Something seems missing
88!$OMP THREADPRIVATE(fpc_max)
89  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:)        :: PFTpresent           !! PFT exists (equivalent to veget > 0 for natural PFTs)
90!$OMP THREADPRIVATE(PFTpresent)
91  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:)        :: senescence           !! The PFT is senescent
92!$OMP THREADPRIVATE(senescence)
93  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:)        :: begin_leaves         !! Signal to start putting leaves on (true/false)
94!$OMP THREADPRIVATE(begin_leaves)
95  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:)        :: need_adjacent        !! This PFT needs to be in present in an adjacent gridbox
96                                                                         !! if it is to be introduced in a new gridbox
97!$OMP THREADPRIVATE(need_adjacent)
98!--
99  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: humrel_daily         !! Daily plant available water -root profile weighted
100                                                                         !! (0-1, unitless)
101!$OMP THREADPRIVATE(humrel_daily)
102  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: humrel_week          !! "Weekly" plant available water -root profile weighted
103                                                                         !! (0-1, unitless)
104!$OMP THREADPRIVATE(humrel_week)
105  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: humrel_month         !! "Monthly" plant available water -root profile weighted
106                                                                         !! (0-1, unitless)
107!$OMP THREADPRIVATE(humrel_month)
108  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxhumrel_lastyear   !! Last year's max plant available water -root profile
109                                                                         !! weighted (0-1, unitless)
110!$OMP THREADPRIVATE(maxhumrel_lastyear)
111  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxhumrel_thisyear   !! This year's max plant available water -root profile
112                                                                         !! weighted (0-1, unitless)
113!$OMP THREADPRIVATE(maxhumrel_thisyear)
114  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: minhumrel_lastyear   !! Last year's min plant available water -root profile
115                                                                         !! weighted (0-1, unitless) 
116!$OMP THREADPRIVATE(minhumrel_lastyear)
117  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: minhumrel_thisyear   !! This year's minimum plant available water -root profile
118                                                                         !! weighted (0-1, unitless)
119!$OMP THREADPRIVATE(minhumrel_thisyear)
120!--- 
121  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_daily            !! Daily air temperature at 2 meter (K)
122!$OMP THREADPRIVATE(t2m_daily)
123
124  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: Tseason              !! "seasonal" 2 meter temperatures (K)
125!$OMP THREADPRIVATE(Tseason)
126  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: Tseason_length       !! temporary variable to calculate Tseason
127!$OMP THREADPRIVATE(Tseason_length)
128  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: Tseason_tmp          !! temporary variable to calculate Tseason
129!$OMP THREADPRIVATE(Tseason_tmp)
130  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: Tmin_spring_time     !! Number of days after begin_leaves (leaf onset)
131!$OMP THREADPRIVATE(Tmin_spring_time)
132  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: onset_date           !! Date in the year at when the leaves started to grow(begin_leaves), only for diagnostics.
133!$OMP THREADPRIVATE(onset_date)
134  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_week             !! Mean "weekly" (default 7 days) air temperature at 2
135                                                                         !! meter (K) 
136!$OMP THREADPRIVATE(t2m_week)
137  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_month            !! Mean "monthly" (default 20 days) air temperature at 2
138                                                                         !! meter (K)
139!$OMP THREADPRIVATE(t2m_month)
140  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_longterm         !! Mean "Long term" (default 3 years) air temperature at
141                                                                         !! 2 meter (K)
142!$OMP THREADPRIVATE(t2m_longterm)
143  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_min_daily        !! Daily minimum air temperature at 2 meter (K)
144!$OMP THREADPRIVATE(t2m_min_daily)
145  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: tsurf_daily          !! Daily surface temperatures (K)
146!$OMP THREADPRIVATE(tsurf_daily)
147!---
148  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: precip_daily         !! Daily precipitations sum @tex $(mm day^{-1})$ @endtex
149!$OMP THREADPRIVATE(precip_daily)
150  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: precip_lastyear      !! Last year's annual precipitation sum
151                                                                         !! @tex $??(mm year^{-1})$ @endtex
152!$OMP THREADPRIVATE(precip_lastyear)
153  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: precip_thisyear      !! This year's annual precipitation sum
154                                                                         !! @tex $??(mm year^{-1})$ @endtex
155!$OMP THREADPRIVATE(precip_thisyear)
156!---
157  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: soilhum_daily        !! Daily soil humidity (0-1, unitless)
158!$OMP THREADPRIVATE(soilhum_daily)
159  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: soilhum_month        !! Soil humidity - integrated over a month (0-1, unitless)
160!$OMP THREADPRIVATE(soilhum_month)
161  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: tsoil_daily          !! Daily soil temperatures (K)
162!$OMP THREADPRIVATE(tsoil_daily)
163  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: tsoil_month          !! Soil temperatures at each soil layer integrated over a
164                                                                         !! month (K)
165!$OMP THREADPRIVATE(tsoil_month)
166!---
167  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: litterhum_daily      !! Daily litter humidity (0-1, unitless)
168!$OMP THREADPRIVATE(litterhum_daily)
169!---
170  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: control_moist        !! Moisture control of heterotrophic respiration
171                                                                         !! (0-1, unitless)
172!$OMP THREADPRIVATE(control_moist)
173  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: control_temp         !! Temperature control of heterotrophic respiration at the
174                                                                         !! different soil levels (0-1, unitless)
175!$OMP THREADPRIVATE(control_temp)
176  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: control_moist_daily  !! Moisture control of heterotrophic respiration daily
177                                                                         !! (0-1, unitless)
178!$OMP THREADPRIVATE(control_moist_daily)
179  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: control_temp_daily   !! Temperature control of heterotrophic respiration, above
180                                                                         !! and below daily (0-1, unitless)
181!$OMP THREADPRIVATE(control_temp_daily)
182!---
183  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gdd_init_date        !! inital date for gdd count
184!$OMP THREADPRIVATE(gdd_init_date)
185
186  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gdd_from_growthinit  !! gdd from beginning of season (C)
187!$OMP THREADPRIVATE(gdd_from_growthinit)
188  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: gdd0_lastyear        !! Last year's annual Growing Degree Days,
189                                                                         !! threshold 0 deg C (K)
190!$OMP THREADPRIVATE(gdd0_lastyear)
191  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: gdd0_thisyear        !! This year's annual Growing Degree Days,
192                                                                         !! threshold 0 deg C (K)
193!$OMP THREADPRIVATE(gdd0_thisyear)
194  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gdd_m5_dormance      !! Growing degree days for onset of growing season,
195                                                                         !! threshold -5 deg C (K)
196!$OMP THREADPRIVATE(gdd_m5_dormance)
197  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gdd_midwinter        !! Growing degree days for onset of growing season,
198                                                                         !! since midwinter (K)
199!$OMP THREADPRIVATE(gdd_midwinter)
200  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: ncd_dormance         !! Number of chilling days since leaves were lost (days)
201!$OMP THREADPRIVATE(ncd_dormance)
202  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: ngd_minus5           !! Number of growing days, threshold -5 deg C (days)
203!$OMP THREADPRIVATE(ngd_minus5)
204  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: hum_min_dormance     !! Minimum moisture during dormance (0-1, unitless)
205!$OMP THREADPRIVATE(hum_min_dormance)
206!---
207  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gpp_daily            !! Daily gross primary productivity per ground area
208                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
209!$OMP THREADPRIVATE(gpp_daily)
210  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gpp_week             !! Mean "weekly" (default 7 days) GPP 
211                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
212!$OMP THREADPRIVATE(gpp_week)
213  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxgppweek_lastyear  !! Last year's maximum "weekly" GPP 
214                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
215!$OMP THREADPRIVATE(maxgppweek_lastyear)
216  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxgppweek_thisyear  !! This year's maximum "weekly" GPP 
217                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex 
218!$OMP THREADPRIVATE(maxgppweek_thisyear)
219!---
220  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: npp_daily            !! Daily net primary productivity per ground area
221                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
222!$OMP THREADPRIVATE(npp_daily)
223  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: npp_longterm         !! "Long term" (default 3 years) net primary productivity
224                                                                         !! per ground area 
225                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex   
226!$OMP THREADPRIVATE(npp_longterm)
227  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: npp_equil            !! Equilibrium NPP written to forcesoil
228                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
229!$OMP THREADPRIVATE(npp_equil)
230  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: npp_tot              !! Total NPP written to forcesoil
231                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
232!$OMP THREADPRIVATE(npp_tot)
233!---
234  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: resp_maint_part_radia!! Maintenance respiration of different plant parts per
235                                                                         !! total ground area at Sechiba time step 
236                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
237!$OMP THREADPRIVATE(resp_maint_part_radia)
238  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: resp_maint_part      !! Maintenance respiration of different plant parts per
239                                                                         !! total ground area at Stomate time step
240                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
241!$OMP THREADPRIVATE(resp_maint_part)
242  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_maint_radia     !! Maintenance respiration per ground area at Sechiba time
243                                                                         !! step   
244                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
245!$OMP THREADPRIVATE(resp_maint_radia)
246  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_maint_d         !! Maintenance respiration per ground area at Stomate time
247                                                                         !! step 
248                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
249!$OMP THREADPRIVATE(resp_maint_d)
250  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_growth_d        !! Growth respiration per ground area
251                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
252!$OMP THREADPRIVATE(resp_growth_d)
253  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_hetero_d        !! Heterotrophic respiration per ground area
254                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
255!$OMP THREADPRIVATE(resp_hetero_d)
256  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_hetero_radia    !! Heterothrophic respiration per ground area at Sechiba
257                                                                         !! time step
258                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
259!$OMP THREADPRIVATE(resp_hetero_radia)
260!---
261  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)     :: turnover_time       !! Turnover time of grasses
262                                                                         !! @tex $(dt_stomate^{-1})$ @endtex
263!$OMP THREADPRIVATE(turnover_time)
264  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: turnover_daily      !! Senescence-driven turnover (better: mortality) of
265                                                                         !! leaves and roots 
266                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
267!$OMP THREADPRIVATE(turnover_daily)
268  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: turnover_littercalc !! Senescence-driven turnover (better: mortality) of
269                                                                         !! leaves and roots at Sechiba time step
270                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
271!$OMP THREADPRIVATE(turnover_littercalc)
272  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: turnover_longterm   !! "Long term" (default 3 years) senescence-driven
273                                                                         !! turnover (better: mortality) of leaves and roots
274                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
275!$OMP THREADPRIVATE(turnover_longterm)
276  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: bm_to_litter        !! Background (not senescence-driven) mortality of biomass
277                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
278!$OMP THREADPRIVATE(bm_to_litter)
279  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: bm_to_littercalc    !! conversion of biomass to litter per ground area at
280                                                                         !! Sechiba time step
281                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
282!$OMP THREADPRIVATE(bm_to_littercalc)
283  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: dead_leaves          !! Metabolic and structural pools of dead leaves on ground
284                                                                         !! per PFT @tex $(gC m^{-2})$ @endtex
285!$OMP THREADPRIVATE(dead_leaves)
286  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:):: litter             !! Above and below ground metabolic and structural litter
287                                                                         !! per ground area
288                                                                         !! @tex $(gC m^{-2})$ @endtex
289!$OMP THREADPRIVATE(litter)
290  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: litterpart           !! Fraction of litter above the ground belonging to
291                                                                         !! different litter pools (unitless)
292!$OMP THREADPRIVATE(litterpart)
293  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: firelitter           !! Total litter above the ground that could potentially
294                                                                         !! burn @tex $(gC m^{-2})$ @endtex
295!$OMP THREADPRIVATE(firelitter)
296  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: soilcarbon_input     !! Quantity of carbon going into carbon pools from litter
297                                                                         !! decomposition per ground area  at Sechiba time step
298                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
299!$OMP THREADPRIVATE(soilcarbon_input)
300  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: soilcarbon_input_daily !! Daily quantity of carbon going into carbon pools from
301                                                                           !! litter decomposition per ground area
302                                                                           !! @tex $(gC m^{-2} day^{-1})$ @endtex
303!$OMP THREADPRIVATE(soilcarbon_input_daily)
304  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: carbon               !! Soil carbon pools per ground area: active, slow, or
305                                                                         !! passive, @tex $(gC m^{-2})$ @endtex
306!$OMP THREADPRIVATE(carbon)
307  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: lignin_struc         !! Ratio Lignine/Carbon in structural litter for above and
308                                                                         !! below ground compartments (unitless)
309!$OMP THREADPRIVATE(lignin_struc)
310  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: lm_lastyearmax       !! Last year's maximum leaf mass per ground area for each
311                                                                         !! PFT @tex $(gC m^{-2})$ @endtex 
312!$OMP THREADPRIVATE(lm_lastyearmax)
313  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: lm_thisyearmax       !! This year's maximum leaf mass per ground area for each
314                                                                         !! PFT @tex $(gC m^{-2})$ @endtex 
315!$OMP THREADPRIVATE(lm_thisyearmax)
316  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxfpc_lastyear      !! Last year's maximum fpc for each natural PFT, on ground
317                                                                         !! [??CHECK] fpc but this ones look ok (computed in
318                                                                         !! season, used in light)??
319!$OMP THREADPRIVATE(maxfpc_lastyear)
320  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxfpc_thisyear      !! This year's maximum fpc for each PFT, on ground (see
321                                                                         !! stomate_season), [??CHECK] fpc but this ones look ok
322                                                                         !! (computed in season, used in light)??
323!$OMP THREADPRIVATE(maxfpc_thisyear)
324!---
325  REAL(r_std), ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: leaf_age             !! Age of different leaf classes (days)
326!$OMP THREADPRIVATE(leaf_age)
327  REAL(r_std), ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: leaf_frac            !! PFT fraction of leaf mass in leaf age class (0-1,
328                                                                         !! unitless)
329!$OMP THREADPRIVATE(leaf_frac)
330  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: when_growthinit      !! Days since beginning of growing season (days)
331!$OMP THREADPRIVATE(when_growthinit)
332  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: herbivores           !! Time constant of probability of a leaf to be eaten by a
333                                                                         !! herbivore (days)
334!$OMP THREADPRIVATE(herbivores)
335  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: RIP_time             !! How much time ago was the PFT eliminated for the last
336                                                                         !! time (year)
337!$OMP THREADPRIVATE(RIP_time)
338  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: time_hum_min         !! Time elapsed since strongest moisture limitation (days)
339!$OMP THREADPRIVATE(time_hum_min)
340!---
341  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: clay_fm              !! Soil clay content (0-1, unitless), parallel computing
342  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: clay_fm_g            !! Soil clay content (0-1, unitless), parallel computing
343  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: precip_fm            !! Daily precipitations sum @tex $(mm day^{-1})$ @endtex,
344                                                                         !! parallel computing
345  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: precip_fm_g          !! Daily precipitations sum @tex $(mm day^{-1})$ @endtex,
346                                                                         !! parallel computing
347  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: litterhum_daily_fm   !! Daily relative humidity of litter (0-1, unitless),
348                                                                         !! parallel computing
349  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: litterhum_daily_fm_g !! Daily relative humidity of litter (0-1, unitless),
350                                                                         !! parallel computing
351  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: t2m_daily_fm         !! Daily air temperature at 2 meter (K), parallel
352                                                                         !! computing
353  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: t2m_daily_fm_g       !! Daily air temperature at 2 meter (K), parallel
354                                                                         !! computing
355  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: t2m_min_daily_fm     !! Daily minimum air temperature at 2 meter (K),
356                                                                         !! parallel computing
357  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: t2m_min_daily_fm_g   !! Daily minimum air temperature at 2 meter (K),
358                                                                         !! parallel computing
359  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: tsurf_daily_fm       !! Daily surface temperatures (K), parallel
360                                                                         !! computing
361  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: tsurf_daily_fm_g     !! Daily surface temperatures (K), parallel
362                                                                         !! computing
363  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: tsoil_daily_fm       !! Daily soil temperatures (K), parallel computing
364  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: tsoil_daily_fm_g     !! Daily soil temperatures (K), parallel computing
365  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: soilhum_daily_fm     !! Daily soil humidity (0-1, unitless), parallel computing
366  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: soilhum_daily_fm_g   !! Daily soil humidity (0-1, unitless), parallel computing
367  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: humrel_daily_fm      !! Daily relative humidity of atmosphere (0-1, unitless),
368                                                                         !! parallel computing
369  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: humrel_daily_fm_g    !! Daily relative humidity of atmosphere (0-1, unitless),
370                                                                         !! parallel computing
371  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: gpp_daily_fm         !! Daily gross primary productivity per ground area
372                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex,
373                                                                         !! parallel computing
374  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: gpp_daily_fm_g       !! Daily gross primary productivity per ground area
375                                                                         !! @tex $(gC m^{-2} day^{-1})$ @endtex,
376                                                                         !! parallel computing
377  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: veget_fm             !! Vegetation coverage taking into account non-biological
378                                                                         !! coverage (unitless), parallel computing
379  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: veget_fm_g           !! Vegetation coverage taking into account non-biological
380                                                                         !! coverage (unitless), parallel computing
381  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: veget_max_fm         !! Maximum vegetation coverage taking into account
382                                                                         !! non-biological coverage (unitless), parallel computing
383  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: veget_max_fm_g       !! Maximum vegetation coverage taking into account none
384                                                                         !! biological coverage (unitless), parallel computing
385  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: lai_fm               !! Leaf area index @tex $@tex $(m^2 m^{-2})$ @endtex$ @endtex,
386                                                                         !! parallel computing
387  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: lai_fm_g             !! Leaf area index @tex $@tex $(m^2 m^{-2})$ @endtex$ @endtex,
388                                                                         !! parallel computing
389!---
390  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: co2_fire             !! Carbon emitted to the atmosphere by burning living
391                                                                         !! and dead biomass
392                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
393!$OMP THREADPRIVATE(co2_fire)
394  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: co2_to_bm_dgvm       !! Psuedo-photosynthesis,C used to provide seedlings with
395                                                                         !! an initial biomass, arbitrarily removed from the
396                                                                         !! atmosphere 
397                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
398!$OMP THREADPRIVATE(co2_to_bm_dgvm)
399  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: co2_flux_daily       !! Daily net CO2 flux between atmosphere and biosphere
400                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
401                                                                         !! [??CHECK] sign convention?
402!$OMP THREADPRIVATE(co2_flux_daily)
403  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: co2_flux_monthly     !! Monthly net CO2 flux between atmosphere and biosphere
404                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
405                                                                         !! [??CHECK] sign convention?
406!$OMP THREADPRIVATE(co2_flux_monthly)
407  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: prod10               !! Wood products remaining in the 10 year-turnover pool
408                                                                         !! after the annual release for each compartment
409                                                                         !! @tex $(gC m^{-2})$ @endtex   
410                                                                         !! (0:10 input from year of land cover change),
411                                                                         !! dimension(#pixels,0:10 years
412!$OMP THREADPRIVATE(prod10)
413  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: prod100              !! Wood products remaining in the 100 year-turnover pool
414                                                                         !! after the annual release for each compartment
415                                                                         !! @tex $(gC m^{-2})$ @endtex 
416                                                                         !! (0:100 input from year of land cover change),
417                                                                         !! dimension(#pixels,0:100 years)
418!$OMP THREADPRIVATE(prod100)
419  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: flux10               !! Wood decomposition from the 10 year-turnover pool
420                                                                         !! compartments
421                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
422                                                                         !! dimension(#pixels,0:10) 
423!$OMP THREADPRIVATE(flux10)
424  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: flux100              !! Wood decomposition from the 100 year-turnover pool
425                                                                         !! compartments
426                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
427                                                                         !! dimension(#pixels,0:100)
428!$OMP THREADPRIVATE(flux100)
429  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: convflux             !! Release during first year following land cover change
430                                                                         !! (paper, burned, etc...)
431                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex 
432!$OMP THREADPRIVATE(convflux)
433  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: cflux_prod10         !! Total annual release from the 10 year-turnover pool
434                                                                         !! sum of flux10 
435                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
436!$OMP THREADPRIVATE(cflux_prod10)
437  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: cflux_prod100        !! Total annual release from the 100 year-turnover pool
438                                                                         !! sum of flux100
439                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
440!$OMP THREADPRIVATE(cflux_prod100)
441  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: harvest_above        !! Harvest of above ground biomass for agriculture -not
442                                                                         !! just from land use change
443                                                                         !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex
444!$OMP THREADPRIVATE(harvest_above)
445  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: carb_mass_total      !! Total on-site and off-site C pool
446                                                                         !! @tex $(??gC m^{-2})$ @endtex                       
447!$OMP THREADPRIVATE(carb_mass_total)
448!---
449  REAL(r_std), SAVE                              :: tau_longterm
450!$OMP THREADPRIVATE(tau_longterm)
451  REAL(r_std),SAVE                               :: dt_days=zero         !! Time step of STOMATE (days)
452!$OMP THREADPRIVATE(dt_days)
453  INTEGER(i_std),SAVE                            :: date=0               !! Date (days)
454!$OMP THREADPRIVATE(date)
455  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:)   :: nforce               !! Number of states calculated for the soil forcing
456                                                                         !! variables (unitless), dimension(::nparan*::nbyear) both
457                                                                         !! given in the run definition file   
458!$OMP THREADPRIVATE(nforce)
459  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:)   :: isf                  !! Index for number of time steps that can be stored in
460                                                                         !! memory (unitless), dimension (#nsfm)
461!$OMP THREADPRIVATE(isf)
462  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:)   :: nf_cumul             !! Number of years over which the average is calculated in
463                                                                         !! forcesoil when cumul flag is set, dimension (#nsft)
464                                                                         !! [??CHECK] definition the dimension is number of
465                                                                         !! timesteps in a year?
466!$OMP THREADPRIVATE(nf_cumul)
467  INTEGER(i_std), SAVE                           :: spinup_period        !! Period of years used to calculate the resolution of the system for spinup analytic.
468                                                                         !! This period correspond in most cases to the period of years of forcing data used
469  INTEGER,PARAMETER                              :: r_typ = nf90_real4   !! Specify data format (server dependent)
470  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:)          :: nf_written           !! Flag indicating whether the forcing data have been
471                                                                         !! written
472!$OMP THREADPRIVATE(nf_written)
473!---
474  LOGICAL, SAVE                                  :: do_slow=.FALSE.      !! Flag that determines whether stomate_accu calculates
475                                                                         !! the sum(do_slow=.FALSE.) or the mean
476                                                                         !! (do_slow=.TRUE.)
477!$OMP THREADPRIVATE(do_slow)
478  LOGICAL, SAVE                                  :: EndOfYear=.FALSE.    !! Update annual variables? This variable must be
479                                                                         !! .TRUE. once a year
480!$OMP THREADPRIVATE(EndOfYear)
481  LOGICAL, SAVE                                  :: EndOfMonth=.FALSE.   !! Update monthly variables? This variable must be
482                                                                         !!.TRUE. once a month
483!$OMP THREADPRIVATE(EndOfMonth)
484  LOGICAL, SAVE                                  :: l_first_stomate = .TRUE.!! Is this the first call of stomate?
485!$OMP THREADPRIVATE(l_first_stomate)
486  LOGICAL, SAVE                                  :: cumul_forcing=.FALSE.!! flag for cumul of forcing if teststomate
487!$OMP THREADPRIVATE(cumul_forcing)
488  LOGICAL, SAVE                                  :: cumul_Cforcing=.FALSE.  !! Flag, if internal parameter cumul_Cforcing is
489                                                                            !! TRUE then ::nbyear (defined in run definition
490                                                                            !! file will be forced to 1 later in this module. If
491                                                                            !! FALSE the mean over ::nbyear is written in forcesoil
492!$OMP THREADPRIVATE(cumul_Cforcing)
493!---   
494  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: harvest_above_monthly   !! [??CHECK] post-processing - should be removed?
495!$OMP THREADPRIVATE(harvest_above_monthly)
496  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: cflux_prod_monthly      !! [??CHECK] post-processing - should be removed?
497!$OMP THREADPRIVATE(cflux_prod_monthly)
498!---
499  INTEGER(i_std), SAVE                               :: global_years        !! Global counter of years (year)
500!$OMP THREADPRIVATE(global_years)
501  LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:)           :: ok_equilibrium      !! Logical array marking the points where the resolution is ok
502                                                                            !! (true/false)
503!$OMP THREADPRIVATE(ok_equilibrium)
504  LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:)           :: carbon_eq           !! Logical array to mark the carbon pools at equilibrium ?
505                                                                            !! If true, the job stops. (true/false)
506!$OMP THREADPRIVATE(carbon_eq)
507  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)       :: nbp_accu            !! Accumulated Net Biospheric Production over the year (gC.m^2 )
508!$OMP THREADPRIVATE(nbp_accu)
509  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)       :: nbp_flux            !! Net Biospheric Production (gC.m^2.day^{-1})
510!$OMP THREADPRIVATE(nbp_flux)
511  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:)       :: matrixA             !! Matrix containing the fluxes between the carbon pools
512                                                                            !! per sechiba time step
513                                                                            !! @tex $(gC.m^2.day^{-1})$ @endtex
514!$OMP THREADPRIVATE(matrixA)
515  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)         :: vectorB             !! Vector containing the litter increase per sechiba time step
516                                                                            !! @tex $(gC m^{-2})$ @endtex
517!$OMP THREADPRIVATE(vectorB)
518  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: MatrixV             !! Matrix containing the accumulated values of matrixA
519!$OMP THREADPRIVATE(MatrixV)
520  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: VectorU             !! Matrix containing the accumulated values of VectorB
521!$OMP THREADPRIVATE(VectorU)
522  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: MatrixW             !! Matrix containing the opposite of matrixA
523!$OMP THREADPRIVATE(MatrixW)
524  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: previous_stock      !! Array containing the carbon stock calculated by the analytical
525                                                                            !! method in the previous resolution
526!$OMP THREADPRIVATE(previous_stock)
527  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: current_stock       !! Array containing the carbon stock calculated by the analytical
528                                                                            !! method in the current resolution
529!$OMP THREADPRIVATE(current_stock)
530  REAL(r_std), SAVE                                  :: eps_carbon          !! Stopping criterion for carbon pools (unitless,0-1)
531!$OMP THREADPRIVATE(eps_carbon)
532! gmjc
533!gmjc for grazing litter
534  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: litter_avail
535  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: litter_not_avail
536  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: litter_avail_frac
537!end gmjc
538!  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: resp_hetero_soil_part
539!  ! heterotrophic respiration for soil pool (gC/day/m**2 of total ground)
540!  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)    :: resp_hetero_soil_d
541!  ! heterotrophic respiration for litter (gC/day/m**2 of total ground)
542!  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_hetero_litter_d
543  ! N fertilization factor controlling Vcmax and SLA
544  REAL(r_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: N_limfert
545  ! specific leaf area (m2/gC)
546  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: sla_calc
547  ! "14days" 2 meter temperatures (K)
548  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)        :: t2m_14
549  REAL(r_std ), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: wshtotsum !! yield from autogestion=1 offer import_yield for autogestion=2
550  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sr_ugb !!Optimised stocking density (animal m-2)
551  REAL(r_std)     , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: compt_ugb !!Counter of grazing days (d)
552  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: nb_ani
553  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: grazed_frac
554  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: import_yield
555  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)      :: sla_age1
556  REAL(r_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: when_growthinit_cut
557  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: nb_grazingdays
558  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)        :: snowfall_daily
559  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)        :: snowmass_daily
560!JCADD top 5 layer grassland soil moisture for grazing
561  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: tmc_topgrass_daily      !! Daily top 5 layer soil moisture (m^3 m^-3)
562  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: after_snow
563  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: after_wet
564  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: wet1day
565  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: wet2day
566! end gmjc
567  REAL(r_std),SAVE                                   :: dt_forcesoil        !! Time step of soil forcing file (days)
568!$OMP THREADPRIVATE(dt_forcesoil)
569  INTEGER(i_std),PARAMETER                           :: nparanmax=366       !! Maximum number of time steps per year for forcesoil
570  INTEGER(i_std),SAVE                                :: nparan              !! Number of time steps per year for forcesoil read from run definition (unitless)
571!$OMP THREADPRIVATE(nparan)
572  INTEGER(i_std),SAVE                                :: nbyear=1            !! Number of years saved for forcesoil (unitless)
573!$OMP THREADPRIVATE(nbyear)
574  INTEGER(i_std),SAVE                                :: iatt                !! Time step of forcing of soil processes (iatt = 1 to ::nparan*::nbyear)
575!$OMP THREADPRIVATE(iatt)
576  INTEGER(i_std),SAVE                                :: iatt_old=1          !! Previous ::iatt
577!$OMP THREADPRIVATE(iatt_old)
578  INTEGER(i_std),SAVE                                :: nsfm                !! Number of time steps that can be stored in memory (unitless)
579!$OMP THREADPRIVATE(nsfm)
580  INTEGER(i_std),SAVE                                :: nsft                !! Number of time steps in a year (unitless)
581!$OMP THREADPRIVATE(nsft)
582  INTEGER(i_std),SAVE                                :: iisf                !! Current pointer for teststomate (unitless)
583!$OMP THREADPRIVATE(iisf)
584  CHARACTER(LEN=100), SAVE                           :: forcing_name        !! Name of forcing file 1
585!$OMP THREADPRIVATE(forcing_name)
586  CHARACTER(LEN=100), SAVE                           :: Cforcing_name       !! Name of forcing file 2
587!$OMP THREADPRIVATE(Cforcing_name)
588  INTEGER(i_std),SAVE                                :: Cforcing_id         !! File identifer of file 2
589!$OMP THREADPRIVATE(Cforcing_id)   
590  INTEGER(i_std),PARAMETER                           :: ndm = 10            !! Maximum number of dimensions (unitless)
591
592 
593PUBLIC clay_fm, humrel_daily_fm, litterhum_daily_fm, t2m_daily_fm, &
594   & t2m_min_daily_fm, tsurf_daily_fm, tsoil_daily_fm, soilhum_daily_fm, &
595   & precip_fm, gpp_daily_fm, veget_fm, veget_max_fm, lai_fm
596PUBLIC  dt_days, date, do_slow, EndOfYear
597PUBLIC isf, nf_written
598
599CONTAINS
600 
601
602!! ================================================================================================================================
603!! SUBROUTINE   : stomate_initialize
604!!
605!>\BRIEF        Initialization routine for stomate module.
606!!
607!! DESCRIPTION  : Initialization routine for stomate module. Read options from parameter file, allocate variables, read variables
608!!                from restart file and initialize variables if necessary.
609!!               
610!! \n
611!_ ================================================================================================================================
612
613SUBROUTINE stomate_initialize &
614        (kjit,           kjpij,             kjpindex,                        &
615         rest_id_stom,   hist_id_stom,      hist_id_stom_IPCC,               &
616         index,          lalo,              neighbours,   resolution,        &
617         contfrac,       totfrac_nobio,     clay,         t2m,               &
618         lai,            veget,             veget_max,                       &
619         co2_flux,       fco2_lu,           deadleaf_cover,  assim_param, temp_growth )
620
621    IMPLICIT NONE
622    !! 0. Variable and parameter declaration
623    !! 0.1 Input variables
624    INTEGER(i_std),INTENT(in)                       :: kjit              !! Time step number (unitless)
625    INTEGER(i_std),INTENT(in)                       :: kjpij             !! Total size of the un-compressed grid (unitless)
626    INTEGER(i_std),INTENT(in)                       :: kjpindex          !! Domain size - terrestrial pixels only (unitless)
627    INTEGER(i_std),INTENT(in)                       :: rest_id_stom      !! STOMATE's _Restart_ file identifier (unitless)
628    INTEGER(i_std),INTENT(in)                       :: hist_id_stom      !! STOMATE's _history_ file identifier (unitless)
629    INTEGER(i_std),INTENT(in)                       :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier(unitless)
630    INTEGER(i_std),DIMENSION(kjpindex),INTENT(in)   :: index             !! The indices of the terrestrial pixels only (unitless)
631    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in)    :: lalo              !! Geographical coordinates (latitude,longitude) for pixels (degrees)
632    INTEGER(i_std),DIMENSION(kjpindex,NbNeighb),INTENT(in) :: neighbours !! Neighoring grid points if land for the DGVM (unitless)
633    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in)    :: resolution        !! Size in x an y of the grid (m) - surface area of the gridbox
634    REAL(r_std),DIMENSION (kjpindex), INTENT (in)   :: contfrac          !! Fraction of continent in the grid cell (unitless)
635    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: totfrac_nobio     !! Fraction of grid cell covered by lakes, land ice, cities, ... (unitless)
636    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
637    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: t2m               !! 2 m air temperature (K)
638    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: lai               !! Leaf area inex @tex $(m^2 m^{-2})$ @endtex
639    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: veget             !! Fraction of vegetation type including
640                                                                         !! non-biological fraction (unitless)
641    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: veget_max         !! Maximum fraction of vegetation type including
642                                                                         !! non-biological fraction (unitless)
643
644    !! 0.2 Output variables
645    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: co2_flux          !! CO2 flux between atmosphere and biosphere
646    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: fco2_lu           !! CO2 flux between atmosphere and biosphere from land-use (without forest management) 
647    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: deadleaf_cover    !! Fraction of soil covered by dead leaves (unitless)
648    REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(out) :: assim_param !! min+max+opt temperatures (K) & vmax for photosynthesis 
649                                                                         !! @tex $(\mu mol m^{-2}s^{-1})$ @endtex 
650    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: temp_growth       !! Growth temperature (°C) 
651                                                                         !! Is equal to t2m_month
652    !! 0.3 Local variables
653    REAL(r_std)                                   :: dt_days_read             !! STOMATE time step read in restart file (days)
654    INTEGER(i_std)                                :: l,k,ji, jv, i, j, m      !! indices   
655    REAL(r_std),PARAMETER                         :: max_dt_days = 5.         !! Maximum STOMATE time step (days)
656    REAL(r_std),DIMENSION(kjpindex,nvm)           :: rprof                    !! Coefficient of the exponential functions that
657                                                                              !! relates root density to soil depth (unitless)
658    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_daily_x              !! "Daily" gpp for teststomate 
659                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex
660    REAL(r_std),DIMENSION(kjpindex,nvm)           :: veget_cov                !! Fractional coverage: actually share of the pixel
661                                                                              !! covered by a PFT (fraction of ground area),
662                                                                              !! taking into account LAI ??(= grid scale fpc)??
663    INTEGER(i_std)                                :: ier                      !! Check errors in netcdf call (unitless)
664
665    INTEGER(i_std)                                :: max_totsize              !! Memory management - maximum memory size (Mb)
666    INTEGER(i_std)                                :: totsize_1step            !! Memory management - memory required to store one
667                                                                              !! time step on one processor (Mb)
668    INTEGER(i_std)                                :: totsize_tmp              !! Memory management - memory required to store one
669                                                                              !! time step on all processors(Mb)
670    INTEGER(i_std)                                :: vid                      !! Variable identifer of netCDF (unitless)
671    INTEGER(i_std)                                :: nneigh                   !! Number of neighbouring pixels
672    INTEGER(i_std)                                :: direct                   !!
673    INTEGER(i_std),DIMENSION(ndm)                 :: d_id                     !!
674
675
676!_ ================================================================================================================================
677   
678    !! 1. Initialize variable
679    !! Update flag
680    l_first_stomate = .FALSE.
681
682    !! 1.1 Store current time step in a common variable
683    itime = kjit
684   
685    !![DISPENSABLE] 1.2 Copy the depth of the different soil layers from diaglev specified in slow_proc
686   
687    !! 1.3 PFT rooting depth across pixels, humescte is pre-defined
688    ! (constantes_veg.f90). It is defined as the coefficient of an exponential
689    ! function relating root density to depth
690    DO j=1,nvm
691       rprof(:,j) = 1./humcste(j)
692    ENDDO
693   
694    !! 1.4.0 Parameters for spinup
695    !
696    eps_carbon = 0.01
697    !Config Key   = EPS_CARBON
698    !Config Desc  = Allowed error on carbon stock
699    !Config If    = SPINUP_ANALYTIC
700    !Config Def   = 0.01
701    !Config Help  =
702    !Config Units = [%]   
703    CALL getin_p('EPS_CARBON',eps_carbon)       
704   
705   
706    !Config Key   = SPINUP_PERIOD
707    !Config Desc  = Period to calulcate equilibrium during spinup analytic
708    !Config If    = SPINUP_ANALYTIC
709    !Config Def   = -1
710    !Config Help  = Period corresponds in most cases to the number of years of forcing data used in the spinup.
711    !Config Units = [years]   
712    spinup_period = -1
713    CALL getin_p('SPINUP_PERIOD',spinup_period)       
714   
715    ! Check spinup_period values.
716    ! For periods uptil 6 years, to obtain equilibrium, a bigger period have to be used
717    ! and therefore spinup_period is adjusted to 10 years.
718    IF (spinup_analytic) THEN
719       IF (spinup_period <= 0) THEN
720          WRITE(numout,*) 'Error in parameter spinup_period. This parameter must be > 0 : spinup_period=',spinup_period
721          CALL ipslerr_p (3,'stomate_initialize', &
722               'Parameter spinup_period must be set to a positive integer.', &
723               'Set this parameter to the number of years of forcing data used for the spinup.', &
724               '')
725       ELSE IF (spinup_period <= 6) THEN
726          ! Adjust to bigger period. The period must be a multiple of the original period.
727          WRITE(numout,*) 'Initial spinup_period =',spinup_period,' will be adjusted.'
728          spinup_period = spinup_period*(INT(9/spinup_period)+1)
729       END IF
730       WRITE(numout,*) 'Spinup analytic is activated using eps_carbon=',eps_carbon, ' and spinup_period=',spinup_period
731    END IF
732   
733    !! 1.4.0 Initialization of PFT specific parameters
734    ! Initialization of PFT specific parameters that have no value
735    ! for the bare soil PFT i.e. fire resistance, flamability, maximum lai,
736    ! settings for growing degree days (GDD), settings for senescence,
737    ! respiration coefficients, photosynthesis, etc.
738    ! [DISPENSABLE]
739   
740    !! 1.4.1 Allocate memory for all variables in stomate
741    ! Allocate memory for all variables in stomate, build new index
742    ! tables accounting for the PFTs, read and check flags and set file
743    ! identifier for restart and history files.
744    CALL stomate_init (kjpij, kjpindex, index, lalo, &
745         rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
746   
747    !! 1.4.2 Initialization of PFT specific parameters
748    ! Initialization of PFT specific parameters i.e. sla from leaf life,
749    ! sapling characteristics (biomass), migration speed, critical diameter,
750    ! coldest tolerable temperature, critical values for phenology, maximum
751    ! life time of leaves, respiration coefficients and photosynthesis.
752    ! The subroutine also communicates settings read by stomate_constant_init.
753    CALL data (kjpindex, lalo)
754   
755    !! 1.4.3 Initial conditions
756   
757    !! 1.4.3.1 Read initial values for STOMATE's variables from the _restart_ file
758    ! ??Shouldn't this be included in stomate_init?? Looks like an initialization!
759    co2_flux(:,:) = zero
760    fco2_lu(:) = zero
761   
762    ! Get values from _restart_ file. Note that only ::kjpindex, ::index, ::lalo
763    ! and ::resolution are input variables, all others are output variables.
764    CALL readstart &
765         (kjpindex, index, lalo, resolution, t2m, &
766         dt_days_read, date, &
767         ind, adapted, regenerate, &
768         humrel_daily, gdd_init_date, litterhum_daily, &
769         t2m_daily, t2m_min_daily, tsurf_daily, tsoil_daily, &
770         soilhum_daily, precip_daily, &
771         gpp_daily, npp_daily, turnover_daily, &
772         humrel_month, humrel_week, &
773         t2m_longterm, tau_longterm, t2m_month, t2m_week, &
774         tsoil_month, soilhum_month, fireindex, firelitter, &
775         maxhumrel_lastyear, maxhumrel_thisyear, &
776         minhumrel_lastyear, minhumrel_thisyear, &
777         maxgppweek_lastyear, maxgppweek_thisyear, &
778         gdd0_lastyear, gdd0_thisyear, &
779         precip_lastyear, precip_thisyear, &
780         gdd_m5_dormance,  gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
781         PFTpresent, npp_longterm, lm_lastyearmax, lm_thisyearmax, &
782         maxfpc_lastyear, maxfpc_thisyear, &
783         turnover_longterm, gpp_week, biomass, resp_maint_part, &
784         leaf_age, leaf_frac, &
785         senescence, when_growthinit, age, &
786         resp_hetero_d, resp_maint_d, resp_growth_d, co2_fire, co2_to_bm_dgvm, &
787         veget_lastlight, everywhere, need_adjacent, RIP_time, &
788         time_hum_min, hum_min_dormance, &
789         litterpart, litter, dead_leaves, &
790         carbon, lignin_struc,turnover_time,&
791         prod10,prod100,flux10, flux100, &
792         convflux, cflux_prod10, cflux_prod100, bm_to_litter, carb_mass_total, &
793         Tseason, Tseason_length, Tseason_tmp, &
794         Tmin_spring_time, begin_leaves, onset_date, &
795         global_years, ok_equilibrium, nbp_accu, nbp_flux, &
796         MatrixV, VectorU, previous_stock, current_stock, assim_param, &
797!gmjc
798        & wshtotsum, sr_ugb, sla_calc, nb_ani, grazed_frac, &
799        & import_yield, t2m_14, litter_not_avail, nb_grazingdays, &
800        & after_snow, after_wet, wet1day, wet2day)
801!end gmjc
802   
803    !! 1.4.5 Check time step
804       
805    !! 1.4.5.1 Allow STOMATE's time step to change although this is dangerous
806    IF (dt_days /= dt_days_read) THEN
807       WRITE(numout,*) 'slow_processes: STOMATE time step changes:', &
808            & dt_days_read,' -> ',dt_days
809    ENDIF
810   
811    !! 1.4.5.2 Time step has to be a multiple of a full day
812    IF ( ( dt_days-REAL(NINT(dt_days),r_std) ) > min_stomate ) THEN
813       WRITE(numout,*) 'slow_processes: STOMATE time step is not a mutiple of a full day:', &
814            & dt_days,' days.'
815       STOP
816    ENDIF
817   
818    !! 1.4.5.3 upper limit to STOMATE's time step
819    IF ( dt_days > max_dt_days ) THEN
820       WRITE(numout,*) 'slow_processes: STOMATE time step exceeds the maximum value:', &
821            & dt_days,' days > ', max_dt_days, ' days.' 
822       STOP
823    ENDIF
824   
825    !! 1.4.5.4 STOMATE time step must not be less than the forcing time step
826    IF ( dt_sechiba > dt_days*one_day ) THEN
827       WRITE(numout,*) &
828            & 'slow_processes: STOMATE time step ::dt_days smaller than forcing time step ::dt_sechiba'
829       STOP
830    ENDIF
831   
832    !! 1.4.5.6 Final message on time step
833    WRITE(numout,*) 'Slow_processes, STOMATE time step (days): ', dt_days
834   
835    !! 1.4.6 Write forcing file for teststomate
836    IF (ok_co2 .AND. allow_forcing_write) THEN
837       
838       !Config Key   = STOMATE_FORCING_NAME
839       !Config Desc  = Name of STOMATE's forcing file
840       !Config If    = OK_STOMATE
841       !Config Def   = NONE
842       !Config Help  = Name that will be given
843       !Config         to STOMATE's offline forcing file
844       !Config         Compatible with Nicolas Viovy's driver
845       !Config Units = [FILE]
846       forcing_name = stomate_forcing_name
847       CALL getin_p('STOMATE_FORCING_NAME',forcing_name)
848       
849       IF ( TRIM(forcing_name) /= 'NONE' ) THEN
850         
851          !! 1.4.6.1 Calculate steps that can be stored in memory
852          ! Action for the root processor only (parallel computing) 
853          IF (is_root_prc) CALL SYSTEM ('rm -f '//TRIM(forcing_name))
854          WRITE(numout,*) 'writing a forcing file for STOMATE.'
855         
856          !Config Key   = STOMATE_FORCING_MEMSIZE
857          !Config Desc  = Size of STOMATE forcing data in memory
858          !Config If    = OK_STOMATE
859          !Config Def   = 50
860          !Config Help  = This variable determines how many
861          !Config         forcing states will be kept in memory.
862          !Config         Must be a compromise between memory
863          !Config         use and frequeny of disk access.
864          !Config Units = [MegaBytes]
865          max_totsize = 50
866          CALL getin_p('STOMATE_FORCING_MEMSIZE', max_totsize)     
867          max_totsize = max_totsize*1000000
868         
869          totsize_1step = &
870               SIZE(clay)*KIND(clay) &
871               +SIZE(humrel_daily)*KIND(humrel_daily) &
872               +SIZE(litterhum_daily)*KIND(litterhum_daily) &
873               +SIZE(t2m_daily)*KIND(t2m_daily) &
874               +SIZE(t2m_min_daily)*KIND(t2m_min_daily) &
875               +SIZE(tsurf_daily)*KIND(tsurf_daily) &
876               +SIZE(tsoil_daily)*KIND(tsoil_daily) &
877               +SIZE(soilhum_daily)*KIND(soilhum_daily) &
878               +SIZE(precip_daily)*KIND(precip_daily) &
879               +SIZE(gpp_daily_x)*KIND(gpp_daily_x) &
880               +SIZE(veget)*KIND(veget) &
881               +SIZE(veget_max)*KIND(veget_max) &
882               +SIZE(lai)*KIND(lai)
883         
884          ! Totsize_1step is the size on a single processor, sum
885          ! all processors and send to all processors
886          CALL reduce_sum(totsize_1step,totsize_tmp)
887          CALL bcast(totsize_tmp)
888          totsize_1step=totsize_tmp
889         
890          ! Total number of forcing steps
891          nsft = INT(one_year/(dt_stomate/one_day))
892         
893          ! Number of forcing steps in memory
894          nsfm = MIN(nsft, &
895               MAX(1,NINT( REAL(max_totsize,r_std) &
896               /REAL(totsize_1step,r_std))))
897           
898             
899          !! 1.6.4.2 Allocate memory for variables containing forcing data 
900          ! and initialize variables (set to zero).
901          CALL init_forcing (kjpindex,nsfm,nsft)
902         
903          ! Indexing for writing forcing file
904          isf(:) = (/ (i,i=1,nsfm) /)
905          nf_written(:) = .FALSE.
906          nf_cumul(:) = 0
907          iisf = 0
908         
909          !! 1.6.4.3 Create netcdf file
910          ! Create, define and populate a netcdf file containing the forcing data.
911          ! For the root processor only (parallel computing). NF90_ are functions
912          ! from and external library. 
913          IF (is_root_prc) THEN
914             
915             ! Create new netCDF dataset
916             ier = NF90_CREATE (TRIM(forcing_name),NF90_SHARE,forcing_id)
917             
918             ! Add variable attribute
919             ! Note ::iim_g and ::jjm_g are dimensions of the global field and
920             ! ::nbp_glo is the number of global continental points
921             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL,'dt_sechiba',dt_sechiba)
922             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL,'dt_stomate',dt_stomate)
923             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
924                  'nsft',REAL(nsft,r_std))
925             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
926                  'kjpij',REAL(iim_g*jjm_g,r_std))
927             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
928                  'kjpindex',REAL(nbp_glo,r_std))
929             
930             ! Add new dimension
931             ier = NF90_DEF_DIM (forcing_id,'points',nbp_glo,d_id(1))
932             ier = NF90_DEF_DIM (forcing_id,'layers',nbdl,d_id(2))
933             ier = NF90_DEF_DIM (forcing_id,'pft',nvm,d_id(3))
934             direct=2
935             ier = NF90_DEF_DIM (forcing_id,'direction',direct,d_id(4))
936             nneigh=8
937             ier = NF90_DEF_DIM (forcing_id,'nneigh',nneigh,d_id(5))
938             ier = NF90_DEF_DIM (forcing_id,'time',NF90_UNLIMITED,d_id(6))
939             ier = NF90_DEF_DIM (forcing_id,'nbparts',nparts,d_id(7))
940             
941             ! Add new variable
942             ier = NF90_DEF_VAR (forcing_id,'points',    r_typ,d_id(1),vid)
943             ier = NF90_DEF_VAR (forcing_id,'layers',    r_typ,d_id(2),vid)
944             ier = NF90_DEF_VAR (forcing_id,'pft',       r_typ,d_id(3),vid)
945             ier = NF90_DEF_VAR (forcing_id,'direction', r_typ,d_id(4),vid)
946             ier = NF90_DEF_VAR (forcing_id,'nneigh',    r_typ,d_id(5),vid)
947             ier = NF90_DEF_VAR (forcing_id,'time',      r_typ,d_id(6),vid)
948             ier = NF90_DEF_VAR (forcing_id,'nbparts',   r_typ,d_id(7),vid)
949             ier = NF90_DEF_VAR (forcing_id,'index',     r_typ,d_id(1),vid)
950             ier = NF90_DEF_VAR (forcing_id,'contfrac',  r_typ,d_id(1),vid) 
951             ier = NF90_DEF_VAR (forcing_id,'lalo', &
952                  r_typ,(/ d_id(1),d_id(4) /),vid)
953             ier = NF90_DEF_VAR (forcing_id,'neighbours', &
954                  r_typ,(/ d_id(1),d_id(5) /),vid)
955             ier = NF90_DEF_VAR (forcing_id,'resolution', &
956                  r_typ,(/ d_id(1),d_id(4) /),vid)
957             ier = NF90_DEF_VAR (forcing_id,'clay', &
958                  r_typ,(/ d_id(1),d_id(6) /),vid)
959             ier = NF90_DEF_VAR (forcing_id,'humrel', &
960                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
961             ier = NF90_DEF_VAR (forcing_id,'litterhum', &
962                  r_typ,(/ d_id(1),d_id(6) /),vid)
963             ier = NF90_DEF_VAR (forcing_id,'t2m', &
964                  r_typ,(/ d_id(1),d_id(6) /),vid)
965             ier = NF90_DEF_VAR (forcing_id,'t2m_min', &
966                  r_typ,(/ d_id(1),d_id(6) /),vid)
967             ier = NF90_DEF_VAR (forcing_id,'tsurf', &
968                  r_typ,(/ d_id(1),d_id(6) /),vid)
969             ier = NF90_DEF_VAR (forcing_id,'tsoil', &
970                  r_typ,(/ d_id(1),d_id(2),d_id(6) /),vid)
971             ier = NF90_DEF_VAR (forcing_id,'soilhum', &
972                  r_typ,(/ d_id(1),d_id(2),d_id(6) /),vid)
973             ier = NF90_DEF_VAR (forcing_id,'precip', &
974                  r_typ,(/ d_id(1),d_id(6) /),vid)
975             ier = NF90_DEF_VAR (forcing_id,'gpp', &
976                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
977             ier = NF90_DEF_VAR (forcing_id,'veget', &
978                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
979             ier = NF90_DEF_VAR (forcing_id,'veget_max', &
980                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
981             ier = NF90_DEF_VAR (forcing_id,'lai', &
982                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
983             ier = NF90_ENDDEF (forcing_id)
984             
985             ! Given the name of a varaible, nf90_inq_varid finds the variable
986             ! ID (::vid). Put data value(s) into variable ::vid
987             ier = NF90_INQ_VARID (forcing_id,'points',vid)
988             ier = NF90_PUT_VAR (forcing_id,vid, &
989                  (/(REAL(i,r_std),i=1,nbp_glo) /))
990             ier = NF90_INQ_VARID (forcing_id,'layers',vid)
991             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nbdl)/))
992             ier = NF90_INQ_VARID (forcing_id,'pft',vid)
993             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nvm)/))
994             ier = NF90_INQ_VARID (forcing_id,'direction',vid)
995             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,2)/))
996             ier = NF90_INQ_VARID (forcing_id,'nneigh',vid)
997             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,8)/))
998             ier = NF90_INQ_VARID (forcing_id,'time',vid)
999             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nsft)/))
1000             ier = NF90_INQ_VARID (forcing_id,'nbparts',vid)
1001             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nparts)/))
1002             ier = NF90_INQ_VARID (forcing_id,'index',vid) 
1003             ier = NF90_PUT_VAR (forcing_id,vid,REAL(index_g,r_std))
1004             ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
1005             ier = NF90_PUT_VAR (forcing_id,vid,REAL(contfrac_g,r_std))
1006             ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
1007             ier = NF90_PUT_VAR (forcing_id,vid,lalo_g)
1008             !ym attention a neighbours, a modifier plus tard     
1009             ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
1010             ier = NF90_PUT_VAR (forcing_id,vid,REAL(neighbours_g,r_std))
1011             ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
1012             ier = NF90_PUT_VAR (forcing_id,vid,resolution_g)
1013          ENDIF ! is_root_prc
1014       ENDIF ! (forcing_name) /= 'NONE'
1015    ENDIF ! ok_co2 =.TRUE.
1016   
1017    !! 1.4.7 write forcing file for forcesoil
1018    !! 1.4.7.1 Initialize
1019    !Config Key   = STOMATE_CFORCING_NAME
1020    !Config Desc  = Name of STOMATE's carbon forcing file
1021    !Config If    = OK_STOMATE
1022    !Config Def   = NONE
1023    !Config Help  = Name that will be given to STOMATE's carbon
1024    !Config         offline forcing file
1025    !Config         Compatible with Nicolas Viovy's driver
1026    !Config Units = [FILE]
1027    Cforcing_name = stomate_Cforcing_name
1028    CALL getin_p('STOMATE_CFORCING_NAME',Cforcing_name)
1029   
1030    IF ( TRIM(Cforcing_name) /= 'NONE' ) THEN
1031       
1032       ! Time step of forcesoil
1033       !Config Key   = FORCESOIL_STEP_PER_YEAR
1034       !Config Desc  = Number of time steps per year for carbon spinup.
1035       !Config If    = OK_STOMATE
1036       !Config Def   = 365
1037       !Config Help  = Number of time steps per year for carbon spinup.
1038       !Config Units = [days, months, year]
1039       nparan = 365
1040       CALL getin_p('FORCESOIL_STEP_PER_YEAR', nparan)
1041       
1042       ! Correct if setting is out of bounds
1043       IF ( nparan < 1 ) nparan = 1
1044       
1045       !Config Key   = FORCESOIL_NB_YEAR
1046       !Config Desc  = Number of years saved for carbon spinup.
1047       !Config If    = OK_STOMATE
1048       !Config Def   = 1
1049       !Config Help  = Number of years saved for carbon spinup. If internal parameter cumul_Cforcing is TRUE in stomate.f90
1050       !Config         Then this parameter is forced to one.
1051       !Config Units = [years]
1052       CALL getin_p('FORCESOIL_NB_YEAR', nbyear)
1053       
1054       ! Set ::nbyear to 1. if ::cumul_Cforcing=.TRUE.
1055       IF ( cumul_Cforcing ) THEN
1056          CALL ipslerr_p (1,'stomate', &
1057               'Internal parameter cumul_Cforcing is TRUE in stomate.f90', &
1058               'Parameter FORCESOIL_NB_YEAR is therefore forced to 1.', &
1059               '::nbyear is thus set to 1.')
1060          nbyear=1
1061       ENDIF
1062       
1063       ! Make use of ::nparan to calculate ::dt_forcesoil
1064       dt_forcesoil = zero
1065       nparan = nparan+1
1066       DO WHILE ( dt_forcesoil < dt_stomate/one_day )
1067          nparan = nparan-1
1068          IF ( nparan < 1 ) THEN
1069             STOP 'Problem with number of soil forcing time steps ::nparan < 1.'
1070          ENDIF
1071          dt_forcesoil = one_year/REAL(nparan,r_std)
1072       ENDDO
1073       IF ( nparan > nparanmax ) THEN
1074          STOP 'Problem with number of soil forcing time steps ::nparan > ::nparanmax'
1075       ENDIF
1076       WRITE(numout,*) 'Time step of soil forcing (d): ',dt_forcesoil
1077       
1078       ! Allocate memory for the forcing variables of soil dynamics
1079       ALLOCATE( nforce(nparan*nbyear))
1080       nforce(:) = 0
1081       ALLOCATE(control_moist(kjpindex,nlevs,nparan*nbyear))
1082       ALLOCATE(npp_equil(kjpindex,nparan*nbyear))
1083       ALLOCATE(npp_tot(kjpindex))
1084       ALLOCATE(control_temp(kjpindex,nlevs,nparan*nbyear))
1085       ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear)) 
1086       
1087       ! Initialize variables, set to zero
1088       control_moist(:,:,:) = zero
1089       npp_equil(:,:) = zero
1090       npp_tot(:) = zero
1091       control_temp(:,:,:) = zero
1092       soilcarbon_input(:,:,:,:) = zero
1093       
1094    ENDIF ! Cforcing_name) /= 'NONE'
1095   
1096    !! 1.4.8 Calculate STOMATE's vegetation fractions from veget, veget_max
1097    DO j=1,nvm
1098       WHERE ((1.-totfrac_nobio(:)) > min_sechiba)       
1099          ! Pixels with vegetation
1100          veget_cov(:,j) = veget(:,j)/( 1.-totfrac_nobio(:) )
1101          veget_cov_max(:,j) = veget_max(:,j)/( 1.-totfrac_nobio(:) )
1102       ELSEWHERE
1103          ! Pixels without vegetation
1104          veget_cov(:,j) = zero
1105          veget_cov_max(:,j) = zero
1106       ENDWHERE
1107    ENDDO ! Loop over PFTs
1108
1109    !! 1.4.9 Initialize non-zero variables
1110    CALL stomate_var_init &
1111         (kjpindex, veget_cov_max, leaf_age, leaf_frac, &
1112         dead_leaves, &
1113         veget, lai, deadleaf_cover, assim_param, &
1114!gmjc
1115           &          N_limfert)
1116!end gmjc
1117   
1118    ! Initialize land cover change variable
1119    ! ??Should be integrated in the subroutine??
1120    harvest_above(:) = zero
1121   
1122    ! Initialize temp_growth
1123    temp_growth(:)=t2m_month(:)-tp_00 
1124
1125     
1126  END SUBROUTINE stomate_initialize
1127 
1128
1129!! ================================================================================================================================
1130!! SUBROUTINE   : stomate_main
1131!!
1132!>\BRIEF        Manages variable initialisation, reading and writing forcing
1133!! files, aggregating data at stomate's time step (dt_stomate), aggregating data
1134!! at longer time scale (i.e. for phenology) and uses these forcing to calculate
1135!! CO2 fluxes (NPP and respirations) and C-pools (litter, soil, biomass, ...)
1136!!
1137!! DESCRIPTION  : The subroutine manages
1138!! divers tasks:
1139!! (1) Initializing all variables of stomate (first call)
1140!! (2) Reading and writing forcing data (last call)
1141!! (3) Adding CO2 fluxes to the IPCC history files
1142!! (4) Converting the time steps of variables to maintain consistency between
1143!! sechiba and stomate
1144!! (5) Use these variables to call stomate_lpj, maint_respiration, littercalc,
1145!! soilcarbon. The called subroutines handle: climate constraints
1146!! for PFTs, PFT dynamics, Phenology, Allocation, NPP (based on GPP and
1147!! authothropic respiration), fire, mortality, vmax, assimilation temperatures,
1148!! all turnover processes, light competition, sapling establishment, lai, 
1149!! land cover change and litter and soil dynamics.
1150!! (6) Use the spin-up method developed by Lardy (2011)(only if SPINUP_ANALYTIC
1151!! is set to TRUE).
1152!!
1153!! RECENT CHANGE(S) : None
1154!!
1155!! MAIN OUTPUT VARIABLE(S): deadleaf_cover, assim_param, lai, height, veget,
1156!! veget_max, resp_maint,
1157!! resp_hetero,resp_growth, co2_flux, fco2_lu.
1158!!
1159!! REFERENCES   :
1160!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
1161!! Environmental Modelling & Software (2011), doi:10.1016|j.envsoft.2011.05.016
1162!!
1163!! FLOWCHART    :
1164!! \latexonly
1165!! \includegraphics[scale=0.5]{stomatemainflow.png}
1166!! \endlatexonly
1167!! \n
1168!_ ================================================================================================================================
1169 
1170SUBROUTINE stomate_main &
1171       & (kjit, kjpij, kjpindex, &
1172       &  index, lalo, neighbours, resolution, contfrac, totfrac_nobio, clay, &
1173       &  t2m, t2m_min, temp_sol, stempdiag, &
1174       &  humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, &
1175       &  gpp, deadleaf_cover, assim_param, &
1176       &  lai, frac_age, height, veget, veget_max, &
1177       &  veget_max_new, totfrac_nobio_new, &
1178       &  hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
1179       &  co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth,temp_growth,&
1180       &  EndOfYear, &
1181!JCADD top 5 layer grassland soil moisture for grazing
1182            tmc_topgrass,fc_grazing, snow)
1183!ENDJCADD
1184   
1185    IMPLICIT NONE
1186
1187   
1188  !! 0. Variable and parameter declaration
1189
1190    !! 0.1 Input variables
1191
1192    INTEGER(i_std),INTENT(in)                       :: kjit              !! Time step number (unitless)
1193    INTEGER(i_std),INTENT(in)                       :: kjpindex          !! Domain size - terrestrial pixels only (unitless)
1194    INTEGER(i_std),INTENT(in)                       :: kjpij             !! Total size of the un-compressed grid (unitless)
1195    INTEGER(i_std),INTENT(in)                       :: rest_id_stom      !! STOMATE's _Restart_ file identifier (unitless)
1196    INTEGER(i_std),INTENT(in)                       :: hist_id_stom      !! STOMATE's _history_ file identifier (unitless)
1197    INTEGER(i_std),INTENT(in)                       :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier
1198                                                                         !! (unitless)
1199    INTEGER(i_std),DIMENSION(kjpindex),INTENT(in)   :: index             !! Indices of the pixels on the map. Stomate uses a
1200                                                                         !! reduced grid excluding oceans. ::index contains
1201                                                                         !! the indices of the terrestrial pixels only
1202                                                                         !! (unitless)
1203    INTEGER(i_std),DIMENSION(kjpindex,NbNeighb),INTENT(in) :: neighbours !! Neighoring grid points if land for the DGVM
1204                                                                         !! (unitless)
1205    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in)    :: lalo              !! Geographical coordinates (latitude,longitude)
1206                                                                         !! for pixels (degrees)
1207    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in)    :: resolution        !! Size in x an y of the grid (m) - surface area of
1208                                                                         !! the gridbox
1209    REAL(r_std),DIMENSION (kjpindex), INTENT (in)   :: contfrac          !! Fraction of continent in the grid cell (unitless)
1210    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: totfrac_nobio     !! Fraction of grid cell covered by lakes, land
1211                                                                         !! ice, cities, ... (unitless)
1212    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
1213    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: humrel            !! Relative humidity ("moisture availability")
1214                                                                         !! (0-1, unitless)
1215    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: t2m               !! 2 m air temperature (K)
1216    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: t2m_min           !! Minimum 2 m air temp. during forcing time step
1217                                                                         !! (K)
1218    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: temp_sol          !! Surface temperature (K)
1219    REAL(r_std),DIMENSION(kjpindex,nbdl),INTENT(in) :: stempdiag         !! Soil temperature (K)
1220    REAL(r_std),DIMENSION(kjpindex,nbdl),INTENT(in) :: shumdiag          !! Relative soil moisture (0-1, unitless)
1221    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: litterhumdiag     !! Litter humidity (0-1, unitless)
1222    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: precip_rain       !! Rain precipitation 
1223                                                                         !! @tex $(mm dt_stomate^{-1})$ @endtex
1224    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: precip_snow       !! Snow precipitation 
1225                                                                         !! @tex $(mm dt_stomate^{-1})$ @endtex
1226    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: gpp               !! GPP of total ground area 
1227                                                                         !! @tex $(gC m^{-2} time step^{-1})$ @endtex
1228                                                                         !! Calculated in sechiba, account for vegetation
1229                                                                         !! cover and effective time step to obtain ::gpp_d
1230    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: veget_max_new     !! New "maximal" coverage fraction of a PFT: only if
1231                                                                         !! vegetation is updated in slowproc
1232    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: totfrac_nobio_new !! New fraction of nobio per gridcell
1233    INTEGER(i_std),INTENT(in)                       :: hist_id           !! ?? [DISPENSABLE] SECHIBA's _history_ file
1234                                                                         !! identifier
1235    INTEGER(i_std),INTENT(in)                       :: hist2_id          !! ?? [DISPENSABLE] SECHIBA's _history_ file 2
1236                                                                         !! identifier
1237    LOGICAL, INTENT(in)                             :: EndOfYear       !! Flag set to true for the first sechiba time step on the year.
1238!JCADD top 5 layer grassland soil moisture for grazing
1239    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: tmc_topgrass
1240    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: fc_grazing
1241    REAL(r_std),DIMENSION(kjpindex),INTENT(in)         :: snow ! snow mass
1242!ENDJCADD
1243    !! 0.2 Output variables
1244
1245    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: co2_flux          !! CO2 flux between atmosphere and biosphere per
1246                                                                         !! average ground area
1247                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex 
1248                                                                         !! [??CHECK] sign convention?
1249    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: fco2_lu           !! CO2 flux between atmosphere and biosphere from
1250                                                                         !! land-use (without forest management) 
1251                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex 
1252                                                                         !! [??CHECK] sign convention?
1253    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_maint        !! Maitenance component of autotrophic respiration in
1254                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
1255    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_growth       !! Growth component of autotrophic respiration in
1256                                                                         !! @tex ($gC m^{-2} dt_stomate^{-1}$) @endtex
1257    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_hetero       !! Heterotrophic respiration in 
1258                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex 
1259    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: temp_growth       !! Growth temperature (°C) 
1260                                                                         !! Is equal to t2m_month
1261
1262    !! 0.3 Modified
1263   
1264    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout)       :: lai            !! Leaf area inex @tex $(m^2 m^{-2})$ @endtex
1265    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)          :: veget          !! Fraction of vegetation type including
1266                                                                              !! non-biological fraction (unitless)
1267    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout)       :: veget_max      !! Maximum fraction of vegetation type including
1268                                                                              !! non-biological fraction (unitless)
1269    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout)       :: height         !! Height of vegetation (m)
1270    REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(inout) :: assim_param    !! min+max+opt temperatures (K) & vmax for
1271                                                                              !! photosynthesis 
1272                                                                              !! @tex $(\mu mol m^{-2}s^{-1})$ @endtex 
1273    REAL(r_std),DIMENSION(kjpindex),INTENT(inout)           :: deadleaf_cover !! Fraction of soil covered by dead leaves
1274                                                                              !! (unitless)
1275    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages),INTENT(inout):: frac_age    !! Age efficacity from STOMATE
1276
1277    !! 0.4 local variables
1278   
1279    REAL(r_std)                                   :: dt_days_read             !! STOMATE time step read in restart file (days)
1280    INTEGER(i_std)                                :: l,k,ji, jv, i, j, m      !! indices   
1281    REAL(r_std),PARAMETER                         :: max_dt_days = 5.         !! Maximum STOMATE time step (days)
1282    REAL(r_std)                                   :: hist_days                !! Writing frequency for history file (days)
1283    REAL(r_std),DIMENSION(0:nbdl)                 :: z_soil                   !! Variable to store depth of the different soil
1284                                                                              !! layers (m)
1285    REAL(r_std),DIMENSION(kjpindex,nvm)           :: rprof                    !! Coefficient of the exponential functions that
1286                                                                              !! relates root density to soil depth (unitless)
1287    REAL(r_std),DIMENSION(kjpindex)               :: cvegtot                  !! Total "vegetation" cover (unitless)
1288    REAL(r_std),DIMENSION(kjpindex)               :: precip                   !! Total liquid and solid precipitation 
1289                                                                              !! @tex $(??mm dt_stomate^{-1})$ @endtex
1290    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_d                    !! Gross primary productivity per ground area
1291                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex 
1292    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_daily_x              !! "Daily" gpp for teststomate 
1293                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex
1294    REAL(r_std),DIMENSION(kjpindex,nvm)           :: resp_hetero_litter       !! Litter heterotrophic respiration per ground area
1295                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex 
1296                                                                              !! ??Same variable is also used to
1297                                                                              !! store heterotrophic respiration per ground area
1298                                                                              !! over ::dt_sechiba??
1299    REAL(r_std),DIMENSION(kjpindex,nvm)           :: resp_hetero_soil         !! soil heterotrophic respiration 
1300                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex
1301    REAL(r_std),DIMENSION(kjpindex,nvm)           :: veget_cov                !! Fractional coverage: actually share of the pixel
1302                                                                              !! covered by a PFT (fraction of ground area),
1303                                                                              !! taking into account LAI ??(= grid scale fpc)??
1304    REAL(r_std),DIMENSION(kjpindex,nvm)           :: veget_cov_max_new        !! New value for maximal fractional coverage (unitless)
1305    REAL(r_std),DIMENSION(kjpindex,nvm)           :: vcmax                    !! Maximum rate of carboxylation
1306                                                                              !! @tex $(\mumol m^{-2} s^{-1})$ @endtex
1307    REAL(r_std),DIMENSION(kjpindex,nlevs)         :: control_moist_inst       !! Moisture control of heterotrophic respiration
1308                                                                              !! (0-1, unitless)
1309    REAL(r_std),DIMENSION(kjpindex,nlevs)         :: control_temp_inst        !! Temperature control of heterotrophic
1310                                                                              !! respiration, above and below (0-1, unitless)
1311    REAL(r_std),DIMENSION(kjpindex,ncarb,nvm)     :: soilcarbon_input_inst    !! Quantity of carbon going into carbon pools from
1312                                                                              !! litter decomposition
1313                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex
1314   
1315    INTEGER(i_std)                                :: ier                      !! Check errors in netcdf call (unitless)
1316    REAL(r_std)                                   :: sf_time                  !! Intermediate variable to calculate current time
1317                                                                              !! step
1318    INTEGER(i_std)                                :: max_totsize              !! Memory management - maximum memory size (Mb)
1319    INTEGER(i_std)                                :: totsize_1step            !! Memory management - memory required to store one
1320                                                                              !! time step on one processor (Mb)
1321    INTEGER(i_std)                                :: totsize_tmp              !! Memory management - memory required to store one
1322                                                                              !! time step on all processors(Mb)
1323    REAL(r_std)                                   :: xn                       !! How many times have we treated in this forcing
1324    REAL(r_std), DIMENSION(kjpindex)              :: vartmp                   !! Temporary variable
1325    INTEGER(i_std)                                :: vid                      !! Variable identifer of netCDF (unitless)
1326    INTEGER(i_std)                                :: nneigh                   !! Number of neighbouring pixels
1327    INTEGER(i_std)                                :: direct                   !! ??
1328    INTEGER(i_std),DIMENSION(ndm)                 :: d_id                     !! ??
1329    REAL(r_std)                                   :: net_co2_flux_monthly     !! ??[DISPENSABLE]
1330    REAL(r_std)                                   :: net_co2_flux_monthly_sum !! ??[DISPENSABLE]
1331    REAL(r_std),DIMENSION(nbp_glo)                :: clay_g                   !! Clay fraction of soil (0-1, unitless), parallel
1332                                                                              !! computing
1333    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:,:)    :: soilcarbon_input_g       !! Quantity of carbon going into carbon pools from
1334                                                                              !! litter decomposition 
1335                                                                              !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex, parallel
1336                                                                              !! computing
1337    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)      :: control_moist_g          !! Moisture control of heterotrophic respiration
1338                                                                              !! (0-1, unitless), parallel computing
1339    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)      :: control_temp_g           !! Temperature control of heterotrophic respiration
1340                                                                              !! (0-1, unitless), parallel computing
1341    REAL(r_std),ALLOCATABLE,DIMENSION(:,:)        :: npp_equil_g              !! Equilibrium NPP written to forcesoil
1342                                                                              !! @tex $(gC m^{-2} year^{-1})$ @endtex, parallel
1343                                                                              !! computing
1344
1345    REAL(r_std)                                   :: net_cflux_prod_monthly_sum    !! AR5 output?? gC m2 month-1 (one variable for
1346                                                                                   !! reduce_sum and one for bcast??), parallel
1347                                                                                   !! computing
1348    REAL(r_std)                                   :: net_cflux_prod_monthly_tot    !! AR5 output?? gC m2 month-1 (one variable for
1349                                                                                   !! reduce_sum and one for bcast??), parallel
1350                                                                                   !! computing
1351    REAL(r_std)                                   :: net_harvest_above_monthly_sum !! AR5 output?? gC m2 month-1 (one variable for
1352                                                                                   !! reduce_sum and one for bcast??), parallel
1353                                                                                   !! computing
1354    REAL(r_std)                                   :: net_harvest_above_monthly_tot !! AR5 output?? gC m2 month-1 (one variable for
1355                                                                                   !! reduce_sum and one for bcast??), parallel
1356                                                                                   !! computing
1357    REAL(r_std)                                   :: net_biosp_prod_monthly_sum    !! AR5 output?? gC m2 month-1 (one variable for
1358                                                                                   !! reduce_sum and one for bcast??), parallel
1359                                                                                   !! computing
1360    REAL(r_std)                                   :: net_biosp_prod_monthly_tot    !! AR5 output?? gC m2 month-1 (one variable for
1361                                                                                   !! reduce_sum and one for bcast??), parallel
1362                                                                                   !! computing
1363    REAL(r_std), DIMENSION(kjpindex,nvm,nbpools)  :: carbon_stock                  !! Array containing the carbon stock for each pool
1364                                                                                   !! used by ORCHIDEE
1365
1366!_ ================================================================================================================================
1367   
1368  !! 1. Initialize variables
1369
1370    !! 1.1 Store current time step in a common variable
1371    itime = kjit
1372   
1373    !![DISPENSABLE] 1.2 Copy the depth of the different soil layers from diaglev specified in slow_proc
1374                                                                                                                       
1375    !! 1.3 PFT rooting depth across pixels, humescte is pre-defined
1376    ! (constantes_veg.f90). It is defined as the coefficient of an exponential
1377    ! function relating root density to depth
1378    DO j=1,nvm
1379       rprof(:,j) = 1./humcste(j)
1380    ENDDO
1381   
1382    !! 1.4 Initialize first call
1383    ! Set growth respiration to zero
1384    resp_growth=zero
1385
1386    ! Check that initialization is done
1387    IF (l_first_stomate) CALL ipslerr_p(3,'stomate_main','Initialization not yet done.','','')
1388   
1389    IF (printlev >= 4) THEN
1390       WRITE(numout,*) 'stomate_main: date=',date,' ymds=', year, month, day, sec, ' itime=', itime, ' do_slow=',do_slow
1391    ENDIF
1392
1393!! 3. Special treatment for some input arrays.
1394   
1395    !! 3.1 Sum of liquid and solid precipitation
1396    precip(:) = ( precip_rain(:) + precip_snow(:) )*one_day/dt_sechiba
1397   
1398    !! 3.2 Calculate STOMATE's vegetation fractions from veget and veget_max
1399    DO j=1,nvm 
1400       WHERE ((1.-totfrac_nobio(:)) > min_sechiba)
1401          ! Pixels with vegetation
1402          veget_cov(:,j) = veget(:,j)/( 1.-totfrac_nobio(:) )
1403          veget_cov_max(:,j) = veget_max(:,j)/( 1.-totfrac_nobio(:) )
1404       ELSEWHERE
1405          ! Pixels without vegetation
1406          veget_cov(:,j) = zero
1407          veget_cov_max(:,j) = zero
1408       ENDWHERE
1409    ENDDO
1410
1411    IF ( do_now_stomate_lcchange ) THEN
1412       DO j=1,nvm
1413          WHERE ((1.-totfrac_nobio_new(:)) > min_sechiba)
1414             ! Pixels with vegetation
1415             veget_cov_max_new(:,j) = veget_max_new(:,j)/( 1.-totfrac_nobio_new(:) )
1416          ELSEWHERE
1417             ! Pixels without vegetation
1418             veget_cov_max_new(:,j) = zero
1419          ENDWHERE
1420       ENDDO
1421    ENDIF
1422
1423    !! 3.3 Adjust time step of GPP
1424    ! No GPP for bare soil
1425    gpp_d(:,1) = zero
1426    ! GPP per PFT
1427    DO j = 2,nvm   
1428       WHERE (veget_cov_max(:,j) > min_stomate)
1429          ! The PFT is available on the pixel
1430          gpp_d(:,j) =  gpp(:,j)/ veget_cov_max(:,j)* one_day/dt_sechiba 
1431       ELSEWHERE
1432          ! The PFT is absent on the pixel
1433          gpp_d(:,j) = zero
1434       ENDWHERE
1435    ENDDO
1436
1437    !! 3.4 The first time step of the first day of the month 
1438    ! implies that the month is over
1439    IF ( day == 1 .AND. sec .LT. dt_sechiba ) THEN
1440       EndOfMonth=.TRUE.
1441    ELSE
1442       EndOfMonth=.FALSE.
1443    ENDIF
1444   
1445
1446  !! 4. Calculate variables for dt_stomate (i.e. "daily")
1447
1448    ! Note: If dt_days /= 1, then variables 'xx_daily' (eg. half-daily or bi-daily) are by definition
1449    ! not expressed on a daily basis. This is not a problem but could be
1450    ! confusing
1451
1452    !! 4.1 Accumulate instantaneous variables (do_slow=.FALSE.)
1453    ! Accumulate instantaneous variables (do_slow=.FALSE.) and eventually
1454    ! calculate daily mean value (do_slow=.TRUE.)
1455    CALL stomate_accu (kjpindex, nvm,  do_slow, humrel,        humrel_daily)
1456    CALL stomate_accu (kjpindex,    1, do_slow, litterhumdiag, litterhum_daily)
1457    CALL stomate_accu (kjpindex,    1, do_slow, t2m,           t2m_daily)
1458    CALL stomate_accu (kjpindex,    1, do_slow, temp_sol,      tsurf_daily)
1459    CALL stomate_accu (kjpindex, nbdl, do_slow, stempdiag,     tsoil_daily)
1460    CALL stomate_accu (kjpindex, nbdl, do_slow, shumdiag,      soilhum_daily)
1461    CALL stomate_accu (kjpindex,    1, do_slow, precip,        precip_daily)
1462    CALL stomate_accu (kjpindex, nvm,  do_slow, gpp_d,         gpp_daily)
1463!gmjc # trunk do not introduce snow in stomate
1464    CALL stomate_accu (kjpindex,    1, do_slow, precip_snow,        snowfall_daily)
1465    CALL stomate_accu (kjpindex,    1, do_slow, snow,        snowmass_daily)
1466    CALL stomate_accu (kjpindex,    1, do_slow, tmc_topgrass, tmc_topgrass_daily)
1467!end gmjc
1468   
1469    !! 4.2 Daily minimum temperature
1470    t2m_min_daily(:) = MIN( t2m_min(:), t2m_min_daily(:) )
1471
1472    !! 4.3 Calculate maintenance respiration
1473    ! Note: lai is passed as output argument to overcome previous problems with
1474    ! natural and agricultural vegetation types.
1475    CALL maint_respiration &
1476         & (kjpindex,lai,t2m,t2m_longterm,stempdiag,height,veget_cov_max, &
1477         & rprof,biomass,resp_maint_part_radia, &
1478!gmjc
1479         & sla_calc)
1480!end gmjc   
1481    ! Aggregate maintenance respiration across the different plant parts
1482    resp_maint_radia(:,:) = zero
1483    DO j=2,nvm
1484       DO k= 1, nparts
1485          resp_maint_radia(:,j) = resp_maint_radia(:,j) &
1486               & + resp_maint_part_radia(:,j,k)
1487       ENDDO
1488    ENDDO
1489   
1490    ! Maintenance respiration separated by plant parts
1491    resp_maint_part(:,:,:) = resp_maint_part(:,:,:) &
1492         & + resp_maint_part_radia(:,:,:)
1493   
1494    !! 4.4 Litter dynamics and litter heterothropic respiration
1495    ! Including: litter update, lignin content, PFT parts, litter decay,
1496    ! litter heterotrophic respiration, dead leaf soil cover.
1497    ! Note: there is no vertical discretisation in the soil for litter decay.
1498    turnover_littercalc(:,:,:,:) = turnover_daily(:,:,:,:) * dt_sechiba/one_day
1499    bm_to_littercalc(:,:,:,:)    = bm_to_litter(:,:,:,:) * dt_sechiba/one_day       
1500    CALL littercalc (kjpindex, &
1501         turnover_littercalc, bm_to_littercalc, &
1502         veget_cov_max, temp_sol, stempdiag, shumdiag, litterhumdiag, &
1503         litterpart, litter, &
1504!gmjc
1505         litter_avail, litter_not_avail, litter_avail_frac, &
1506!end gmjc
1507         dead_leaves, lignin_struc, &
1508         deadleaf_cover, resp_hetero_litter, &
1509         soilcarbon_input_inst, control_temp_inst, control_moist_inst, &
1510         matrixA, vectorB, &!)
1511!gmjc
1512         sla_calc,do_slow)
1513!end gmjc
1514   
1515    ! Heterothropic litter respiration during time step ::dt_sechiba @tex $(gC m^{-2})$ @endtex
1516    resp_hetero_litter(:,:) = resp_hetero_litter(:,:) * dt_sechiba/one_day
1517   
1518    !! 4.5 Soil carbon dynamics and soil heterotrophic respiration
1519    ! Note: there is no vertical discretisation in the soil for litter decay.
1520    CALL soilcarbon (kjpindex, clay, &
1521         soilcarbon_input_inst, control_temp_inst, control_moist_inst, &
1522         carbon, resp_hetero_soil, matrixA)
1523   
1524    ! Heterothropic soil respiration during time step ::dt_sechiba @tex $(gC m^{-2})$ @endtex
1525    resp_hetero_soil(:,:) = resp_hetero_soil(:,:) * dt_sechiba/one_day
1526
1527    ! Total heterothrophic respiration during time step ::dt_sechiba @tex $(gC m^{-2})$ @endtex
1528    resp_hetero_radia(:,:) = resp_hetero_litter(:,:) + resp_hetero_soil(:,:)
1529    resp_hetero_d(:,:) = resp_hetero_d(:,:) + resp_hetero_radia(:,:)
1530   
1531    !! 4.6 Accumulate instantaneous variables (do_slow=.FALSE.)
1532    ! Accumulate instantaneous variables (do_slow=.FALSE.) and eventually
1533    ! calculate daily mean value (do_slow=.TRUE.)
1534    CALL stomate_accu (kjpindex, nlevs, do_slow, control_moist_inst, control_moist_daily)
1535    CALL stomate_accu (kjpindex, nlevs, do_slow, control_temp_inst,  control_temp_daily)
1536
1537    DO i = 1,ncarb
1538       CALL stomate_accu (kjpindex, nvm, do_slow, soilcarbon_input_inst(:,i,:), soilcarbon_input_daily(:,i,:))
1539    ENDDO
1540   
1541!! 5. Daily processes - performed at the end of the day
1542   
1543    IF (do_slow) THEN
1544
1545       !! 5.1 Update lai
1546       ! Use lai from stomate
1547       ! ?? check if this is the only time ok_pheno is used??
1548       ! ?? Looks like it is the only time. But this variables probably is defined
1549       ! in stomate_constants or something, in which case, it is difficult to track.
1550       IF (ok_pheno) THEN
1551          !! 5.1.1 Update LAI
1552          ! Set lai of bare soil to zero
1553          lai(:,ibare_sechiba) = zero
1554          ! lai for all PFTs
1555          DO j = 2, nvm
1556             lai(:,j) = biomass(:,j,ileaf,icarbon)*sla_calc(:,j)
1557          ENDDO
1558          frac_age(:,:,:) = leaf_frac(:,:,:)
1559       ELSE 
1560          ! 5.1.2 Use a prescribed lai
1561          ! WARNING: code in setlai is identical to the lines above
1562          ! Update subroutine if LAI has to be forced
1563          CALL  setlai(kjpindex,lai) 
1564          frac_age(:,:,:) = zero
1565       ENDIF
1566
1567       !! 5.2 Calculate long-term "meteorological" and biological parameters
1568       ! mainly in support of calculating phenology. If ::EndOfYear=.TRUE.
1569       ! annual values are update (i.e. xx_lastyear).
1570       CALL season &
1571            &          (kjpindex, dt_days, EndOfYear, &
1572            &           veget_cov, veget_cov_max, &
1573            &           humrel_daily, t2m_daily, tsoil_daily, soilhum_daily, lalo, &
1574            &           precip_daily, npp_daily, biomass, &
1575            &           turnover_daily, gpp_daily, when_growthinit, &
1576            &           maxhumrel_lastyear, maxhumrel_thisyear, &
1577            &           minhumrel_lastyear, minhumrel_thisyear, &
1578            &           maxgppweek_lastyear, maxgppweek_thisyear, &
1579            &           gdd0_lastyear, gdd0_thisyear, &
1580            &           precip_lastyear, precip_thisyear, &
1581            &           lm_lastyearmax, lm_thisyearmax, &
1582            &           maxfpc_lastyear, maxfpc_thisyear, &
1583            &           humrel_month, humrel_week, t2m_longterm, tau_longterm, &
1584            &           t2m_month, t2m_week, tsoil_month, soilhum_month, &
1585            &           npp_longterm, turnover_longterm, gpp_week, &
1586            &           gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
1587            &           time_hum_min, hum_min_dormance, gdd_init_date, &
1588            &           gdd_from_growthinit, herbivores, &
1589            &           Tseason, Tseason_length, Tseason_tmp, &
1590            &           Tmin_spring_time, t2m_min_daily, begin_leaves, onset_date, &!)
1591!gmjc
1592            &           t2m_14, sla_calc)
1593!end gmjc
1594       
1595       !! 5.3 Use all processes included in stomate
1596
1597       !! 5.3.1  Activate stomate processes
1598       ! Activate stomate processes (the complete list of processes depends
1599       ! on whether the DGVM is used or not). Processes include: climate constraints
1600       ! for PFTs, PFT dynamics, Phenology, Allocation, NPP (based on GPP and
1601       ! authothropic respiration), fire, mortality, vmax, assimilation temperatures,
1602       ! all turnover processes, light competition, sapling establishment, lai and
1603       ! land cover change.
1604       CALL StomateLpj &
1605            &            (kjpindex, dt_days, &
1606            &             neighbours, resolution, &
1607            &             clay, herbivores, &
1608            &             tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
1609            &             litterhum_daily, soilhum_daily, &
1610            &             maxhumrel_lastyear, minhumrel_lastyear, &
1611            &             gdd0_lastyear, precip_lastyear, &
1612            &             humrel_month, humrel_week, t2m_longterm, t2m_month, t2m_week, &
1613            &             tsoil_month, soilhum_month, &
1614            &             gdd_m5_dormance,  gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
1615            &             turnover_longterm, gpp_daily, &
1616            &             time_hum_min, maxfpc_lastyear, resp_maint_part,&
1617            &             PFTpresent, age, fireindex, firelitter, &
1618            &             leaf_age, leaf_frac, biomass, ind, adapted, regenerate, &
1619            &             senescence, when_growthinit, litterpart, litter, &
1620!gmjc
1621       litter_avail, litter_not_avail, litter_avail_frac, &
1622!end gmjc
1623            &             dead_leaves, carbon, lignin_struc, &
1624            &             veget_cov_max, veget_cov_max_new, npp_longterm, lm_lastyearmax, &
1625            &             veget_lastlight, everywhere, need_adjacent, RIP_time, &
1626            &             lai, rprof,npp_daily, turnover_daily, turnover_time,&
1627            &             control_moist_inst, control_temp_inst, soilcarbon_input_inst, &
1628            &             co2_to_bm_dgvm, co2_fire, &
1629            &             resp_hetero_d, resp_maint_d, resp_growth_d, &
1630            &             height, deadleaf_cover, vcmax, &
1631            &             bm_to_litter,&
1632            &             prod10, prod100, flux10, flux100, &
1633            &             convflux, cflux_prod10, cflux_prod100, harvest_above, carb_mass_total, &
1634            &             fpc_max, matrixA, &
1635            &             Tseason, Tmin_spring_time, begin_leaves, onset_date, &
1636!gmjc
1637            &             wshtotsum, sr_ugb, compt_ugb, nb_ani, grazed_frac,&
1638            &             import_yield, sla_age1, t2m_14, sla_calc,snowfall_daily,MODULO(date,365),&
1639            &             N_limfert, when_growthinit_cut, nb_grazingdays, &
1640            &             EndOfYear, &          !! Flag set to true for the first sechiba
1641!JCADD top 5 layer grassland soil moisture for grazing
1642            &             humrel_daily,tmc_topgrass_daily,fc_grazing,snowmass_daily, &
1643            &             after_snow, after_wet, wet1day, wet2day)
1644!end gmjc
1645       
1646       !! 5.3.2 Calculate the total CO2 flux from land use change
1647       fco2_lu(:) = convflux(:) &
1648            &             + cflux_prod10(:)  &
1649            &             + cflux_prod100(:) &
1650            &             + harvest_above(:)
1651       
1652       !! 5.4 Calculate veget and veget_max
1653       veget_max(:,:) = zero 
1654       DO j = 1, nvm
1655          veget_max(:,j) = veget_max(:,j) + &
1656               & veget_cov_max(:,j) * ( 1.-totfrac_nobio(:) )
1657       ENDDO
1658       
1659       !! 5.5 Photosynthesis parameters
1660       assim_param(:,:,ivcmax) = zero
1661       DO j = 2,nvm
1662          assim_param(:,j,ivcmax) = vcmax(:,j)
1663       ENDDO
1664       
1665       !! 5.6 Update forcing variables for soil carbon
1666       IF (TRIM(Cforcing_name) /= 'NONE') THEN
1667          npp_tot(:) = 0
1668          DO j=2,nvm
1669             npp_tot(:) = npp_tot(:) + npp_daily(:,j)
1670          ENDDO
1671          ! ::nbyear Number of years saved for carbon spinup
1672          sf_time = MODULO(REAL(date,r_std)-1,one_year*REAL(nbyear,r_std))
1673          iatt=FLOOR(sf_time/dt_forcesoil) + 1
1674          IF (iatt == 0) iatt = iatt_old + 1
1675          IF ((iatt<iatt_old) .and. (.not. cumul_Cforcing)) THEN
1676             nforce(:)=0
1677             soilcarbon_input(:,:,:,:) = zero
1678             control_moist(:,:,:) = zero
1679             control_temp(:,:,:) = zero
1680             npp_equil(:,:) = zero
1681          ENDIF
1682          iatt_old = iatt
1683          ! Update forcing
1684          nforce(iatt) = nforce(iatt)+1
1685          soilcarbon_input(:,:,:,iatt) = soilcarbon_input(:,:,:,iatt) + soilcarbon_input_daily(:,:,:)
1686          control_moist(:,:,iatt) = control_moist(:,:,iatt) + control_moist_daily(:,:)
1687          control_temp(:,:,iatt) = control_temp(:,:,iatt) + control_temp_daily(:,:)
1688          npp_equil(:,iatt) = npp_equil(:,iatt) + npp_tot(:)
1689       ENDIF
1690       
1691       !! 5.8 Write forcing file if ::ok_co2=.TRUE.
1692       ! Note: if STOMATE is run in coupled mode the forcing file is written
1693       ! If run in stand-alone mode, the forcing file is read!
1694       IF ( ok_co2 .AND. TRIM(forcing_name) /= 'NONE' ) THEN
1695         
1696          !! 5.8.1 Convert GPP to sechiba time steps
1697          ! GPP is multiplied by coverage to obtain forcing @tex $(gC m^{-2} dt_stomate^{-1})$\f \end@tex $(m^2 m^{-2})$ @endtexonly
1698          ! @tex$ m^{-2}$ @endtex remains in the units because ::veget_cov_max is a fraction, not a
1699          ! surface area. In sechiba values are ponderated by surface and frac_no_bio.
1700          ! At the beginning of stomate, the units are converted.
1701          ! When we use forcesoil we call sechiba_main and so we need the have the same units as in sechiba.
1702          gpp_daily_x(:,:) = zero
1703          DO j = 2, nvm             
1704             gpp_daily_x(:,j) = gpp_daily_x(:,j) + &
1705              & gpp_daily(:,j) * dt_stomate / one_day * veget_cov_max(:,j)
1706          ENDDO
1707         
1708          ! Bare soil moisture availability has not been treated
1709          ! in STOMATE, update it here
1710          humrel_daily(:,ibare_sechiba) = humrel(:,ibare_sechiba)   
1711
1712          ! Update index to store the next forcing step in memory
1713          iisf = iisf+1
1714
1715          ! How many times have we treated this forcing state
1716          xn = REAL(nf_cumul(isf(iisf)),r_std)
1717         
1718          !! 5.8.2 Cumulate forcing variables
1719          ! Cumulate forcing variables (calculate average)
1720          ! Note: precipitation is multiplied by dt_stomate/one_day to be consistent with
1721          ! the units in sechiba
1722          IF (cumul_forcing) THEN
1723             clay_fm(:,iisf) = (xn*clay_fm(:,iisf)+clay(:))/(xn+1.)
1724             humrel_daily_fm(:,:,iisf) = &
1725                  & (xn*humrel_daily_fm(:,:,iisf) + humrel_daily(:,:))/(xn+1.)
1726             litterhum_daily_fm(:,iisf) = &
1727                  & (xn*litterhum_daily_fm(:,iisf)+litterhum_daily(:))/(xn+1.)
1728             t2m_daily_fm(:,iisf) = &
1729                  & (xn*t2m_daily_fm(:,iisf)+t2m_daily(:))/(xn+1.)
1730             t2m_min_daily_fm(:,iisf) = &
1731                  & (xn*t2m_min_daily_fm(:,iisf)+t2m_min_daily(:))/(xn+1.)
1732             tsurf_daily_fm(:,iisf) = &
1733                  & (xn*tsurf_daily_fm(:,iisf)+tsurf_daily(:))/(xn+1.)
1734             tsoil_daily_fm(:,:,iisf) = &
1735                  & (xn*tsoil_daily_fm(:,:,iisf)+tsoil_daily(:,:))/(xn+1.)
1736             soilhum_daily_fm(:,:,iisf) = &
1737                  & (xn*soilhum_daily_fm(:,:,iisf)+soilhum_daily(:,:))/(xn+1.)
1738             precip_fm(:,iisf) = &
1739                  & (xn*precip_fm(:,iisf)+precip_daily(:)*dt_stomate/one_day)/(xn+1.)
1740             gpp_daily_fm(:,:,iisf) = &
1741                  & (xn*gpp_daily_fm(:,:,iisf) + gpp_daily_x(:,:))/(xn+1.)
1742             veget_fm(:,:,iisf) = &
1743                  & (xn*veget_fm(:,:,iisf) + veget(:,:) )/(xn+1.)
1744             veget_max_fm(:,:,iisf) = &
1745                  & (xn*veget_max_fm(:,:,iisf) + veget_max(:,:) )/(xn+1.)
1746             lai_fm(:,:,iisf) = &
1747                  & (xn*lai_fm(:,:,iisf) + lai(:,:) )/(xn+1.)
1748          ELSE
1749             ! Here we just calculate the values
1750             clay_fm(:,iisf) = clay(:)
1751             humrel_daily_fm(:,:,iisf) = humrel_daily(:,:)
1752             litterhum_daily_fm(:,iisf) = litterhum_daily(:)
1753             t2m_daily_fm(:,iisf) = t2m_daily(:)
1754             t2m_min_daily_fm(:,iisf) =t2m_min_daily(:)
1755             tsurf_daily_fm(:,iisf) = tsurf_daily(:)
1756             tsoil_daily_fm(:,:,iisf) =tsoil_daily(:,:)
1757             soilhum_daily_fm(:,:,iisf) =soilhum_daily(:,:)
1758             precip_fm(:,iisf) = precip_daily(:)
1759             gpp_daily_fm(:,:,iisf) =gpp_daily_x(:,:)
1760             veget_fm(:,:,iisf) = veget(:,:)
1761             veget_max_fm(:,:,iisf) =veget_max(:,:)
1762             lai_fm(:,:,iisf) =lai(:,:)
1763          ENDIF
1764          nf_cumul(isf(iisf)) = nf_cumul(isf(iisf))+1
1765
1766          ! 5.8.3 Do we have to write the forcing states?
1767          IF (iisf == nsfm) THEN
1768
1769             !! 5.8.3.1 Write these forcing states
1770             CALL forcing_write(forcing_id,1,nsfm)
1771             ! determine which forcing states must be read
1772             isf(1) = isf(nsfm)+1
1773             IF ( isf(1) > nsft ) isf(1) = 1
1774             DO iisf = 2, nsfm
1775                isf(iisf) = isf(iisf-1)+1
1776                IF (isf(iisf) > nsft)  isf(iisf) = 1
1777             ENDDO
1778
1779             ! Read forcing variables - for debug use only
1780             ! CALL forcing_read(forcing_id,nsfm)
1781             iisf = 0
1782
1783          ENDIF
1784
1785       ENDIF
1786
1787
1788       !! 5.9 Compute daily CO2 flux (AR5 output - not essential)
1789       ! CO2 flux in @tex $gC m^{-2} s^{-1}$ @endtex (positive towards the atmosphere) is sum of:
1790       ! (1) heterotrophic respiration from ground + (2) maintenance respiration
1791       ! from the plants + (3) growth respiration from the plants + (4) co2
1792       ! emissions from fire - (5) co2 taken up in the DGVM to establish
1793       ! saplings - (6) co2 taken up by photosyntyhesis
1794       co2_flux_daily(:,:)=   &
1795            & resp_maint_d(:,:) + resp_growth_d(:,:) + resp_hetero_d(:,:) + &
1796            & co2_fire(:,:) - co2_to_bm_dgvm(:,:) - gpp_daily(:,:)
1797
1798     CALL xios_orchidee_send_field("nep",SUM(co2_flux_daily*veget_cov_max,dim=2)/1e3/one_day*contfrac)
1799
1800       IF ( hist_id_stom_IPCC > 0 ) THEN
1801          vartmp(:) = SUM(co2_flux_daily*veget_cov_max,dim=2)/1e3/one_day*contfrac
1802          CALL histwrite_p (hist_id_stom_IPCC, "nep", itime, &
1803               vartmp, kjpindex, hori_index)
1804       ENDIF
1805
1806       ! See 5.9 for details on NEP + fire. At the monthly time step also
1807       ! harvest and land use change are calculated
1808       co2_flux_monthly(:,:) = co2_flux_monthly(:,:) + co2_flux_daily(:,:)
1809       harvest_above_monthly(:) = harvest_above_monthly(:) + harvest_above(:)
1810       cflux_prod_monthly(:) = cflux_prod_monthly(:) + convflux(:) + & 
1811        & cflux_prod10(:) + cflux_prod100(:)
1812     
1813       !! 5.10 Compute monthly CO2 fluxes
1814       IF ( EndOfMonth ) THEN
1815          !! 5.10.1 Write history file for monthly fluxes
1816          CALL histwrite_p (hist_id_stomate, 'CO2FLUX', itime, &
1817               co2_flux_monthly, kjpindex*nvm, horipft_index)
1818         
1819          !?? I (=VB) translated the French, but the whole stuff does not make sense to me.
1820          ! If one deletes the montly cumulation,
1821          ! one should not forget this change in resolution(:,1)*resolution(:,2)*contfrac(:)
1822          ! Si on supprimer le cumul par mois,
1823          ! il ne faut pas oublier cette modif resolution(:,1)*resolution(:,2)*contfrac(:)
1824          ! Should be supressed, this is post-processing
1825          DO j=2, nvm
1826             co2_flux_monthly(:,j) = co2_flux_monthly(:,j)* &
1827                  resolution(:,1)*resolution(:,2)*contfrac(:)
1828          ENDDO
1829
1830          ! Should be supressed, this is post-processing
1831          ! ?? How does it differ from co2_flux_monthly??
1832          net_co2_flux_monthly = zero
1833          DO ji=1,kjpindex
1834             DO j=2,nvm
1835                net_co2_flux_monthly = net_co2_flux_monthly + &
1836                     &  co2_flux_monthly(ji,j)*veget_cov_max(ji,j)
1837             ENDDO
1838          ENDDO
1839
1840     
1841          !! 5.10.2 Cumulative fluxes of land use cover change, harvest and net biosphere production
1842          ! Parallel processing, gather the information from different processors. first argument is the lo
1843          ! local variable, the second argument is the global variable. bcast send it to all processors.
1844          net_cflux_prod_monthly_sum = &
1845              &  SUM(cflux_prod_monthly(:)*resolution(:,1)*resolution(:,2)*contfrac(:))*1e-15
1846          CALL reduce_sum(net_cflux_prod_monthly_sum,net_cflux_prod_monthly_tot)
1847          CALL bcast(net_cflux_prod_monthly_tot)
1848          net_harvest_above_monthly_sum = &
1849             &   SUM(harvest_above_monthly(:)*resolution(:,1)*resolution(:,2)*contfrac(:))*1e-15
1850          CALL reduce_sum(net_harvest_above_monthly_sum,net_harvest_above_monthly_tot)
1851          CALL bcast(net_harvest_above_monthly_tot)
1852          net_co2_flux_monthly = net_co2_flux_monthly*1e-15
1853          CALL reduce_sum(net_co2_flux_monthly,net_co2_flux_monthly_sum)
1854          CALL bcast(net_co2_flux_monthly_sum)
1855          net_biosp_prod_monthly_tot =  &
1856             & ( net_co2_flux_monthly_sum + net_cflux_prod_monthly_tot + &
1857             & net_harvest_above_monthly_tot )
1858         
1859          WRITE(numout,9010) 'GLOBAL net_cflux_prod_monthly    (Peta gC/month)  = ',net_cflux_prod_monthly_tot
1860          WRITE(numout,9010) 'GLOBAL net_harvest_above_monthly (Peta gC/month)  = ',net_harvest_above_monthly_tot
1861          WRITE(numout,9010) 'GLOBAL net_co2_flux_monthly      (Peta gC/month)  = ',net_co2_flux_monthly_sum
1862          WRITE(numout,9010) 'GLOBAL net_biosp_prod_monthly    (Peta gC/month)  = ',net_biosp_prod_monthly_tot
1863
18649010  FORMAT(A52,F17.14)
1865
1866          ! Reset Monthly values
1867          co2_flux_monthly(:,:) = zero
1868          harvest_above_monthly(:) = zero
1869          cflux_prod_monthly(:)    = zero
1870
1871       ENDIF ! Monthly processes - at the end of the month
1872       
1873       IF (spinup_analytic) THEN
1874          nbp_accu(:) = nbp_accu(:) + (-SUM(co2_flux_daily(:,:) * veget_max(:,:),dim=2) - (convflux(:) + cflux_prod10(:) + &
1875                    cflux_prod100(:))  - harvest_above(:))/1e3 
1876       ENDIF
1877
1878       !! 5.11 Reset daily variables
1879       humrel_daily(:,:) = zero
1880       litterhum_daily(:) = zero
1881       t2m_daily(:) = zero
1882       t2m_min_daily(:) = large_value
1883       tsurf_daily(:) = zero
1884       tsoil_daily(:,:) = zero
1885       soilhum_daily(:,:) = zero
1886       precip_daily(:) = zero
1887       gpp_daily(:,:) = zero
1888       resp_maint_part(:,:,:)=zero
1889       resp_hetero_d=zero
1890       IF (printlev >= 3) THEN
1891          WRITE(numout,*) 'stomate_main: daily processes done'
1892       ENDIF
1893
1894    ENDIF  ! Daily processes - at the end of the day
1895   
1896  !! 6. Outputs from Stomate
1897
1898    ! co2_flux receives a value from STOMATE only if STOMATE is activated.
1899    ! Otherwise, the calling hydrological module must do this itself.
1900
1901    !! 6.1 Respiration and fluxes
1902    resp_maint(:,:) = resp_maint_radia(:,:)*veget_cov_max(:,:)
1903    resp_maint(:,ibare_sechiba) = zero
1904    resp_growth(:,:)= resp_growth_d(:,:)*veget_cov_max(:,:)*dt_sechiba/one_day
1905    resp_hetero(:,:) = resp_hetero_radia(:,:)*veget_cov_max(:,:)
1906   
1907    !! 6.2 Derived CO2 fluxes
1908    ! CO2 flux in gC m^{-2} s^{-1} (positive towards the atmosphere) is sum of:
1909    ! (1) heterotrophic respiration from ground + (2) maintenance respiration
1910    ! from the plants + (3) growth respiration from the plants + (4) co2
1911    ! emissions from fire - (5) co2 taken up in the DGVM to establish
1912    ! saplings - (6) co2 taken up by photosyntyhesis
1913    co2_flux(:,:) = resp_hetero(:,:) + resp_maint(:,:) + resp_growth(:,:) &
1914         & + (co2_fire(:,:)-co2_to_bm_dgvm(:,:))*veget_cov_max(:,:)/one_day &
1915         & - gpp(:,:)
1916   
1917    temp_growth(:)=t2m_month(:)-tp_00 
1918   
1919  !! 7. Analytical spinup
1920
1921    IF (spinup_analytic) THEN
1922
1923       !! 7.1. Update V and U at sechiba time step
1924       DO m = 2,nvm
1925          DO j = 1,kjpindex 
1926             ! V <- A * V
1927             MatrixV(j,m,:,:) = MATMUL(matrixA(j,m,:,:),MatrixV(j,m,:,:))
1928             ! U <- A*U + B
1929             VectorU(j,m,:) = MATMUL(matrixA(j,m,:,:),VectorU(j,m,:)) + vectorB(j,m,:)
1930          ENDDO ! loop pixels
1931       ENDDO ! loop PFTS
1932
1933
1934       !! 7.2. What happened at the end of the year ?
1935       IF (EndOfYear) THEN
1936
1937          !
1938          ! 7.2.1 Increase the years counter every EndOfyear
1939          !
1940          global_years = global_years + 1 
1941
1942
1943          !
1944          ! 7.2.3 Is global_years is a multiple of the period time ?
1945          !
1946
1947          !
1948          ! 3.2.1 When global_years is a multiple of the spinup_period, we calculate :
1949          !       1) the mean nbp flux over the period. This value is restarted
1950          !       2) we solve the matrix system by Gauss Jordan method
1951          !       3) We test if a point is at equilibrium : if yes, we mark the point (ok_equilibrium array)
1952          !       4) Then we reset the matrix
1953          !       5) We erase the carbon_stock calculated by ORCHIDEE by the one found by the method
1954          IF( MOD(global_years, spinup_period) == 0 ) THEN
1955             WRITE(numout,*) 'Spinup analytic : Calculate if system is in equlibrium. global_years=',global_years
1956             ! The number total of days during the forcing period is given by :
1957             !    spinup_period*365 (we consider only the noleap calendar)
1958             nbp_flux(:) = nbp_accu(:) / ( spinup_period * 365.)
1959             ! Reset the values
1960             nbp_accu(:) = zero
1961
1962             carbon_stock(:,ibare_sechiba,:) = zero
1963             ! Prepare the matrix for the resolution
1964             ! Add a temporary matrix W which contains I-MatrixV
1965             ! we should take the opposite of matrixV and add the identitiy : we solve (I-MatrixV)*C = VectorU
1966             MatrixW(:,:,:,:) = moins_un * MatrixV(:,:,:,:)
1967             DO jv = 1,nbpools
1968                MatrixW(:,:,jv,jv) =  MatrixW(:,:,jv,jv) + un
1969             ENDDO
1970             carbon_stock(:,:,:) = VectorU(:,:,:)
1971
1972             !
1973             !  Solve the linear system
1974             !
1975             DO m = 2,nvm
1976                DO j = 1,kjpindex
1977                   ! the solution will be stored in VectorU : so it should be restarted before
1978                   ! loop over npts and nvm, so we solved npts*(nvm-1) (7,7) linear systems
1979                   CALL gauss_jordan_method(nbpools,MatrixW(j,m,:,:),carbon_stock(j,m,:))
1980                ENDDO ! loop pixels
1981             ENDDO ! loop PFTS
1982
1983             ! Reset temporary matrixW
1984             MatrixW(:,:,:,:) = zero 
1985
1986
1987             previous_stock(:,:,:) = current_stock(:,:,:)
1988             current_stock(:,:,:) = carbon_stock(:,:,:) 
1989             ! The relative error is calculated over the passive carbon pool (sum over the pfts) over the pixel.
1990             CALL error_L1_passive(kjpindex,nvm, nbpools, current_stock, previous_stock, veget_max, &
1991                  &                eps_carbon, carbon_eq)   
1992
1993             !! ok_equilibrium is saved,
1994             WHERE( carbon_eq(:) .AND. .NOT.(ok_equilibrium(:)) )
1995                ok_equilibrium(:) = .TRUE. 
1996             ENDWHERE
1997
1998             ! Reset matrixV for the pixel to the identity matrix and vectorU to zero
1999             MatrixV(:,:,:,:) = zero
2000             VectorU(:,:,:) = zero
2001             DO jv = 1,nbpools
2002                MatrixV(:,:,jv,jv) = un
2003             END DO
2004             WRITE(numout,*) 'Reset for matrixV and VectorU done'   
2005
2006             !! Write the values found in the standard outputs of ORCHIDEE
2007             litter(:,istructural,:,iabove,icarbon) = carbon_stock(:,:,istructural_above)
2008             litter(:,istructural,:,ibelow,icarbon) = carbon_stock(:,:,istructural_below)
2009             litter(:,imetabolic,:,iabove,icarbon) = carbon_stock(:,:,imetabolic_above)
2010             litter(:,imetabolic,:,ibelow,icarbon) = carbon_stock(:,:,imetabolic_below)
2011             carbon(:,iactive,:) = carbon_stock(:,:,iactive_pool)
2012             carbon(:,islow,:) = carbon_stock(:,:,islow_pool)
2013             carbon(:,ipassive,:) = carbon_stock(:,:,ipassive_pool) 
2014
2015             ! Final step, test if all points at the local domain are at equilibrium
2016             ! The simulation can be stopped when all local domains have reached the equilibrium
2017             IF(ALL(ok_equilibrium)) THEN
2018                WRITE(numout,*) 'Spinup analytic : Equilibrium for carbon pools is reached for current local domain'
2019             ELSE
2020                WRITE(numout,*) 'Spinup analytic : Equilibrium for carbon pools is not yet reached for current local domain'
2021             END IF
2022          ENDIF ! ( MOD(global_years,spinup_period) == 0)
2023       ENDIF ! (EndOfYear)
2024
2025    ENDIF !(spinup_analytic)
2026   
2027    IF (printlev >= 4) WRITE(numout,*) 'Leaving stomate_main'
2028
2029  END SUBROUTINE stomate_main
2030
2031!! ================================================================================================================================
2032!! SUBROUTINE   : stomate_finalize
2033!!
2034!>\BRIEF        Write variables to restart file
2035!!
2036!! DESCRIPTION  : Write variables to restart file
2037!! RECENT CHANGE(S) : None
2038!!
2039!! MAIN OUTPUT VARIABLE(S):
2040!!
2041!! REFERENCES   :
2042!!
2043!! \n
2044!_ ================================================================================================================================
2045
2046  SUBROUTINE stomate_finalize (kjit, kjpindex, index, clay, assim_param) 
2047   
2048    IMPLICIT NONE
2049   
2050    !! 0. Variable and parameter declaration
2051    !! 0.1 Input variables
2052    INTEGER(i_std),INTENT(in)                       :: kjit              !! Time step number (unitless)
2053    INTEGER(i_std),INTENT(in)                       :: kjpindex          !! Domain size - terrestrial pixels only (unitless)
2054    INTEGER(i_std),DIMENSION(kjpindex),INTENT(in)   :: index             !! Indices of the terrestrial pixels only (unitless)
2055    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
2056    REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(in) :: assim_param    !! min+max+opt temperatures (K) & vmax for
2057                                                                            !! photosynthesis 
2058    !! 0.4 Local variables
2059    REAL(r_std)                                   :: dt_days_read             !! STOMATE time step read in restart file (days)
2060    INTEGER(i_std)                                :: l,k,ji, jv, i, j, m      !! indices   
2061    REAL(r_std),PARAMETER                         :: max_dt_days = 5.         !! Maximum STOMATE time step (days)
2062    REAL(r_std)                                   :: hist_days                !! Writing frequency for history file (days)
2063    REAL(r_std),DIMENSION(0:nbdl)                 :: z_soil                   !! Variable to store depth of the different soil layers (m)
2064    REAL(r_std),DIMENSION(kjpindex)               :: cvegtot                  !! Total "vegetation" cover (unitless)
2065    REAL(r_std),DIMENSION(kjpindex)               :: precip                   !! Total liquid and solid precipitation 
2066                                                                              !! @tex $(??mm dt_stomate^{-1})$ @endtex
2067    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_d                    !! Gross primary productivity per ground area
2068                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex 
2069    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_daily_x              !! "Daily" gpp for teststomate 
2070                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex
2071    REAL(r_std),DIMENSION(kjpindex,nvm)           :: resp_hetero_litter       !! Litter heterotrophic respiration per ground area
2072                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex 
2073                                                                              !! ??Same variable is also used to
2074                                                                              !! store heterotrophic respiration per ground area
2075                                                                              !! over ::dt_sechiba??
2076    REAL(r_std),DIMENSION(kjpindex,nvm)           :: resp_hetero_soil         !! soil heterotrophic respiration 
2077                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex
2078    REAL(r_std),DIMENSION(kjpindex,nvm)           :: veget_cov                !! Fractional coverage: actually share of the pixel
2079                                                                              !! covered by a PFT (fraction of ground area),
2080                                                                              !! taking into account LAI ??(= grid scale fpc)??
2081    REAL(r_std),DIMENSION(kjpindex,nvm)           :: vcmax                    !! Maximum rate of carboxylation
2082                                                                              !! @tex $(\mumol m^{-2} s^{-1})$ @endtex
2083    REAL(r_std),DIMENSION(kjpindex,nlevs)         :: control_moist_inst       !! Moisture control of heterotrophic respiration
2084                                                                              !! (0-1, unitless)
2085    REAL(r_std),DIMENSION(kjpindex,nlevs)         :: control_temp_inst        !! Temperature control of heterotrophic
2086                                                                              !! respiration, above and below (0-1, unitless)
2087    REAL(r_std),DIMENSION(kjpindex,ncarb,nvm)     :: soilcarbon_input_inst    !! Quantity of carbon going into carbon pools from
2088                                                                              !! litter decomposition
2089                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex
2090   
2091    INTEGER(i_std)                                :: ier                      !! Check errors in netcdf call (unitless)
2092    REAL(r_std)                                   :: sf_time                  !! Intermediate variable to calculate current time
2093                                                                              !! step
2094    INTEGER(i_std)                                :: max_totsize              !! Memory management - maximum memory size (Mb)
2095    INTEGER(i_std)                                :: totsize_1step            !! Memory management - memory required to store one
2096                                                                              !! time step on one processor (Mb)
2097    INTEGER(i_std)                                :: totsize_tmp              !! Memory management - memory required to store one
2098                                                                              !! time step on all processors(Mb)
2099    REAL(r_std)                                   :: xn                       !! How many times have we treated in this forcing
2100    REAL(r_std), DIMENSION(kjpindex)              :: vartmp                   !! Temporary variable
2101    INTEGER(i_std)                                :: vid                      !! Variable identifer of netCDF (unitless)
2102    INTEGER(i_std)                                :: nneigh                   !! Number of neighbouring pixels
2103    INTEGER(i_std)                                :: direct                   !! ??
2104    INTEGER(i_std),DIMENSION(ndm)                 :: d_id                     !! ??
2105    REAL(r_std),DIMENSION(nbp_glo)                :: clay_g                   !! Clay fraction of soil (0-1, unitless), parallel
2106                                                                              !! computing
2107    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:,:)    :: soilcarbon_input_g       !! Quantity of carbon going into carbon pools from
2108                                                                              !! litter decomposition 
2109                                                                              !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex, parallel
2110                                                                              !! computing
2111    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)      :: control_moist_g          !! Moisture control of heterotrophic respiration
2112                                                                              !! (0-1, unitless), parallel computing
2113    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)      :: control_temp_g           !! Temperature control of heterotrophic respiration
2114                                                                              !! (0-1, unitless), parallel computing
2115    REAL(r_std),ALLOCATABLE,DIMENSION(:,:)        :: npp_equil_g              !! Equilibrium NPP written to forcesoil
2116                                                                              !! @tex $(gC m^{-2} year^{-1})$ @endtex, parallel
2117                                                                              !! computing
2118
2119    REAL(r_std)                                   :: net_cflux_prod_monthly_sum    !! AR5 output?? gC m2 month-1 (one variable for
2120                                                                                   !! reduce_sum and one for bcast??), parallel
2121                                                                                   !! computing
2122    REAL(r_std)                                   :: net_cflux_prod_monthly_tot    !! AR5 output?? gC m2 month-1 (one variable for
2123                                                                                   !! reduce_sum and one for bcast??), parallel
2124                                                                                   !! computing
2125    REAL(r_std)                                   :: net_harvest_above_monthly_sum !! AR5 output?? gC m2 month-1 (one variable for
2126                                                                                   !! reduce_sum and one for bcast??), parallel
2127                                                                                   !! computing
2128    REAL(r_std)                                   :: net_harvest_above_monthly_tot !! AR5 output?? gC m2 month-1 (one variable for
2129                                                                                   !! reduce_sum and one for bcast??), parallel
2130                                                                                   !! computing
2131    REAL(r_std)                                   :: net_biosp_prod_monthly_sum    !! AR5 output?? gC m2 month-1 (one variable for
2132                                                                                   !! reduce_sum and one for bcast??), parallel
2133                                                                                   !! computing
2134    REAL(r_std)                                   :: net_biosp_prod_monthly_tot    !! AR5 output?? gC m2 month-1 (one variable for
2135                                                                                   !! reduce_sum and one for bcast??), parallel
2136                                                                                   !! computing
2137    REAL(r_std), DIMENSION(kjpindex,nvm,nbpools)  :: carbon_stock                  !! Array containing the carbon stock for each pool
2138                                                                                   !! used by ORCHIDEE
2139
2140!_ ================================================================================================================================
2141   
2142    !! 1. Write restart file for stomate
2143    IF (printlev>=3) WRITE (numout,*) 'Write restart file for STOMATE'
2144       
2145    CALL writerestart &
2146         (kjpindex, index, &
2147         dt_days, date, &
2148         ind, adapted, regenerate, &
2149         humrel_daily, gdd_init_date, litterhum_daily, &
2150         t2m_daily, t2m_min_daily, tsurf_daily, tsoil_daily, &
2151         soilhum_daily, precip_daily, &
2152         gpp_daily, npp_daily, turnover_daily, &
2153         humrel_month, humrel_week, &
2154         t2m_longterm, tau_longterm, t2m_month, t2m_week, &
2155         tsoil_month, soilhum_month, fireindex, firelitter, &
2156         maxhumrel_lastyear, maxhumrel_thisyear, &
2157         minhumrel_lastyear, minhumrel_thisyear, &
2158         maxgppweek_lastyear, maxgppweek_thisyear, &
2159         gdd0_lastyear, gdd0_thisyear, &
2160         precip_lastyear, precip_thisyear, &
2161         gdd_m5_dormance, gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
2162         PFTpresent, npp_longterm, lm_lastyearmax, lm_thisyearmax, &
2163         maxfpc_lastyear, maxfpc_thisyear, &
2164         turnover_longterm, gpp_week, biomass, resp_maint_part, &
2165         leaf_age, leaf_frac, &
2166         senescence, when_growthinit, age, &
2167         resp_hetero_d, resp_maint_d, resp_growth_d, co2_fire, co2_to_bm_dgvm, &
2168         veget_lastlight, everywhere, need_adjacent, &
2169         RIP_time, &
2170         time_hum_min, hum_min_dormance, &
2171         litterpart, litter, dead_leaves, &
2172         carbon, lignin_struc,turnover_time,&
2173         prod10,prod100,flux10, flux100, &
2174         convflux, cflux_prod10, cflux_prod100, bm_to_litter, carb_mass_total, &
2175         Tseason, Tseason_length, Tseason_tmp, &
2176         Tmin_spring_time, begin_leaves, onset_date, &
2177         global_years, ok_equilibrium, nbp_accu, nbp_flux, &
2178         MatrixV, VectorU, previous_stock, current_stock, assim_param, &
2179    !gmjc
2180        & wshtotsum, sr_ugb, sla_calc, nb_ani, grazed_frac, &
2181        & import_yield, t2m_14, litter_not_avail, nb_grazingdays, &
2182        & after_snow, after_wet, wet1day, wet2day)
2183    !end gmjc)
2184   
2185    !! 2. Write file with variables that force general processes in stomate
2186    IF (ok_co2 .AND. allow_forcing_write ) THEN
2187       IF ( TRIM(forcing_name) /= 'NONE' ) THEN 
2188          CALL forcing_write(forcing_id,1,iisf)
2189          ! Close forcing file
2190          IF (is_root_prc) ier = NF90_CLOSE (forcing_id)
2191          forcing_id=-1
2192       END IF
2193    END IF
2194   
2195    !! 3. Collect variables that force the soil processes in stomate
2196    IF (TRIM(Cforcing_name) /= 'NONE' ) THEN 
2197       
2198       !! Collet variables
2199       WRITE(numout,*) &
2200            &      'stomate: writing the forcing file for carbon spinup'
2201       DO iatt = 1, nparan*nbyear
2202          IF ( nforce(iatt) > 0 ) THEN
2203             soilcarbon_input(:,:,:,iatt) = &
2204                  & soilcarbon_input(:,:,:,iatt)/REAL(nforce(iatt),r_std)
2205             control_moist(:,:,iatt) = &
2206                  & control_moist(:,:,iatt)/REAL(nforce(iatt),r_std)
2207             control_temp(:,:,iatt) = &
2208                  & control_temp(:,:,iatt)/REAL(nforce(iatt),r_std)
2209             npp_equil(:,iatt) = &
2210                  & npp_equil(:,iatt)/REAL(nforce(iatt),r_std)
2211          ELSE
2212             WRITE(numout,*) &
2213                  &         'We have no soil carbon forcing data for this time step:', &
2214                  &         iatt
2215             WRITE(numout,*) ' -> we set them to zero'
2216             soilcarbon_input(:,:,:,iatt) = zero
2217             control_moist(:,:,iatt) = zero
2218             control_temp(:,:,iatt) = zero
2219             npp_equil(:,iatt) = zero
2220          ENDIF
2221       ENDDO
2222       
2223       ! Allocate memory for parallel computing
2224       IF (is_root_prc) THEN
2225          ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan*nbyear))
2226          ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan*nbyear))
2227          ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan*nbyear))
2228          ALLOCATE(npp_equil_g(nbp_glo,nparan*nbyear))
2229       ENDIF
2230       
2231       ! Gather distributed variables
2232       CALL gather(clay,clay_g)
2233       CALL gather(soilcarbon_input,soilcarbon_input_g)
2234       CALL gather(control_moist,control_moist_g)
2235       CALL gather(control_temp,control_temp_g)
2236       CALL gather(npp_equil,npp_equil_g)
2237       
2238       !! Create netcdf
2239       ! Create, define and populate a netcdf file containing the forcing data.
2240       ! For the root processor only (parallel computing). NF90_ are functions
2241       ! from and external library. 
2242       IF (is_root_prc) THEN
2243          WRITE (numout,*) 'Create Cforcing file : ',TRIM(Cforcing_name)
2244          ! Create new netCDF dataset
2245          ier = NF90_CREATE (TRIM(Cforcing_name),NF90_64BIT_OFFSET ,Cforcing_id)
2246          IF (ier /= NF90_NOERR) THEN
2247             WRITE (numout,*) 'Error in creating Cforcing file : ',TRIM(Cforcing_name)
2248             CALL ipslerr_p (3,'stomate_finalize', &
2249                  &        'PROBLEM creating Cforcing file', &
2250                  &        NF90_STRERROR(ier),'')
2251          END IF
2252         
2253          ! Add variable attribute
2254          ! Note ::nbp_glo is the number of global continental points
2255          ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
2256               &                        'kjpindex',REAL(nbp_glo,r_std))
2257          ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
2258               &                        'nparan',REAL(nparan,r_std))
2259          ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
2260               &                        'nbyear',REAL(nbyear,r_std))
2261         
2262          ! Add new dimension
2263          ier = NF90_DEF_DIM (Cforcing_id,'points',nbp_glo,d_id(1))
2264          ier = NF90_DEF_DIM (Cforcing_id,'carbtype',ncarb,d_id(2))
2265          ier = NF90_DEF_DIM (Cforcing_id,'vegtype',nvm,d_id(3))
2266          ier = NF90_DEF_DIM (Cforcing_id,'level',nlevs,d_id(4))
2267          ier = NF90_DEF_DIM (Cforcing_id,'time_step',NF90_UNLIMITED,d_id(5))
2268         
2269          ! Add new variable
2270          ier = NF90_DEF_VAR (Cforcing_id,'points',    r_typ,d_id(1),vid)
2271          ier = NF90_DEF_VAR (Cforcing_id,'carbtype',  r_typ,d_id(2),vid)
2272          ier = NF90_DEF_VAR (Cforcing_id,'vegtype',   r_typ,d_id(3),vid)
2273          ier = NF90_DEF_VAR (Cforcing_id,'level',     r_typ,d_id(4),vid)
2274          ier = NF90_DEF_VAR (Cforcing_id,'time_step', r_typ,d_id(5),vid)
2275          ier = NF90_DEF_VAR (Cforcing_id,'index',     r_typ,d_id(1),vid)
2276          ier = NF90_DEF_VAR (Cforcing_id,'clay',      r_typ,d_id(1),vid)
2277          ier = NF90_DEF_VAR (Cforcing_id,'soilcarbon_input',r_typ, &
2278               &                        (/ d_id(1),d_id(2),d_id(3),d_id(5) /),vid)
2279          ier = NF90_DEF_VAR (Cforcing_id,'control_moist',r_typ, &
2280               &                        (/ d_id(1),d_id(4),d_id(5) /),vid)
2281          ier = NF90_DEF_VAR (Cforcing_id,'control_temp',r_typ, &
2282               &                        (/ d_id(1),d_id(4),d_id(5) /),vid)
2283          ier = NF90_DEF_VAR (Cforcing_id,'npp_equil',r_typ, &
2284               &                        (/ d_id(1),d_id(5) /),vid)
2285          ier = NF90_ENDDEF (Cforcing_id)
2286         
2287          ! Given the name of a varaible, nf90_inq_varid finds the variable
2288          ! ID (::vid). Put data value(s) into variable ::vid
2289          ier = NF90_INQ_VARID (Cforcing_id,'points',vid)
2290          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2291               &                          (/(REAL(i,r_std),i=1,nbp_glo)/))
2292          ier = NF90_INQ_VARID (Cforcing_id,'carbtype',vid)
2293          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2294               &                        (/(REAL(i,r_std),i=1,ncarb)/))
2295          ier = NF90_INQ_VARID (Cforcing_id,'vegtype',vid)
2296          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2297               &                            (/(REAL(i,r_std),i=1,nvm)/))
2298          ier = NF90_INQ_VARID (Cforcing_id,'level',vid)
2299          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2300               &                          (/(REAL(i,r_std),i=1,nlevs)/))
2301          ier = NF90_INQ_VARID (Cforcing_id,'time_step',vid)
2302          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2303               &                          (/(REAL(i,r_std),i=1,nparan*nbyear)/))
2304          ier = NF90_INQ_VARID (Cforcing_id,'index',vid)
2305          ier = NF90_PUT_VAR (Cforcing_id,vid, REAL(index_g,r_std) )
2306          ier = NF90_INQ_VARID (Cforcing_id,'clay',vid)
2307          ier = NF90_PUT_VAR (Cforcing_id,vid, clay_g )
2308          ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',vid)
2309          ier = NF90_PUT_VAR (Cforcing_id,vid, soilcarbon_input_g )
2310          ier = NF90_INQ_VARID (Cforcing_id,'control_moist',vid)
2311          ier = NF90_PUT_VAR (Cforcing_id,vid, control_moist_g )
2312          ier = NF90_INQ_VARID (Cforcing_id,'control_temp',vid)
2313          ier = NF90_PUT_VAR (Cforcing_id,vid, control_temp_g )
2314          ier = NF90_INQ_VARID (Cforcing_id,'npp_equil',vid)
2315          ier = NF90_PUT_VAR (Cforcing_id,vid, npp_equil_g )
2316         
2317          ! Close netCDF
2318          ier = NF90_CLOSE (Cforcing_id)
2319          IF (ier /= NF90_NOERR) THEN
2320             CALL ipslerr_p (3,'stomate_finalize', &
2321                  &        'PROBLEM in closing Cforcing file', &
2322                  &        NF90_STRERROR(ier),'')
2323          END IF
2324         
2325          Cforcing_id = -1
2326       ENDIF
2327
2328       ! Clear memory
2329       IF (is_root_prc) THEN
2330          DEALLOCATE(soilcarbon_input_g)
2331          DEALLOCATE(control_moist_g)
2332          DEALLOCATE(control_temp_g)
2333          DEALLOCATE(npp_equil_g)
2334       ENDIF
2335       
2336    ENDIF
2337 
2338  END SUBROUTINE stomate_finalize
2339
2340
2341!! ================================================================================================================================
2342!! SUBROUTINE   : stomate_init
2343!!
2344!>\BRIEF        The routine is called only at the first simulation. At that
2345!! time settings and flags are read and checked for internal consistency and
2346!! memory is allocated for the variables in stomate.
2347!!
2348!! DESCRIPTION  : The routine reads the
2349!! following flags from the run definition file:
2350!! -ipd (index of grid point for online diagnostics)\n
2351!! -ok_herbivores (flag to activate herbivores)\n
2352!! -treat_expansion (flag to activate PFT expansion across a pixel\n
2353!! -harvest_agri (flag to harvest aboveground biomass from agricultural PFTs)\n
2354!! \n
2355!! Check for inconsistent setting between the following flags:
2356!! -ok_stomate\n
2357!! -ok_dgvm\n
2358!! -ok_co2\n
2359!! \n
2360!! Memory is allocated for all the variables of stomate and new indexing tables
2361!! are build. New indexing tables are needed because a single pixel can conatin
2362!! several PFTs. The new indexing tables have separate indices for the different
2363!! PFTs. Similar index tables are build for land use cover change.\n
2364!! \n
2365!! Several global variables and land cover change variables are initialized to
2366!! zero.\n
2367!!
2368!! RECENT CHANGE(S) : None
2369!!
2370!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
2371!! variables. However, the routine allocates memory and builds new indexing
2372!! variables for later use.\n
2373!!
2374!! REFERENCE(S) : None
2375!!
2376!! FLOWCHART    : None
2377!! \n
2378!_ ================================================================================================================================
2379
2380  SUBROUTINE stomate_init &
2381       &  (kjpij, kjpindex, index, lalo, &
2382       &   rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
2383
2384  !! 0. Variable and parameter declaration
2385
2386    !! 0.1 Input variables
2387
2388    INTEGER(i_std),INTENT(in)                    :: kjpij             !! Total size of the un-compressed grid, including
2389                                                                      !! oceans (unitless)
2390    INTEGER(i_std),INTENT(in)                    :: kjpindex          !! Domain size - number of terrestrial pixels
2391                                                                      !! (unitless)
2392    INTEGER(i_std),INTENT(in)                    :: rest_id_stom      !! STOMATE's _Restart_ file identifier
2393    INTEGER(i_std),INTENT(in)                    :: hist_id_stom      !! STOMATE's _history_ file identifier
2394    INTEGER(i_std),INTENT(in)                    :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier
2395    INTEGER(i_std),DIMENSION(kjpindex),INTENT(in):: index             !! Indices of the terrestrial pixels on the global
2396                                                                      !! map
2397    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in) :: lalo              !! Geogr. coordinates (latitude,longitude) (degrees)
2398   
2399    !! 0.2 Output variables
2400
2401    !! 0.3 Modified variables
2402
2403    !! 0.4 Local variables
2404
2405    LOGICAL                                      :: l_error           !! Check errors in netcdf call
2406    INTEGER(i_std)                               :: ier               !! Check errors in netcdf call
2407    INTEGER(i_std)                               :: ji,j,ipd,l        !! Indices
2408!_ ================================================================================================================================
2409   
2410  !! 1. Online diagnostics
2411
2412    IF ( kjpindex > 0 ) THEN
2413       !Config  Key  = STOMATE_DIAGPT
2414       !Config  Desc = Index of grid point for online diagnostics
2415       !Config If    = OK_STOMATE
2416       !Config  Def  = 1
2417       !Config  Help = This is the index of the grid point which
2418       !               will be used for online diagnostics.
2419       !Config Units = [-]
2420       ! By default ::ipd is set to 1
2421       ipd = 1
2422       ! Get ::ipd from run definition file
2423       CALL getin_p('STOMATE_DIAGPT',ipd)
2424       ipd = MIN( ipd, kjpindex )
2425       WRITE(numout,*) 'Stomate: '
2426       WRITE(numout,*) '  Index of grid point for online diagnostics: ',ipd
2427       WRITE(numout,*) '  Lon, lat:',lalo(ipd,2),lalo(ipd,1)
2428       WRITE(numout,*) '  Index of this point on GCM grid: ',index(ipd)
2429       !
2430    ENDIF
2431   
2432  !! 2. Check consistency of flags
2433
2434    IF ( ( .NOT. ok_stomate ) .AND. ok_dgvm ) THEN
2435       WRITE(numout,*) 'Cannot do dynamical vegetation without STOMATE.'
2436       WRITE(numout,*) 'Inconsistency between ::ok_stomate and ::ok_dgvm'
2437       WRITE(numout,*) 'Stop: fatal error'
2438       STOP
2439    ENDIF
2440
2441    IF ((.NOT.ok_co2).AND.ok_stomate) THEN
2442       WRITE(numout,*) 'Cannot call STOMATE without GPP.'
2443       WRITE(numout,*) 'Inconsistency between ::ok_stomate and ::ok_co2'
2444       WRITE(numout,*) 'Stop: fatal error'
2445       STOP
2446    ENDIF
2447
2448  !! 3. Communicate settings
2449   
2450    WRITE(numout,*) 'stomate first call - overview of the activated flags:'
2451    WRITE(numout,*) '  Photosynthesis: ', ok_co2
2452    WRITE(numout,*) '  STOMATE: ', ok_stomate
2453    WRITE(numout,*) '  LPJ: ', ok_dgvm
2454   
2455  !! 4. Allocate memory for STOMATE's variables
2456
2457    l_error = .FALSE.
2458
2459    ALLOCATE(veget_cov_max(kjpindex,nvm),stat=ier)
2460    l_error = l_error .OR. (ier /= 0)
2461    IF (l_error) THEN
2462       WRITE(numout,*) 'Memory allocation error for veget_cov_max. We stop. We need kjpindex*nvm words',kjpindex,nvm
2463       STOP 'stomate_init'
2464    ENDIF
2465
2466    ALLOCATE(ind(kjpindex,nvm),stat=ier)
2467    l_error = l_error .OR. (ier /= 0)
2468    IF (l_error) THEN
2469       WRITE(numout,*) 'Memory allocation error for ind. We stop. We need kjpindex*nvm words',kjpindex,nvm
2470       STOP 'stomate_init'
2471    ENDIF
2472
2473    ALLOCATE(adapted(kjpindex,nvm),stat=ier)
2474    l_error = l_error .OR. (ier /= 0)
2475    IF (l_error) THEN
2476       WRITE(numout,*) 'Memory allocation error for adapted. We stop. We need kjpindex*nvm words',kjpindex,nvm
2477       STOP 'stomate_init'
2478    ENDIF
2479
2480    ALLOCATE(regenerate(kjpindex,nvm),stat=ier)
2481    l_error = l_error .OR. (ier /= 0)
2482    IF (l_error) THEN
2483       WRITE(numout,*) 'Memory allocation error for regenerate. We stop. We need kjpindex*nvm words',kjpindex,nvm
2484       STOP 'stomate_init'
2485    ENDIF
2486
2487    ALLOCATE(humrel_daily(kjpindex,nvm),stat=ier)
2488    l_error = l_error .OR. (ier /= 0)
2489    IF (l_error) THEN
2490       WRITE(numout,*) 'Memory allocation error for humrel_daily. We stop. We need kjpindex*nvm words',kjpindex,nvm
2491       STOP 'stomate_init'
2492    ENDIF
2493
2494    ALLOCATE(litterhum_daily(kjpindex),stat=ier)
2495    l_error = l_error .OR. (ier /= 0)
2496    IF (l_error) THEN
2497       WRITE(numout,*) 'Memory allocation error for litterhum_daily. We stop. We need kjpindex words',kjpindex
2498       STOP 'stomate_init'
2499    ENDIF
2500
2501    ALLOCATE(t2m_daily(kjpindex),stat=ier)
2502    l_error = l_error .OR. (ier /= 0)
2503    IF (l_error) THEN
2504       WRITE(numout,*) 'Memory allocation error for t2m_daily. We stop. We need kjpindex words',kjpindex
2505       STOP 'stomate_init'
2506    ENDIF
2507
2508    ALLOCATE(t2m_min_daily(kjpindex),stat=ier)
2509    l_error = l_error .OR. (ier /= 0)
2510    IF (l_error) THEN
2511       WRITE(numout,*) 'Memory allocation error for t2m_min_daily. We stop. We need kjpindex words',kjpindex
2512       STOP 'stomate_init'
2513    ENDIF
2514
2515    ALLOCATE(tsurf_daily(kjpindex),stat=ier)
2516    l_error = l_error .OR. (ier /= 0)
2517    IF (l_error) THEN
2518       WRITE(numout,*) 'Memory allocation error for tsurf_daily. We stop. We need kjpindex words',kjpindex
2519       STOP 'stomate_init'
2520    ENDIF
2521
2522    ALLOCATE(tsoil_daily(kjpindex,nbdl),stat=ier)
2523    l_error = l_error .OR. (ier /= 0)
2524    IF (l_error) THEN
2525       WRITE(numout,*) 'Memory allocation error for tsoil_daily. We stop. We need kjpindex*nbdl words',kjpindex,nbdl
2526       STOP 'stomate_init'
2527    ENDIF
2528
2529    ALLOCATE(soilhum_daily(kjpindex,nbdl),stat=ier)
2530    l_error = l_error .OR. (ier /= 0)
2531    IF (l_error) THEN
2532       WRITE(numout,*) 'Memory allocation error for soilhum_daily. We stop. We need kjpindex*nbdl words',kjpindex,nbdl
2533       STOP 'stomate_init'
2534    ENDIF
2535
2536    ALLOCATE(precip_daily(kjpindex),stat=ier)
2537    l_error = l_error .OR. (ier /= 0)
2538    IF (l_error) THEN
2539       WRITE(numout,*) 'Memory allocation error for precip_daily. We stop. We need kjpindex words',kjpindex,nvm
2540       STOP 'stomate_init'
2541    ENDIF
2542
2543    ALLOCATE(gpp_daily(kjpindex,nvm),stat=ier)
2544    l_error = l_error .OR. (ier /= 0)
2545    IF (l_error) THEN
2546       WRITE(numout,*) 'Memory allocation error for gpp_daily. We stop. We need kjpindex*nvm words',kjpindex,nvm
2547       STOP 'stomate_init'
2548    ENDIF
2549
2550    ALLOCATE(npp_daily(kjpindex,nvm),stat=ier)
2551    l_error = l_error .OR. (ier /= 0)
2552    IF (l_error) THEN
2553       WRITE(numout,*) 'Memory allocation error for npp_daily. We stop. We need kjpindex*nvm words',kjpindex,nvm
2554       STOP 'stomate_init'
2555    ENDIF
2556
2557    ALLOCATE(turnover_daily(kjpindex,nvm,nparts,nelements),stat=ier)
2558    l_error = l_error .OR. (ier /= 0)
2559    IF (l_error) THEN
2560       WRITE(numout,*) 'Memory allocation error for turnover_daily. We stop. We need kjpindex*nvm*nparts*nelements words', &
2561       &   kjpindex,nvm,nparts,nelements
2562       STOP 'stomate_init'
2563    ENDIF
2564
2565    ALLOCATE(turnover_littercalc(kjpindex,nvm,nparts,nelements),stat=ier)
2566    l_error = l_error .OR. (ier /= 0)
2567    IF (l_error) THEN
2568       WRITE(numout,*) 'Memory allocation error for turnover_littercalc. We stop. We need kjpindex*nvm*nparts*nelements words', & 
2569        &  kjpindex,nvm,nparts,nelements
2570       STOP 'stomate_init'
2571    ENDIF
2572
2573    ALLOCATE(humrel_month(kjpindex,nvm),stat=ier)
2574    l_error = l_error .OR. (ier /= 0)
2575    IF (l_error) THEN
2576       WRITE(numout,*) 'Memory allocation error for humrel_month. We stop. We need kjpindex*nvm words',kjpindex,nvm
2577       STOP 'stomate_init'
2578    ENDIF
2579
2580    ALLOCATE(humrel_week(kjpindex,nvm),stat=ier)
2581    l_error = l_error .OR. (ier /= 0)
2582    IF (l_error) THEN
2583       WRITE(numout,*) 'Memory allocation error for humrel_week. We stop. We need kjpindex*nvm words',kjpindex,nvm
2584       STOP 'stomate_init'
2585    ENDIF
2586
2587    ALLOCATE(t2m_longterm(kjpindex),stat=ier)
2588    l_error = l_error .OR. (ier /= 0)
2589    IF (l_error) THEN
2590       WRITE(numout,*) 'Memory allocation error for t2m_longterm. We stop. We need kjpindex*nvm words',kjpindex,nvm
2591       STOP 'stomate_init'
2592    ENDIF
2593
2594    ALLOCATE(t2m_month(kjpindex),stat=ier)
2595    l_error = l_error .OR. (ier /= 0)
2596    IF (l_error) THEN
2597       WRITE(numout,*) 'Memory allocation error for t2m_month. We stop. We need kjpindex words',kjpindex
2598       STOP 'stomate_init'
2599    ENDIF
2600
2601    ALLOCATE(Tseason(kjpindex),stat=ier)
2602    l_error = l_error .OR. (ier /= 0)
2603    IF (l_error) THEN
2604       WRITE(numout,*) 'Memory allocation error for Tseason. We stop. We need kjpindex words',kjpindex
2605       STOP 'stomate_init'
2606    ENDIF
2607
2608    ALLOCATE(Tseason_length(kjpindex),stat=ier)
2609    l_error = l_error .OR. (ier /= 0)
2610    IF (l_error) THEN
2611       WRITE(numout,*) 'Memory allocation error for Tseason_length. We stop. We need kjpindex words',kjpindex
2612       STOP 'stomate_init'
2613    ENDIF
2614
2615    ALLOCATE(Tseason_tmp(kjpindex),stat=ier)
2616    l_error = l_error .OR. (ier /= 0)
2617    IF (l_error) THEN
2618       WRITE(numout,*) 'Memory allocation error for Tseason_tmp. We stop. We need kjpindex words',kjpindex
2619       STOP 'stomate_init'
2620    ENDIF
2621
2622    ALLOCATE(Tmin_spring_time(kjpindex,nvm),stat=ier)
2623    l_error = l_error .OR. (ier /= 0)
2624    IF (l_error) THEN
2625       WRITE(numout,*) 'Memory allocation error for Tmin_spring_time. We stop. We need kjpindex*nvm words',kjpindex,nvm
2626       STOP 'stomate_init'
2627    ENDIF
2628
2629    ALLOCATE(onset_date(kjpindex,nvm),stat=ier)
2630    l_error = l_error .OR. (ier /= 0)
2631    IF (l_error) THEN
2632       WRITE(numout,*) 'Memory allocation error for onset_date. We stop. We need kjpindex*nvm*nparts words',kjpindex,nvm,2
2633       STOP 'stomate_init'
2634    ENDIF
2635
2636    ALLOCATE(t2m_week(kjpindex),stat=ier)
2637    l_error = l_error .OR. (ier /= 0)
2638    IF (l_error) THEN
2639       WRITE(numout,*) 'Memory allocation error for t2m_week. We stop. We need kjpindex words',kjpindex
2640       STOP 'stomate_init'
2641    ENDIF
2642
2643    ALLOCATE(tsoil_month(kjpindex,nbdl),stat=ier)
2644    l_error = l_error .OR. (ier /= 0)
2645    IF (l_error) THEN
2646       WRITE(numout,*) 'Memory allocation error for tsoil_month. We stop. We need kjpindex*nbdl words',kjpindex,nbdl
2647       STOP 'stomate_init'
2648    ENDIF
2649
2650    ALLOCATE(soilhum_month(kjpindex,nbdl),stat=ier)
2651    l_error = l_error .OR. (ier /= 0)
2652    IF (l_error) THEN
2653       WRITE(numout,*) 'Memory allocation error for soilhum_month. We stop. We need kjpindex*nbdl words',kjpindex,nbdl
2654       STOP 'stomate_init'
2655    ENDIF
2656
2657    ALLOCATE(fireindex(kjpindex,nvm),stat=ier) 
2658    l_error = l_error .OR. (ier /= 0)
2659    IF (l_error) THEN
2660       WRITE(numout,*) 'Memory allocation error for fireindex. We stop. We need kjpindex*nvm words',kjpindex,nvm
2661       STOP 'stomate_init'
2662    ENDIF
2663
2664    ALLOCATE(firelitter(kjpindex,nvm),stat=ier)
2665    l_error = l_error .OR. (ier /= 0)
2666    IF (l_error) THEN
2667       WRITE(numout,*) 'Memory allocation error for firelitter. We stop. We need kjpindex*nvm words',kjpindex,nvm
2668       STOP 'stomate_init'
2669    ENDIF
2670
2671    ALLOCATE(maxhumrel_lastyear(kjpindex,nvm),stat=ier)
2672    l_error = l_error .OR. (ier /= 0)
2673    IF (l_error) THEN
2674       WRITE(numout,*) 'Memory allocation error for maxhumrel_lastyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2675       STOP 'stomate_init'
2676    ENDIF
2677
2678    ALLOCATE(maxhumrel_thisyear(kjpindex,nvm),stat=ier)
2679    l_error = l_error .OR. (ier /= 0)
2680    IF (l_error) THEN
2681       WRITE(numout,*) 'Memory allocation error for maxhumrel_thisyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2682       STOP 'stomate_init'
2683    ENDIF
2684
2685    ALLOCATE(minhumrel_lastyear(kjpindex,nvm),stat=ier)
2686    l_error = l_error .OR. (ier /= 0)
2687    IF (l_error) THEN
2688       WRITE(numout,*) 'Memory allocation error for minhumrel_lastyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2689       STOP 'stomate_init'
2690    ENDIF
2691
2692    ALLOCATE(minhumrel_thisyear(kjpindex,nvm),stat=ier)
2693    l_error = l_error .OR. (ier /= 0)
2694    IF (l_error) THEN
2695       WRITE(numout,*) 'Memory allocation error for minhumrel_thisyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2696       STOP 'stomate_init'
2697    ENDIF
2698
2699    ALLOCATE(maxgppweek_lastyear(kjpindex,nvm),stat=ier)
2700    l_error = l_error .OR. (ier /= 0)
2701    IF (l_error) THEN
2702       WRITE(numout,*) 'Memory allocation error for maxgppweek_lastyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2703       STOP 'stomate_init'
2704    ENDIF
2705
2706    ALLOCATE(maxgppweek_thisyear(kjpindex,nvm),stat=ier)
2707    l_error = l_error .OR. (ier /= 0)
2708    IF (l_error) THEN
2709       WRITE(numout,*) 'Memory allocation error for maxgppweek_thisyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2710       STOP 'stomate_init'
2711    ENDIF
2712
2713    ALLOCATE(gdd0_lastyear(kjpindex),stat=ier)
2714    l_error = l_error .OR. (ier /= 0)
2715    IF (l_error) THEN
2716       WRITE(numout,*) 'Memory allocation error for gdd0_lastyear. We stop. We need kjpindex words',kjpindex
2717       STOP 'stomate_init'
2718    ENDIF
2719
2720    ALLOCATE(gdd0_thisyear(kjpindex),stat=ier)
2721    l_error = l_error .OR. (ier /= 0)
2722    IF (l_error) THEN
2723       WRITE(numout,*) 'Memory allocation error for gdd0_thisyear. We stop. We need kjpindex words',kjpindex
2724       STOP 'stomate_init'
2725    ENDIF
2726
2727    ALLOCATE(gdd_init_date(kjpindex,2),stat=ier)
2728    l_error = l_error .OR. (ier /= 0)
2729    IF (l_error) THEN
2730       WRITE(numout,*) 'Memory allocation error for gdd_init_date. We stop. We need kjpindex*2 words',kjpindex,2
2731       STOP 'stomate_init'
2732    ENDIF
2733
2734    ALLOCATE(gdd_from_growthinit(kjpindex,nvm),stat=ier)
2735    l_error = l_error .OR. (ier /= 0)
2736    IF (l_error) THEN
2737       WRITE(numout,*) 'Memory allocation error for gdd_from_growthinit. We stop. We need kjpindex*nvm words',kjpindex,nvm
2738       STOP 'stomate_init'
2739    ENDIF
2740
2741    ALLOCATE(precip_lastyear(kjpindex),stat=ier)
2742    l_error = l_error .OR. (ier /= 0)
2743    IF (l_error) THEN
2744       WRITE(numout,*) 'Memory allocation error for precip_lastyear. We stop. We need kjpindex*nvm words',kjpindex
2745       STOP 'stomate_init'
2746    ENDIF
2747
2748    ALLOCATE(precip_thisyear(kjpindex),stat=ier)
2749    l_error = l_error .OR. (ier /= 0)
2750    IF (l_error) THEN
2751       WRITE(numout,*) 'Memory allocation error for precip_thisyear. We stop. We need kjpindex words',kjpindex
2752       STOP 'stomate_init'
2753    ENDIF
2754
2755    ALLOCATE(gdd_m5_dormance(kjpindex,nvm),stat=ier)
2756    l_error = l_error .OR. (ier /= 0)
2757    IF (l_error) THEN
2758       WRITE(numout,*) 'Memory allocation error for gdd_m5_dormance. We stop. We need kjpindex*nvm words',kjpindex,nvm
2759       STOP 'stomate_init'
2760    ENDIF
2761
2762    ALLOCATE(gdd_midwinter(kjpindex,nvm),stat=ier)
2763    l_error = l_error .OR. (ier /= 0)
2764    IF (l_error) THEN
2765       WRITE(numout,*) 'Memory allocation error for gdd_midwinter. We stop. We need kjpindex*nvm words',kjpindex,nvm
2766       STOP 'stomate_init'
2767    ENDIF
2768
2769    ALLOCATE(ncd_dormance(kjpindex,nvm),stat=ier)
2770    l_error = l_error .OR. (ier /= 0)
2771    IF (l_error) THEN
2772       WRITE(numout,*) 'Memory allocation error for ncd_dormance. We stop. We need kjpindex*nvm words',kjpindex,nvm
2773       STOP 'stomate_init'
2774    ENDIF
2775
2776    ALLOCATE(ngd_minus5(kjpindex,nvm),stat=ier)
2777    l_error = l_error .OR. (ier /= 0)
2778    IF (l_error) THEN
2779       WRITE(numout,*) 'Memory allocation error for ngd_minus5. We stop. We need kjpindex*nvm words',kjpindex,nvm
2780       STOP 'stomate_init'
2781    ENDIF
2782
2783    ALLOCATE(PFTpresent(kjpindex,nvm),stat=ier)
2784    l_error = l_error .OR. (ier /= 0)
2785    IF (l_error) THEN
2786       WRITE(numout,*) 'Memory allocation error for PFTpresent. We stop. We need kjpindex*nvm words',kjpindex,nvm
2787       STOP 'stomate_init'
2788    ENDIF
2789
2790    ALLOCATE(npp_longterm(kjpindex,nvm),stat=ier)
2791    l_error = l_error .OR. (ier /= 0)
2792    IF (l_error) THEN
2793       WRITE(numout,*) 'Memory allocation error for npp_longterm. We stop. We need kjpindex*nvm words',kjpindex,nvm
2794       STOP 'stomate_init'
2795    ENDIF
2796
2797    ALLOCATE(lm_lastyearmax(kjpindex,nvm),stat=ier)
2798    l_error = l_error .OR. (ier /= 0)
2799    IF (l_error) THEN
2800       WRITE(numout,*) 'Memory allocation error for lm_lastyearmax. We stop. We need kjpindex*nvm words',kjpindex,nvm
2801       STOP 'stomate_init'
2802    ENDIF
2803
2804    ALLOCATE(lm_thisyearmax(kjpindex,nvm),stat=ier)
2805    l_error = l_error .OR. (ier /= 0)
2806    IF (l_error) THEN
2807       WRITE(numout,*) 'Memory allocation error for lm_thisyearmax. We stop. We need kjpindex*nvm words',kjpindex,nvm
2808       STOP 'stomate_init'
2809    ENDIF
2810
2811    ALLOCATE(maxfpc_lastyear(kjpindex,nvm),stat=ier)
2812    l_error = l_error .OR. (ier /= 0)
2813    IF (l_error) THEN
2814       WRITE(numout,*) 'Memory allocation error for maxfpc_lastyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2815       STOP 'stomate_init'
2816    ENDIF
2817
2818    ALLOCATE(maxfpc_thisyear(kjpindex,nvm),stat=ier)
2819    l_error = l_error .OR. (ier /= 0)
2820    IF (l_error) THEN
2821       WRITE(numout,*) 'Memory allocation error for maxfpc_thisyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2822       STOP 'stomate_init'
2823    ENDIF
2824
2825    ALLOCATE(turnover_longterm(kjpindex,nvm,nparts,nelements),stat=ier)
2826    l_error = l_error .OR. (ier /= 0)
2827    IF (l_error) THEN
2828       WRITE(numout,*) 'Memory allocation error for turnover_longterm. We stop. We need kjpindex*nvm*nparts*nelements words', & 
2829       &    kjpindex,nvm,nparts,nelements
2830       STOP 'stomate_init'
2831    ENDIF
2832
2833    ALLOCATE(gpp_week(kjpindex,nvm),stat=ier)
2834    l_error = l_error .OR. (ier /= 0)
2835    IF (l_error) THEN
2836       WRITE(numout,*) 'Memory allocation error for gpp_week. We stop. We need kjpindex*nvm words',kjpindex,nvm
2837       STOP 'stomate_init'
2838    ENDIF
2839
2840    ALLOCATE(biomass(kjpindex,nvm,nparts,nelements),stat=ier)
2841    l_error = l_error .OR. (ier /= 0)
2842    IF (l_error) THEN
2843       WRITE(numout,*) 'Memory allocation error for biomass. We stop. We need kjpindex*nvm*nparts*nelements words', &
2844       &    kjpindex,nvm,nparts,nelements
2845       STOP 'stomate_init'
2846    ENDIF
2847
2848    ALLOCATE(senescence(kjpindex,nvm),stat=ier)
2849    l_error = l_error .OR. (ier /= 0)
2850    IF (l_error) THEN
2851       WRITE(numout,*) 'Memory allocation error for senescence. We stop. We need kjpindex*nvm words',kjpindex,nvm
2852       STOP 'stomate_init'
2853    ENDIF
2854
2855    ALLOCATE(begin_leaves(kjpindex,nvm),stat=ier)
2856    l_error = l_error .OR. (ier /= 0)
2857    IF (l_error) THEN
2858       WRITE(numout,*) 'Memory allocation error for begin_leaves. We stop. We need kjpindex*nvm words',kjpindex,nvm
2859       STOP 'stomate_init'
2860    ENDIF
2861
2862    ALLOCATE(when_growthinit(kjpindex,nvm),stat=ier)
2863    l_error = l_error .OR. (ier /= 0)
2864    IF (l_error) THEN
2865       WRITE(numout,*) 'Memory allocation error for when_growthinit. We stop. We need kjpindex*nvm words',kjpindex,nvm
2866       STOP 'stomate_init'
2867    ENDIF
2868
2869    ALLOCATE(age(kjpindex,nvm),stat=ier)
2870    l_error = l_error .OR. (ier /= 0)
2871    IF (l_error) THEN
2872       WRITE(numout,*) 'Memory allocation error for age. We stop. We need kjpindex*nvm words',kjpindex,nvm
2873       STOP 'stomate_init'
2874    ENDIF
2875
2876    ALLOCATE(resp_hetero_d(kjpindex,nvm),stat=ier)
2877    l_error = l_error .OR. (ier /= 0)
2878    IF (l_error) THEN
2879       WRITE(numout,*) 'Memory allocation error for resp_hetero_d. We stop. We need kjpindex*nvm words',kjpindex,nvm
2880       STOP 'stomate_init'
2881    ENDIF
2882
2883    ALLOCATE(resp_hetero_radia(kjpindex,nvm),stat=ier)
2884    l_error = l_error .OR. (ier /= 0)
2885    IF (l_error) THEN
2886       WRITE(numout,*) 'Memory allocation error for resp_hetero_radia. We stop. We need kjpindex*nvm words',kjpindex,nvm
2887       STOP 'stomate_init'
2888    ENDIF
2889
2890    ALLOCATE(resp_maint_d(kjpindex,nvm),stat=ier)
2891    l_error = l_error .OR. (ier /= 0)
2892    IF (l_error) THEN
2893       WRITE(numout,*) 'Memory allocation error for resp_maint_d. We stop. We need kjpindex*nvm words',kjpindex,nvm
2894       STOP 'stomate_init'
2895    ENDIF
2896
2897    ALLOCATE(resp_growth_d(kjpindex,nvm),stat=ier)
2898    l_error = l_error .OR. (ier /= 0)
2899    IF (l_error) THEN
2900       WRITE(numout,*) 'Memory allocation error for resp_growth_d. We stop. We need kjpindex*nvm words',kjpindex,nvm
2901       STOP 'stomate_init'
2902    ENDIF
2903
2904    ALLOCATE(co2_fire(kjpindex,nvm),stat=ier)
2905    l_error = l_error .OR. (ier /= 0)
2906    IF (l_error) THEN
2907       WRITE(numout,*) 'Memory allocation error for co2_fire. We stop. We need kjpindex*nvm words',kjpindex,nvm
2908       STOP 'stomate_init'
2909    ENDIF
2910
2911    ALLOCATE(co2_to_bm_dgvm(kjpindex,nvm),stat=ier)
2912    l_error = l_error .OR. (ier /= 0)
2913    IF (l_error) THEN
2914       WRITE(numout,*) 'Memory allocation error for co2_to_bm_dgvm. We stop. We need kjpindex*nvm words',kjpindex,nvm
2915       STOP 'stomate_init'
2916    ENDIF
2917
2918    ALLOCATE(veget_lastlight(kjpindex,nvm),stat=ier)
2919    l_error = l_error .OR. (ier /= 0)
2920    IF (l_error) THEN
2921       WRITE(numout,*) 'Memory allocation error for veget_lastlight. We stop. We need kjpindex*nvm words',kjpindex,nvm
2922       STOP 'stomate_init'
2923    ENDIF
2924
2925    ALLOCATE(everywhere(kjpindex,nvm),stat=ier)
2926    l_error = l_error .OR. (ier /= 0)
2927    IF (l_error) THEN
2928       WRITE(numout,*) 'Memory allocation error for everywhere. We stop. We need kjpindex*nvm words',kjpindex,nvm
2929       STOP 'stomate_init'
2930    ENDIF
2931
2932    ALLOCATE(need_adjacent(kjpindex,nvm),stat=ier)
2933    l_error = l_error .OR. (ier /= 0)
2934    IF (l_error) THEN
2935       WRITE(numout,*) 'Memory allocation error for need_adjacent. We stop. We need kjpindex*nvm words',kjpindex,nvm
2936       STOP 'stomate_init'
2937    ENDIF
2938
2939    ALLOCATE(leaf_age(kjpindex,nvm,nleafages),stat=ier)
2940    l_error = l_error .OR. (ier /= 0)
2941    IF (l_error) THEN
2942       WRITE(numout,*) 'Memory allocation error for leaf_age. We stop. We need kjpindex*nvm*nleafages words', & 
2943       &      kjpindex,nvm,nleafages
2944       STOP 'stomate_init'
2945    ENDIF
2946
2947    ALLOCATE(leaf_frac(kjpindex,nvm,nleafages),stat=ier)
2948    l_error = l_error .OR. (ier /= 0)
2949    IF (l_error) THEN
2950       WRITE(numout,*) 'Memory allocation error for leaf_frac. We stop. We need kjpindex*nvm*nleafages words', & 
2951       &      kjpindex,nvm,nleafages
2952       STOP 'stomate_init'
2953    ENDIF
2954
2955    ALLOCATE(RIP_time(kjpindex,nvm),stat=ier)
2956    l_error = l_error .OR. (ier /= 0)
2957    IF (l_error) THEN
2958       WRITE(numout,*) 'Memory allocation error for RIP_time. We stop. We need kjpindex*nvm words',kjpindex,nvm
2959       STOP 'stomate_init'
2960    ENDIF
2961
2962    ALLOCATE(time_hum_min(kjpindex,nvm),stat=ier)
2963    l_error = l_error .OR. (ier /= 0)
2964    IF (l_error) THEN
2965       WRITE(numout,*) 'Memory allocation error for time_hum_min. We stop. We need kjpindex*nvm words',kjpindex,nvm
2966       STOP 'stomate_init'
2967    ENDIF
2968
2969    ALLOCATE(hum_min_dormance(kjpindex,nvm),stat=ier)
2970    l_error = l_error .OR. (ier /= 0)
2971    IF (l_error) THEN
2972       WRITE(numout,*) 'Memory allocation error for hum_min_dormance. We stop. We need kjpindex*nvm words',kjpindex,nvm
2973       STOP 'stomate_init'
2974    ENDIF
2975
2976    ALLOCATE(litterpart(kjpindex,nvm,nlitt),stat=ier)
2977    l_error = l_error .OR. (ier /= 0)
2978    IF (l_error) THEN
2979       WRITE(numout,*) 'Memory allocation error for litterpart. We stop. We need kjpindex*nvm*nlitt words',  &
2980       &  kjpindex,nvm,nlitt
2981       STOP 'stomate_init'
2982    ENDIF
2983
2984    ALLOCATE(litter(kjpindex,nlitt,nvm,nlevs,nelements),stat=ier)
2985    l_error = l_error .OR. (ier /= 0)
2986    IF (l_error) THEN
2987       WRITE(numout,*) 'Memory allocation error for litter. We stop. We need kjpindex*nlitt*nvm*nlevs*nelements words', & 
2988       &    kjpindex,nlitt,nvm,nlevs,nelements
2989       STOP 'stomate_init'
2990    ENDIF
2991
2992    ALLOCATE(dead_leaves(kjpindex,nvm,nlitt),stat=ier)
2993    l_error = l_error .OR. (ier /= 0)
2994    IF (l_error) THEN
2995       WRITE(numout,*) 'Memory allocation error for dead_leaves. We stop. We need kjpindex*nvm*nlitt words', & 
2996       &   kjpindex,nvm,nlitt
2997       STOP 'stomate_init'
2998    ENDIF
2999
3000    ALLOCATE(carbon(kjpindex,ncarb,nvm),stat=ier)
3001    l_error = l_error .OR. (ier /= 0)
3002    IF (l_error) THEN
3003       WRITE(numout,*) 'Memory allocation error for carbon. We stop. We need kjpindex*ncarb*nvm words',kjpindex,ncarb,nvm
3004       STOP 'stomate_init'
3005    ENDIF
3006
3007    ALLOCATE(lignin_struc(kjpindex,nvm,nlevs),stat=ier)
3008    l_error = l_error .OR. (ier /= 0)
3009    IF (l_error) THEN
3010       WRITE(numout,*) 'Memory allocation error for lignin_struc. We stop. We need kjpindex*nvm*nlevs words',kjpindex,nvm,nlevs
3011       STOP 'stomate_init'
3012    ENDIF
3013
3014    ALLOCATE(turnover_time(kjpindex,nvm),stat=ier)
3015    l_error = l_error .OR. (ier /= 0)
3016    IF (l_error) THEN
3017       WRITE(numout,*) 'Memory allocation error for turnover_time. We stop. We need kjpindex*nvm words',kjpindex,nvm
3018       STOP 'stomate_init'
3019    ENDIF
3020
3021    ALLOCATE(co2_flux_daily(kjpindex,nvm),stat=ier)
3022    l_error = l_error .OR. (ier /= 0)
3023    IF (l_error) THEN
3024       WRITE(numout,*) 'Memory allocation error for co2_flux_daily. We stop. We need kjpindex*nvm words',kjpindex,nvm
3025       STOP 'stomate_init'
3026    ENDIF
3027
3028    ALLOCATE(co2_flux_monthly(kjpindex,nvm),stat=ier)
3029    l_error = l_error .OR. (ier /= 0)
3030    IF (l_error) THEN
3031       WRITE(numout,*) 'Memory allocation error for co2_flux_monthly. We stop. We need kjpindex*nvm words',kjpindex,nvm
3032       STOP 'stomate_init'
3033    ENDIF
3034
3035    ALLOCATE (cflux_prod_monthly(kjpindex), stat=ier)
3036    l_error = l_error .OR. (ier /= 0)
3037    IF (l_error) THEN
3038       WRITE(numout,*) 'Memory allocation error for cflux_prod_monthly. We stop. We need kjpindex words',kjpindex
3039       STOP 'stomate_init'
3040    ENDIF
3041 
3042    ALLOCATE (harvest_above_monthly(kjpindex), stat=ier)
3043    l_error = l_error .OR. (ier /= 0)
3044    IF (l_error) THEN
3045       WRITE(numout,*) 'Memory allocation error for harvest_above_monthly. We stop. We need kjpindex words',kjpindex
3046       STOP 'stomate_init'
3047    ENDIF
3048
3049    ALLOCATE(bm_to_litter(kjpindex,nvm,nparts,nelements),stat=ier)
3050    l_error = l_error .OR. (ier /= 0)
3051    IF (l_error) THEN
3052       WRITE(numout,*) 'Memory allocation error for bm_to_litter. We stop. We need kjpindex*nvm*nparts*nelements words', & 
3053       &    kjpindex,nvm,nparts,nelements
3054       STOP 'stomate_init'
3055    ENDIF
3056
3057    ALLOCATE(bm_to_littercalc(kjpindex,nvm,nparts,nelements),stat=ier)
3058    l_error = l_error .OR. (ier /= 0)
3059    IF (l_error) THEN
3060       WRITE(numout,*) 'Memory allocation error for bm_to_littercalc. We stop. We need kjpindex*nvm*nparts*nelements words', &
3061       &   kjpindex,nvm,nparts,nelements
3062       STOP 'stomate_init'
3063    ENDIF
3064
3065    ALLOCATE(herbivores(kjpindex,nvm),stat=ier)
3066    l_error = l_error .OR. (ier /= 0)
3067    IF (l_error) THEN
3068       WRITE(numout,*) 'Memory allocation error for herbivores. We stop. We need kjpindex*nvm words',kjpindex,nvm
3069       STOP 'stomate_init'
3070    ENDIF
3071
3072    ALLOCATE(hori_index(kjpindex),stat=ier)
3073    l_error = l_error .OR. (ier /= 0)
3074    IF (l_error) THEN
3075       WRITE(numout,*) 'Memory allocation error for hori_index. We stop. We need kjpindex words',kjpindex
3076       STOP 'stomate_init'
3077    ENDIF
3078
3079    ALLOCATE(horipft_index(kjpindex*nvm),stat=ier)
3080    l_error = l_error .OR. (ier /= 0)
3081    IF (l_error) THEN
3082       WRITE(numout,*) 'Memory allocation error for horipft_index. We stop. We need kjpindex*nvm words',kjpindex*nvm
3083       STOP 'stomate_init'
3084    ENDIF
3085
3086    ALLOCATE(resp_maint_part_radia(kjpindex,nvm,nparts),stat=ier)
3087    l_error = l_error .OR. (ier /= 0)
3088    IF (l_error) THEN
3089       WRITE(numout,*) 'Memory allocation error for resp_maint_part_radia. We stop. We need kjpindex*nvm*nparts words', &
3090       &  kjpindex,nvm,nparts
3091       STOP 'stomate_init'
3092    ENDIF
3093
3094    ALLOCATE(resp_maint_radia(kjpindex,nvm),stat=ier)
3095    l_error = l_error .OR. (ier /= 0)
3096    IF (l_error) THEN
3097       WRITE(numout,*) 'Memory allocation error for resp_maint_radia. We stop. We need kjpindex*nvm words',kjpindex,nvm
3098       STOP 'stomate_init'
3099    ENDIF
3100
3101    ALLOCATE(resp_maint_part(kjpindex,nvm,nparts),stat=ier)
3102    l_error = l_error .OR. (ier /= 0)
3103    IF (l_error) THEN
3104       WRITE(numout,*) 'Memory allocation error for resp_maint_part. We stop. We need kjpindex*nvm*nparts words', &
3105       &    kjpindex,nvm,nparts
3106       STOP 'stomate_init'
3107    ENDIF
3108    resp_maint_part(:,:,:) = zero
3109
3110    ALLOCATE (horip10_index(kjpindex*10), stat=ier)
3111    l_error = l_error .OR. (ier /= 0)
3112    IF (l_error) THEN
3113       WRITE(numout,*) 'Memory allocation error for horip10_index. We stop. We need kjpindex*10 words',kjpindex,10
3114       STOP 'stomate_init'
3115    ENDIF
3116
3117    ALLOCATE (horip100_index(kjpindex*100), stat=ier)
3118    l_error = l_error .OR. (ier /= 0)
3119    IF (l_error) THEN
3120       WRITE(numout,*) 'Memory allocation error for horip100_index. We stop. We need kjpindex*100 words',kjpindex,100
3121       STOP 'stomate_init'
3122    ENDIF
3123
3124    ALLOCATE (horip11_index(kjpindex*11), stat=ier)
3125    l_error = l_error .OR. (ier /= 0)
3126    IF (l_error) THEN
3127       WRITE(numout,*) 'Memory allocation error for horip11_index. We stop. We need kjpindex*11 words',kjpindex,11
3128       STOP 'stomate_init'
3129    ENDIF
3130
3131    ALLOCATE (horip101_index(kjpindex*101), stat=ier)
3132    l_error = l_error .OR. (ier /= 0)
3133    IF (l_error) THEN
3134       WRITE(numout,*) 'Memory allocation error for horip101_index. We stop. We need kjpindex*101 words',kjpindex,101
3135       STOP 'stomate_init'
3136    ENDIF
3137
3138    ALLOCATE (prod10(kjpindex,0:10), stat=ier)
3139    l_error = l_error .OR. (ier /= 0)
3140    IF (l_error) THEN
3141       WRITE(numout,*) 'Memory allocation error for prod10. We stop. We need kjpindex*11 words',kjpindex,11
3142       STOP 'stomate_init'
3143    ENDIF
3144
3145    ALLOCATE (prod100(kjpindex,0:100), stat=ier)
3146    l_error = l_error .OR. (ier /= 0)
3147    IF (l_error) THEN
3148       WRITE(numout,*) 'Memory allocation error for prod100. We stop. We need kjpindex*101 words',kjpindex,101
3149       STOP 'stomate_init'
3150    ENDIF
3151
3152    ALLOCATE (flux10(kjpindex,10), stat=ier)
3153    l_error = l_error .OR. (ier /= 0)
3154    IF (l_error) THEN
3155       WRITE(numout,*) 'Memory allocation error for flux10. We stop. We need kjpindex*10 words',kjpindex,10
3156       STOP 'stomate_init'
3157    ENDIF
3158
3159    ALLOCATE (flux100(kjpindex,100), stat=ier)
3160    l_error = l_error .OR. (ier /= 0)
3161    IF (l_error) THEN
3162       WRITE(numout,*) 'Memory allocation error for flux100. We stop. We need kjpindex*100 words',kjpindex,100
3163       STOP 'stomate_init'
3164    ENDIF
3165
3166    ALLOCATE (convflux(kjpindex), stat=ier)
3167    l_error = l_error .OR. (ier /= 0)
3168    IF (l_error) THEN
3169       WRITE(numout,*) 'Memory allocation error for convflux. We stop. We need kjpindex words',kjpindex
3170       STOP 'stomate_init'
3171    ENDIF
3172
3173    ALLOCATE (cflux_prod10(kjpindex), stat=ier)
3174    l_error = l_error .OR. (ier /= 0)
3175    IF (l_error) THEN
3176       WRITE(numout,*) 'Memory allocation error for cflux_prod10. We stop. We need kjpindex words',kjpindex
3177       STOP 'stomate_init'
3178    ENDIF
3179
3180    ALLOCATE (cflux_prod100(kjpindex), stat=ier)
3181    l_error = l_error .OR. (ier /= 0)
3182    IF (l_error) THEN
3183       WRITE(numout,*) 'Memory allocation error for cflux_prod100. We stop. We need kjpindex words',kjpindex
3184       STOP 'stomate_init'
3185    ENDIF
3186
3187    ALLOCATE (harvest_above(kjpindex), stat=ier)
3188    l_error = l_error .OR. (ier /= 0)
3189    IF (l_error) THEN
3190       WRITE(numout,*) 'Memory allocation error for harvest_above. We stop. We need kjpindex words',kjpindex
3191       STOP 'stomate_init'
3192    ENDIF
3193
3194    ALLOCATE (carb_mass_total(kjpindex), stat=ier)
3195    l_error = l_error .OR. (ier /= 0)
3196    IF (l_error) THEN
3197       WRITE(numout,*) 'Memory allocation error for carb_mass_total. We stop. We need kjpindex words',kjpindex
3198       STOP 'stomate_init'
3199    ENDIF
3200
3201    ALLOCATE (soilcarbon_input_daily(kjpindex,ncarb,nvm), stat=ier)
3202    l_error = l_error .OR. (ier /= 0)
3203    IF (l_error) THEN
3204       WRITE(numout,*) 'Memory allocation error for soilcarbon_input_daily. We stop. We need kjpindex*ncarb*nvm words', & 
3205       &    kjpindex,ncarb,nvm
3206       STOP 'stomate_init'
3207    ENDIF
3208
3209    ALLOCATE (control_temp_daily(kjpindex,nlevs), stat=ier)
3210    l_error = l_error .OR. (ier /= 0)
3211    IF (l_error) THEN
3212       WRITE(numout,*) 'Memory allocation error for control_temp_daily. We stop. We need kjpindex*nlevs words',kjpindex,nlevs
3213       STOP 'stomate_init'
3214    ENDIF
3215
3216    ALLOCATE (control_moist_daily(kjpindex,nlevs), stat=ier)
3217    l_error = l_error .OR. (ier /= 0)
3218    IF (l_error) THEN
3219       WRITE(numout,*) 'Memory allocation error for control_moist_daily. We stop. We need kjpindex*nlevs words',kjpindex,nlevs
3220       STOP 'stomate_init'
3221    ENDIF
3222
3223    ALLOCATE (fpc_max(kjpindex,nvm), stat=ier)
3224    l_error = l_error .OR. (ier /= 0)
3225    IF (l_error) THEN
3226       WRITE(numout,*) 'Memory allocation error for fpc_max. We stop. We need kjpindex*nvm words',kjpindex,nvm
3227       STOP 'stomate_init'
3228    ENDIF
3229
3230    ALLOCATE(ok_equilibrium(kjpindex),stat=ier)
3231    l_error = l_error .OR. (ier /= 0) 
3232    IF (l_error) THEN
3233       WRITE(numout,*) 'Memory allocation error for ok_equilibrium. We stop. We need kjpindex words',kjpindex
3234       STOP 'stomate_init'
3235    ENDIF
3236
3237    ALLOCATE(carbon_eq(kjpindex),stat=ier)
3238    l_error = l_error .OR. (ier /= 0)
3239    IF (l_error) THEN
3240       WRITE(numout,*) 'Memory allocation error for carbon_eq. We stop. We need kjpindex words',kjpindex
3241       STOP 'stomate_init'
3242    ENDIF
3243
3244    ALLOCATE(nbp_accu(kjpindex),stat=ier)
3245    l_error = l_error .OR. (ier /= 0)
3246    IF (l_error) THEN
3247       WRITE(numout,*) 'Memory allocation error for nbp_accu. We stop. We need kjpindex words',kjpindex
3248       STOP 'stomate_init'
3249    ENDIF
3250
3251    ALLOCATE(nbp_flux(kjpindex),stat=ier)
3252    l_error = l_error .OR. (ier /= 0)
3253    IF (l_error) THEN
3254       WRITE(numout,*) 'Memory allocation error for nbp_flux. We stop. We need kjpindex words',kjpindex
3255       STOP 'stomate_init'
3256    ENDIF
3257
3258    ALLOCATE(matrixA(kjpindex,nvm,nbpools,nbpools),stat=ier)
3259    l_error = l_error .OR. (ier /= 0)
3260    IF (l_error) THEN
3261       WRITE(numout,*) 'Memory allocation error for matrixA. We stop. We need kjpindex*nvm*nbpools*nbpools words',  & 
3262       &     kjpindex, nvm, nbpools, nbpools
3263       STOP 'stomate_init'
3264    ENDIF
3265
3266    ALLOCATE(vectorB(kjpindex,nvm,nbpools),stat=ier)
3267    l_error = l_error .OR. (ier /= 0)
3268    IF (l_error) THEN
3269       WRITE(numout,*) 'Memory allocation error for vectorB. We stop. We need kjpindex*nvm*nbpools words',  & 
3270       &     kjpindex, nvm, nbpools
3271       STOP 'stomate_init'
3272    ENDIF
3273
3274    ALLOCATE(VectorU(kjpindex,nvm,nbpools),stat=ier)
3275    l_error = l_error .OR. (ier /= 0)
3276    IF (l_error) THEN
3277       WRITE(numout,*) 'Memory allocation error for VectorU. We stop. We need kjpindex*nvm*nbpools words',  & 
3278       &     kjpindex, nvm, nbpools
3279       STOP 'stomate_init'
3280    ENDIF
3281
3282    ALLOCATE(MatrixV(kjpindex,nvm,nbpools,nbpools),stat=ier)
3283    l_error = l_error .OR. (ier /= 0)
3284    IF (l_error) THEN
3285       WRITE(numout,*) 'Memory allocation error for MatrixV. We stop. We need kjpindex*nvm*nbpools*nbpools words',  & 
3286       &     kjpindex, nvm, nbpools, nbpools
3287       STOP 'stomate_init'
3288    ENDIF
3289
3290    ALLOCATE(MatrixW(kjpindex,nvm,nbpools,nbpools),stat=ier)
3291    l_error = l_error .OR. (ier /= 0)
3292    IF (l_error) THEN
3293       WRITE(numout,*) 'Memory allocation error for MatrixW. We stop. We need kjpindex*nvm*nbpools*nbpools words',  & 
3294       &     kjpindex, nvm, nbpools, nbpools
3295       STOP 'stomate_init'
3296    ENDIF
3297
3298    ALLOCATE(previous_stock(kjpindex,nvm,nbpools),stat=ier)
3299    l_error = l_error .OR. (ier /= 0)
3300    IF (l_error) THEN
3301       WRITE(numout,*) 'Memory allocation error for previous_stock. We stop. We need kjpindex*nvm*nbpools words',  & 
3302       &     kjpindex, nvm, nbpools
3303       STOP 'stomate_init'
3304    ENDIF
3305
3306    ALLOCATE(current_stock(kjpindex,nvm,nbpools),stat=ier)
3307    l_error = l_error .OR. (ier /= 0)
3308    IF (l_error) THEN
3309       WRITE(numout,*) 'Memory allocation error for current_stock. We stop. We need kjpindex*nvm*nbpools words',  & 
3310       &     kjpindex, nvm, nbpools
3311       STOP 'stomate_init'
3312    ENDIF
3313!gmjc for grazing litter
3314    ALLOCATE(litter_avail(kjpindex,nlitt,nvm),stat=ier)
3315    l_error = l_error .OR. (ier /= 0)
3316    IF (l_error) THEN
3317       WRITE(numout,*) 'Memory allocation error for litter_avail. We stop. We need kjpindex*nlitt*nvm words',  &
3318       &  kjpindex,nlitt,nvm
3319       STOP 'stomate_init'
3320    ENDIF
3321
3322    ALLOCATE(litter_not_avail(kjpindex,nlitt,nvm),stat=ier)
3323    l_error = l_error .OR. (ier /= 0)
3324    IF (l_error) THEN
3325       WRITE(numout,*) 'Memory allocation error for litter_not_avail. We stop. We need kjpindex*nlitt*nvm words',  &
3326       &  kjpindex,nlitt,nvm
3327       STOP 'stomate_init'
3328    ENDIF
3329
3330    ALLOCATE(litter_avail_frac(kjpindex,nlitt,nvm),stat=ier)
3331    l_error = l_error .OR. (ier /= 0)
3332    IF (l_error) THEN
3333       WRITE(numout,*) 'Memory allocation error for litter_avail_frac. We stop. We need kjpindex*nlitt*nvm words',  &
3334       &  kjpindex,nlitt,nvm
3335       STOP 'stomate_init'
3336    ENDIF
3337
3338    ALLOCATE(sla_calc(kjpindex,nvm),stat=ier)
3339   l_error = l_error .OR. (ier /= 0)
3340    IF (l_error) THEN
3341       WRITE(numout,*) 'Memory allocation error for sla_calc. We stop. We need kjpindex*nvm words',kjpindex,nvm
3342       STOP 'stomate_init'
3343    ENDIF
3344
3345   ALLOCATE(wshtotsum(kjpindex,nvm),stat=ier)
3346   l_error = l_error .OR. (ier /= 0)
3347    IF (l_error) THEN
3348       WRITE(numout,*) 'Memory allocation error for wshtotsum. We stop. We need kjpindex*nvm words',kjpindex,nvm
3349       STOP 'stomate_init'
3350    ENDIF
3351
3352   ALLOCATE(sr_ugb(kjpindex,nvm),stat=ier)
3353   l_error = l_error .OR. (ier /= 0)
3354    IF (l_error) THEN
3355       WRITE(numout,*) 'Memory allocation error for sr_ugb. We stop. We need kjpindex*nvm words',kjpindex,nvm
3356       STOP 'stomate_init'
3357    ENDIF
3358
3359   ALLOCATE(compt_ugb(kjpindex,nvm),stat=ier)
3360   l_error = l_error .OR. (ier /= 0)
3361    IF (l_error) THEN
3362       WRITE(numout,*) 'Memory allocation error for compt_ugb. We stop. We need kjpindex*nvm words',kjpindex,nvm
3363       STOP 'stomate_init'
3364    ENDIF
3365
3366   ALLOCATE(nb_ani(kjpindex,nvm),stat=ier)
3367   l_error = l_error .OR. (ier /= 0)
3368    IF (l_error) THEN
3369       WRITE(numout,*) 'Memory allocation error for nb_ani. We stop. We need kjpindex*nvm words',kjpindex,nvm
3370       STOP 'stomate_init'
3371    ENDIF
3372
3373   ALLOCATE(grazed_frac(kjpindex,nvm),stat=ier)
3374   l_error = l_error .OR. (ier /= 0)
3375    IF (l_error) THEN
3376       WRITE(numout,*) 'Memory allocation error for grazed_frac. We stop. We need kjpindex*nvm words',kjpindex,nvm
3377       STOP 'stomate_init'
3378    ENDIF
3379
3380   ALLOCATE(import_yield(kjpindex,nvm),stat=ier)
3381   l_error = l_error .OR. (ier /= 0)
3382    IF (l_error) THEN
3383       WRITE(numout,*) 'Memory allocation error for import_yield. We stop. We need kjpindex*nvm words',kjpindex,nvm
3384       STOP 'stomate_init'
3385    ENDIF
3386
3387   ALLOCATE(t2m_14(kjpindex),stat=ier)
3388   l_error = l_error .OR. (ier /= 0)
3389    IF (l_error) THEN
3390       WRITE(numout,*) 'Memory allocation error for t2m_14. We stop. We need kjpindex*nvm words',kjpindex
3391       STOP 'stomate_init'
3392    ENDIF
3393
3394   ALLOCATE (sla_age1(kjpindex,nvm), stat=ier)
3395    l_error = l_error .OR. (ier.NE.0)
3396    IF (l_error) THEN
3397       WRITE(numout,*) 'Memory allocation error for sla_age1. We stop. We need kjpindex*nvm words',kjpindex,nvm
3398       STOP 'stomate_init'
3399    ENDIF
3400
3401   ALLOCATE (N_limfert(kjpindex,nvm), stat=ier)
3402    l_error = l_error .OR. (ier.NE.0)
3403    IF (l_error) THEN
3404       WRITE(numout,*) 'Memory allocation error for N_limfert. We stop. We need kjpindex*nvm words',kjpindex,nvm
3405       STOP 'stomate_init'
3406    ENDIF
3407   ALLOCATE (when_growthinit_cut(kjpindex,nvm), stat=ier)
3408    l_error = l_error .OR. (ier.NE.0)
3409    IF (l_error) THEN
3410       WRITE(numout,*) 'Memory allocation error for when_growthinit_cut. We stop. We need kjpindex*nvm words',kjpindex,nvm
3411       STOP 'stomate_init'
3412    ENDIF
3413
3414   ALLOCATE(nb_grazingdays(kjpindex,nvm),stat=ier)
3415   l_error = l_error .OR. (ier /= 0)
3416    IF (l_error) THEN
3417       WRITE(numout,*) 'Memory allocation error for nb_grazingdays. We stop. We need kjpindex*nvm words',kjpindex,nvm
3418       STOP 'stomate_init'
3419    ENDIF
3420
3421   ALLOCATE(snowfall_daily(kjpindex),stat=ier)
3422   l_error = l_error .OR. (ier /= 0)
3423    IF (l_error) THEN
3424       WRITE(numout,*) 'Memory allocation error for snowfall_daily. We stop. We need kjpindex words',kjpindex
3425       STOP 'stomate_init'
3426    ENDIF
3427
3428   ALLOCATE(snowmass_daily(kjpindex),stat=ier)
3429   l_error = l_error .OR. (ier /= 0)
3430    IF (l_error) THEN
3431       WRITE(numout,*) 'Memory allocation error for snowmass_daily. We stop. We need kjpindex words',kjpindex
3432       STOP 'stomate_init'
3433    ENDIF
3434!JCADD top 5 layer grassland soil moisture for grazing
3435    ALLOCATE(tmc_topgrass_daily(kjpindex),stat=ier)
3436    l_error = l_error .OR. (ier /= 0)
3437    IF (l_error) THEN
3438       WRITE(numout,*) 'Memory allocation error for tmc_topgrass_daily. We stop. We need kjpindex words',kjpindex
3439       STOP 'stomate_init'
3440    ENDIF
3441    ALLOCATE(after_snow(kjpindex),stat=ier)
3442    l_error = l_error .OR. (ier /= 0)
3443    IF (l_error) THEN
3444       WRITE(numout,*) 'Memory allocation error for after_snow. We stop. We need kjpindex words',kjpindex
3445       STOP 'stomate_init'
3446    ENDIF
3447    ALLOCATE(after_wet(kjpindex),stat=ier)
3448    l_error = l_error .OR. (ier /= 0)
3449    IF (l_error) THEN
3450       WRITE(numout,*) 'Memory allocation error for after_wet. We stop. We need kjpindex words',kjpindex
3451       STOP 'stomate_init'
3452    ENDIF
3453    ALLOCATE(wet1day(kjpindex),stat=ier)
3454    l_error = l_error .OR. (ier /= 0)
3455    IF (l_error) THEN
3456       WRITE(numout,*) 'Memory allocation error for wet1day. We stop. We need kjpindex words',kjpindex
3457       STOP 'stomate_init'
3458    ENDIF
3459    ALLOCATE(wet2day(kjpindex),stat=ier)
3460    l_error = l_error .OR. (ier /= 0)
3461    IF (l_error) THEN
3462       WRITE(numout,*) 'Memory allocation error for wet2day. We stop. We need kjpindex words',kjpindex
3463       STOP 'stomate_init'
3464    ENDIF
3465!end gmjc   
3466  !! 5. File definitions
3467
3468    ! Store history and restart files in common variables
3469    hist_id_stomate = hist_id_stom
3470    hist_id_stomate_IPCC = hist_id_stom_IPCC
3471    rest_id_stomate = rest_id_stom
3472   
3473    ! In STOMATE reduced grids are used containing only terrestrial pixels.
3474    ! Build a new indexing table for the vegetation fields separating
3475    ! between the different PFTs. Note that ::index has dimension (kjpindex)
3476    ! wheras ::indexpft has dimension (kjpindex*nvm).
3477
3478    hori_index(:) = index(:)
3479
3480    DO j = 1, nvm
3481       DO ji = 1, kjpindex
3482          horipft_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3483       ENDDO
3484    ENDDO
3485
3486    ! Similar index tables are build for the land cover change variables
3487    DO j = 1, 10
3488       DO ji = 1, kjpindex
3489          horip10_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3490       ENDDO
3491    ENDDO
3492
3493    DO j = 1, 100
3494       DO ji = 1, kjpindex
3495          horip100_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3496       ENDDO
3497    ENDDO
3498
3499    DO j = 1, 11
3500       DO ji = 1, kjpindex
3501          horip11_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3502       ENDDO
3503    ENDDO
3504
3505    DO j = 1, 101
3506       DO ji = 1, kjpindex
3507          horip101_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3508       ENDDO
3509    ENDDO
3510
3511  !! 6. Initialization of global and land cover change variables.
3512
3513    ! All variables are cumulative variables. bm_to_litter is not and is therefore
3514    ! excluded
3515    !   bm_to_litter(:,:,:) = zero
3516    turnover_daily(:,:,:,:) = zero
3517    resp_hetero_d(:,:) = zero
3518    co2_flux_daily(:,:) = zero
3519    co2_flux_monthly(:,:) = zero
3520    cflux_prod_monthly(:) = zero
3521    harvest_above_monthly(:) = zero
3522    control_moist_daily(:,:) = zero
3523    control_temp_daily(:,:) = zero
3524    soilcarbon_input_daily(:,:,:) = zero
3525    ! Land cover change variables
3526    prod10(:,:)  = zero
3527    prod100(:,:) = zero
3528    flux10(:,:)  = zero
3529    flux100(:,:) = zero
3530    convflux(:)  = zero
3531    cflux_prod10(:) = zero
3532    cflux_prod100(:) = zero
3533    fpc_max(:,:)=zero
3534
3535!gmjc
3536   DO j = 1,nvm
3537     sla_calc(:,j) = sla(j)
3538     sla_age1(:,j) = sla_max(j)
3539   ENDDO
3540   wshtotsum(:,:)=zero
3541   sr_ugb(:,:) = zero
3542   compt_ugb(:,:)=zero
3543   nb_ani(:,:)=zero
3544   grazed_frac(:,:)=zero
3545   import_yield(:,:)=zero
3546   N_limfert(:,:) = 1.0
3547   litter_not_avail = 0.0
3548   litter_avail = 0.0
3549   litter_avail_frac = 1.0
3550   when_growthinit_cut(:,:)=20.0
3551   nb_grazingdays(:,:) = zero
3552   snowfall_daily(:) = zero
3553   snowmass_daily(:) = zero
3554!JCADD top 5 layer grassland soil moisture for grazing
3555   tmc_topgrass_daily(:) = zero
3556   after_snow(:) = zero
3557   after_wet(:) = zero
3558   wet1day(:) = 6
3559   wet2day(:) = 6
3560!ENDJCADD
3561!end gmjc
3562   
3563  END SUBROUTINE stomate_init
3564
3565
3566!! ================================================================================================================================
3567!! SUBROUTINE   : stomate_clear
3568!!
3569!>\BRIEF        Deallocate memory of the stomate variables.
3570!!
3571!! DESCRIPTION  : None
3572!!
3573!! RECENT CHANGE(S) : None
3574!!
3575!! MAIN OUTPUT VARIABLE(S): None
3576!!
3577!! REFERENCES   : None
3578!!
3579!! FLOWCHART    : None
3580!! \n
3581!_ ================================================================================================================================
3582 
3583  SUBROUTINE stomate_clear
3584
3585  !! 1. Deallocate all dynamics variables
3586
3587    IF (ALLOCATED(veget_cov_max)) DEALLOCATE(veget_cov_max)
3588    IF (ALLOCATED(ind)) DEALLOCATE(ind)
3589    IF (ALLOCATED(adapted)) DEALLOCATE(adapted)
3590    IF (ALLOCATED(regenerate)) DEALLOCATE(regenerate)
3591    IF (ALLOCATED(humrel_daily)) DEALLOCATE(humrel_daily)
3592    IF (ALLOCATED(gdd_init_date)) DEALLOCATE(gdd_init_date)
3593    IF (ALLOCATED(litterhum_daily)) DEALLOCATE(litterhum_daily)
3594    IF (ALLOCATED(t2m_daily))  DEALLOCATE(t2m_daily)
3595    IF (ALLOCATED(t2m_min_daily))  DEALLOCATE(t2m_min_daily)
3596    IF (ALLOCATED(tsurf_daily))  DEALLOCATE(tsurf_daily)
3597    IF (ALLOCATED(tsoil_daily)) DEALLOCATE(tsoil_daily)
3598    IF (ALLOCATED(soilhum_daily)) DEALLOCATE(soilhum_daily)
3599    IF (ALLOCATED(precip_daily)) DEALLOCATE(precip_daily)
3600    IF (ALLOCATED(gpp_daily)) DEALLOCATE(gpp_daily)
3601    IF (ALLOCATED(npp_daily)) DEALLOCATE(npp_daily)
3602    IF (ALLOCATED(turnover_daily)) DEALLOCATE(turnover_daily)
3603    IF (ALLOCATED(turnover_littercalc)) DEALLOCATE(turnover_littercalc)
3604    IF (ALLOCATED(humrel_month)) DEALLOCATE(humrel_month)
3605    IF (ALLOCATED(humrel_week)) DEALLOCATE(humrel_week)
3606    IF (ALLOCATED(t2m_longterm)) DEALLOCATE(t2m_longterm)
3607    IF (ALLOCATED(t2m_month)) DEALLOCATE(t2m_month)
3608    IF (ALLOCATED(Tseason)) DEALLOCATE(Tseason)
3609    IF (ALLOCATED(Tseason_length)) DEALLOCATE(Tseason_length)
3610    IF (ALLOCATED(Tseason_tmp)) DEALLOCATE(Tseason_tmp)
3611    IF (ALLOCATED(Tmin_spring_time)) DEALLOCATE(Tmin_spring_time)
3612    IF (ALLOCATED(onset_date)) DEALLOCATE(onset_date)
3613    IF (ALLOCATED(begin_leaves)) DEALLOCATE(begin_leaves)
3614    IF (ALLOCATED(t2m_week)) DEALLOCATE(t2m_week)
3615    IF (ALLOCATED(tsoil_month)) DEALLOCATE(tsoil_month)
3616    IF (ALLOCATED(soilhum_month)) DEALLOCATE(soilhum_month)
3617    IF (ALLOCATED(fireindex)) DEALLOCATE(fireindex)
3618    IF (ALLOCATED(firelitter)) DEALLOCATE(firelitter)
3619    IF (ALLOCATED(maxhumrel_lastyear)) DEALLOCATE(maxhumrel_lastyear)
3620    IF (ALLOCATED(maxhumrel_thisyear)) DEALLOCATE(maxhumrel_thisyear)
3621    IF (ALLOCATED(minhumrel_lastyear)) DEALLOCATE(minhumrel_lastyear)
3622    IF (ALLOCATED(minhumrel_thisyear)) DEALLOCATE(minhumrel_thisyear)
3623    IF (ALLOCATED(maxgppweek_lastyear)) DEALLOCATE(maxgppweek_lastyear)
3624    IF (ALLOCATED(maxgppweek_thisyear)) DEALLOCATE(maxgppweek_thisyear)
3625    IF (ALLOCATED(gdd0_lastyear)) DEALLOCATE(gdd0_lastyear)
3626    IF (ALLOCATED(gdd0_thisyear)) DEALLOCATE(gdd0_thisyear)
3627    IF (ALLOCATED(precip_lastyear)) DEALLOCATE(precip_lastyear)
3628    IF (ALLOCATED(precip_thisyear)) DEALLOCATE(precip_thisyear)
3629    IF (ALLOCATED(gdd_m5_dormance)) DEALLOCATE(gdd_m5_dormance)
3630    IF (ALLOCATED(gdd_from_growthinit)) DEALLOCATE(gdd_from_growthinit)
3631    IF (ALLOCATED(gdd_midwinter)) DEALLOCATE(gdd_midwinter)
3632    IF (ALLOCATED(ncd_dormance)) DEALLOCATE(ncd_dormance)
3633    IF (ALLOCATED(ngd_minus5))  DEALLOCATE(ngd_minus5)
3634    IF (ALLOCATED(PFTpresent)) DEALLOCATE(PFTpresent)
3635    IF (ALLOCATED(npp_longterm)) DEALLOCATE(npp_longterm)
3636    IF (ALLOCATED(lm_lastyearmax)) DEALLOCATE(lm_lastyearmax)
3637    IF (ALLOCATED(lm_thisyearmax)) DEALLOCATE(lm_thisyearmax)
3638    IF (ALLOCATED(maxfpc_lastyear)) DEALLOCATE(maxfpc_lastyear)
3639    IF (ALLOCATED(maxfpc_thisyear)) DEALLOCATE(maxfpc_thisyear)
3640    IF (ALLOCATED(turnover_longterm)) DEALLOCATE(turnover_longterm)
3641    IF (ALLOCATED(gpp_week)) DEALLOCATE(gpp_week)
3642    IF (ALLOCATED(biomass)) DEALLOCATE(biomass)
3643    IF (ALLOCATED(senescence)) DEALLOCATE(senescence)
3644    IF (ALLOCATED(when_growthinit)) DEALLOCATE(when_growthinit)
3645    IF (ALLOCATED(age))  DEALLOCATE(age)
3646    IF (ALLOCATED(resp_hetero_d)) DEALLOCATE(resp_hetero_d)
3647    IF (ALLOCATED(resp_hetero_radia)) DEALLOCATE(resp_hetero_radia)
3648    IF (ALLOCATED(resp_maint_d)) DEALLOCATE(resp_maint_d)
3649    IF (ALLOCATED(resp_growth_d)) DEALLOCATE(resp_growth_d)
3650    IF (ALLOCATED(co2_fire)) DEALLOCATE(co2_fire)
3651    IF (ALLOCATED(co2_to_bm_dgvm)) DEALLOCATE(co2_to_bm_dgvm)
3652    IF (ALLOCATED(veget_lastlight)) DEALLOCATE(veget_lastlight)
3653    IF (ALLOCATED(everywhere)) DEALLOCATE(everywhere)
3654    IF (ALLOCATED(need_adjacent)) DEALLOCATE(need_adjacent)
3655    IF (ALLOCATED(leaf_age)) DEALLOCATE(leaf_age)
3656    IF (ALLOCATED(leaf_frac)) DEALLOCATE(leaf_frac)
3657    IF (ALLOCATED(RIP_time)) DEALLOCATE(RIP_time)
3658    IF (ALLOCATED(time_hum_min)) DEALLOCATE(time_hum_min)
3659    IF (ALLOCATED(hum_min_dormance)) DEALLOCATE(hum_min_dormance)
3660    IF (ALLOCATED(litterpart)) DEALLOCATE(litterpart)
3661    IF (ALLOCATED(litter)) DEALLOCATE(litter)
3662    IF (ALLOCATED(dead_leaves)) DEALLOCATE(dead_leaves)
3663    IF (ALLOCATED(carbon)) DEALLOCATE(carbon)
3664    IF (ALLOCATED(lignin_struc)) DEALLOCATE(lignin_struc)
3665    IF (ALLOCATED(turnover_time)) DEALLOCATE(turnover_time)
3666    IF (ALLOCATED(co2_flux_daily)) DEALLOCATE(co2_flux_daily)
3667    IF (ALLOCATED(co2_flux_monthly)) DEALLOCATE(co2_flux_monthly)
3668    IF (ALLOCATED(harvest_above_monthly)) DEALLOCATE (harvest_above_monthly)
3669    IF (ALLOCATED(cflux_prod_monthly)) DEALLOCATE (cflux_prod_monthly)
3670    IF (ALLOCATED(bm_to_litter)) DEALLOCATE(bm_to_litter)
3671    IF (ALLOCATED(bm_to_littercalc)) DEALLOCATE(bm_to_littercalc)
3672    IF (ALLOCATED(herbivores)) DEALLOCATE(herbivores)
3673    IF (ALLOCATED(resp_maint_part_radia)) DEALLOCATE(resp_maint_part_radia)
3674    IF (ALLOCATED(resp_maint_radia)) DEALLOCATE(resp_maint_radia)
3675    IF (ALLOCATED(resp_maint_part)) DEALLOCATE(resp_maint_part)
3676    IF (ALLOCATED(hori_index)) DEALLOCATE(hori_index)
3677    IF (ALLOCATED(horipft_index)) DEALLOCATE(horipft_index)
3678    IF (ALLOCATED(clay_fm)) DEALLOCATE(clay_fm)
3679    IF (ALLOCATED(humrel_daily_fm)) DEALLOCATE(humrel_daily_fm)
3680    IF (ALLOCATED(litterhum_daily_fm))  DEALLOCATE(litterhum_daily_fm)
3681    IF (ALLOCATED(t2m_daily_fm))  DEALLOCATE(t2m_daily_fm)
3682    IF (ALLOCATED(t2m_min_daily_fm))  DEALLOCATE(t2m_min_daily_fm)
3683    IF (ALLOCATED(tsurf_daily_fm)) DEALLOCATE(tsurf_daily_fm)
3684    IF (ALLOCATED(tsoil_daily_fm)) DEALLOCATE(tsoil_daily_fm)
3685    IF (ALLOCATED(soilhum_daily_fm))  DEALLOCATE(soilhum_daily_fm)
3686    IF (ALLOCATED(precip_fm)) DEALLOCATE(precip_fm)
3687    IF (ALLOCATED(gpp_daily_fm))  DEALLOCATE(gpp_daily_fm)
3688    IF (ALLOCATED(veget_fm)) DEALLOCATE(veget_fm)
3689    IF (ALLOCATED(veget_max_fm)) DEALLOCATE(veget_max_fm)
3690    IF (ALLOCATED(lai_fm))  DEALLOCATE(lai_fm)
3691    !
3692    IF (ALLOCATED(ok_equilibrium)) DEALLOCATE(ok_equilibrium)
3693    IF (ALLOCATED(carbon_eq)) DEALLOCATE(carbon_eq)
3694    IF (ALLOCATED(matrixA)) DEALLOCATE(matrixA)
3695    IF (ALLOCATED(vectorB)) DEALLOCATE(vectorB)
3696    IF (ALLOCATED(MatrixV)) DEALLOCATE(MatrixV)
3697    IF (ALLOCATED(VectorU)) DEALLOCATE(VectorU)
3698    IF (ALLOCATED(MatrixW)) DEALLOCATE(MatrixW)
3699    IF (ALLOCATED(previous_stock)) DEALLOCATE(previous_stock)
3700    IF (ALLOCATED(current_stock)) DEALLOCATE(current_stock) 
3701    IF (ALLOCATED(nbp_accu)) DEALLOCATE(nbp_accu)
3702    IF (ALLOCATED(nbp_flux)) DEALLOCATE(nbp_flux)
3703
3704    IF (ALLOCATED(clay_fm_g)) DEALLOCATE(clay_fm_g)
3705    IF (ALLOCATED(humrel_daily_fm_g)) DEALLOCATE(humrel_daily_fm_g)
3706    IF (ALLOCATED(litterhum_daily_fm_g))  DEALLOCATE(litterhum_daily_fm_g)
3707    IF (ALLOCATED(t2m_daily_fm_g))  DEALLOCATE(t2m_daily_fm_g)
3708    IF (ALLOCATED(t2m_min_daily_fm_g))  DEALLOCATE(t2m_min_daily_fm_g)
3709    IF (ALLOCATED(tsurf_daily_fm_g)) DEALLOCATE(tsurf_daily_fm_g)
3710    IF (ALLOCATED(tsoil_daily_fm_g)) DEALLOCATE(tsoil_daily_fm_g)
3711    IF (ALLOCATED(soilhum_daily_fm_g))  DEALLOCATE(soilhum_daily_fm_g)
3712    IF (ALLOCATED(precip_fm_g)) DEALLOCATE(precip_fm_g)
3713    IF (ALLOCATED(gpp_daily_fm_g))  DEALLOCATE(gpp_daily_fm_g)
3714    IF (ALLOCATED(veget_fm_g)) DEALLOCATE(veget_fm_g)
3715    IF (ALLOCATED(veget_max_fm_g)) DEALLOCATE(veget_max_fm_g)
3716    IF (ALLOCATED(lai_fm_g))  DEALLOCATE(lai_fm_g)
3717   
3718    IF (ALLOCATED(isf)) DEALLOCATE(isf)
3719    IF (ALLOCATED(nf_written)) DEALLOCATE(nf_written)
3720    IF (ALLOCATED(nf_cumul)) DEALLOCATE(nf_cumul)
3721    IF (ALLOCATED(nforce)) DEALLOCATE(nforce)
3722    IF (ALLOCATED(control_moist)) DEALLOCATE(control_moist)
3723    IF (ALLOCATED(control_temp)) DEALLOCATE(control_temp)
3724    IF (ALLOCATED(soilcarbon_input)) DEALLOCATE(soilcarbon_input)
3725    IF ( ALLOCATED (horip10_index)) DEALLOCATE (horip10_index)
3726    IF ( ALLOCATED (horip100_index)) DEALLOCATE (horip100_index)
3727    IF ( ALLOCATED (horip11_index)) DEALLOCATE (horip11_index)
3728    IF ( ALLOCATED (horip101_index)) DEALLOCATE (horip101_index)
3729    IF ( ALLOCATED (prod10)) DEALLOCATE (prod10)
3730    IF ( ALLOCATED (prod100)) DEALLOCATE (prod100)
3731    IF ( ALLOCATED (flux10)) DEALLOCATE (flux10)
3732    IF ( ALLOCATED (flux100)) DEALLOCATE (flux100)
3733    IF ( ALLOCATED (convflux)) DEALLOCATE (convflux)
3734    IF ( ALLOCATED (cflux_prod10)) DEALLOCATE (cflux_prod10)
3735    IF ( ALLOCATED (cflux_prod100)) DEALLOCATE (cflux_prod100)
3736    IF ( ALLOCATED (harvest_above)) DEALLOCATE (harvest_above)
3737    IF ( ALLOCATED (soilcarbon_input_daily)) DEALLOCATE (soilcarbon_input_daily)
3738    IF ( ALLOCATED (control_temp_daily)) DEALLOCATE (control_temp_daily)
3739    IF ( ALLOCATED (control_moist_daily)) DEALLOCATE (control_moist_daily)
3740
3741    IF ( ALLOCATED (fpc_max)) DEALLOCATE (fpc_max)
3742!gmjc
3743    IF (ALLOCATED(litter_avail)) DEALLOCATE(litter_avail)
3744    IF (ALLOCATED(litter_not_avail)) DEALLOCATE(litter_not_avail)
3745    IF (ALLOCATED(litter_avail_frac)) DEALLOCATE(litter_avail_frac)
3746!    IF (ALLOCATED(resp_hetero_litter_d)) DEALLOCATE(resp_hetero_litter_d)
3747!    IF (ALLOCATED(resp_hetero_soil_d)) DEALLOCATE(resp_hetero_soil_d)
3748!    IF (ALLOCATED(resp_hetero_soil_part)) DEALLOCATE(resp_hetero_soil_part)
3749    IF (ALLOCATED(sla_calc)) DEALLOCATE(sla_calc)
3750    IF (ALLOCATED(wshtotsum)) DEALLOCATE(wshtotsum)
3751    IF (ALLOCATED(sr_ugb)) DEALLOCATE(sr_ugb)
3752    IF (ALLOCATED(compt_ugb)) DEALLOCATE(compt_ugb)
3753    IF (ALLOCATED(nb_ani)) DEALLOCATE(nb_ani)
3754    IF (ALLOCATED(grazed_frac)) DEALLOCATE(grazed_frac)
3755    IF (ALLOCATED(import_yield)) DEALLOCATE(import_yield)
3756    IF (ALLOCATED(t2m_14)) DEALLOCATE(t2m_14)
3757    IF (ALLOCATED(sla_age1)) DEALLOCATE(sla_age1)
3758    IF (ALLOCATED(N_limfert)) DEALLOCATE(N_limfert)
3759    IF (ALLOCATED(when_growthinit_cut)) DEALLOCATE(when_growthinit_cut)
3760    IF (ALLOCATED(nb_grazingdays)) DEALLOCATE(nb_grazingdays)
3761    IF (ALLOCATED(snowfall_daily)) DEALLOCATE(snowfall_daily)
3762    IF (ALLOCATED(snowmass_daily)) DEALLOCATE(snowmass_daily)
3763!JCADD top 5 layer grassland soil moisture for grazing
3764    IF (ALLOCATED(tmc_topgrass_daily)) DEALLOCATE(tmc_topgrass_daily)
3765    IF (ALLOCATED(after_snow)) DEALLOCATE(after_snow)
3766    IF (ALLOCATED(after_wet)) DEALLOCATE(after_wet)
3767    IF (ALLOCATED(wet1day)) DEALLOCATE(wet1day)
3768    IF (ALLOCATED(wet2day)) DEALLOCATE(wet2day)
3769!end gmjc
3770 !! 2. reset l_first
3771
3772    l_first_stomate=.TRUE.
3773
3774 !! 3. call to clear functions
3775
3776    CALL season_clear
3777    CALL stomatelpj_clear
3778    CALL littercalc_clear
3779    CALL vmax_clear
3780 
3781  END SUBROUTINE stomate_clear
3782
3783
3784!! ================================================================================================================================
3785!! SUBROUTINE   : stomate_var_init
3786!!
3787!>\BRIEF        Initialize variables of stomate with a none-zero initial value.
3788!! Subroutine is called only if ::ok_stomate = .TRUE. STOMATE diagnoses some
3789!! variables for SECHIBA : assim_param, deadleaf_cover, etc. These variables can
3790!! be recalculated from STOMATE's prognostic variables. Note that height is
3791!! saved in SECHIBA.
3792!!
3793!! DESCRIPTION  : None
3794!!
3795!! RECENT CHANGE(S) : None
3796!!
3797!! MAIN OUTPUT VARIABLE(S): leaf age (::leaf_age) and fraction of leaves in leaf
3798!! age class (::leaf_frac). The maximum water on vegetation available for
3799!! interception, fraction of soil covered by dead leaves
3800!! (::deadleaf_cover) and assimilation parameters (:: assim_param).
3801!!
3802!! REFERENCE(S) : None
3803!!
3804!! FLOWCHART    : None
3805!! \n
3806!_ ================================================================================================================================
3807 
3808  SUBROUTINE stomate_var_init &
3809       &  (kjpindex, veget_cov_max, leaf_age, leaf_frac, &
3810       &   dead_leaves, &
3811       &   veget, lai, deadleaf_cover, assim_param, &!)
3812!gmjc
3813       &    N_limfert)
3814!end gmjc
3815
3816  !! 0. Variable and parameter declaration
3817
3818    !! 0.1 Input variables
3819
3820    INTEGER(i_std),INTENT(in)                             :: kjpindex        !! Domain size - terrestrial pixels only
3821    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)        :: veget           !! Fraction of pixel covered by PFT. Fraction
3822                                                                             !! accounts for none-biological land covers
3823                                                                             !! (unitless)
3824    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)        :: veget_cov_max   !! Fractional coverage: maximum share of the pixel
3825                                                                             !! covered by a PFT (unitless)
3826    REAL(r_std),DIMENSION(kjpindex,nvm,nlitt),INTENT(in)  :: dead_leaves     !! Metabolic and structural fraction of dead leaves
3827                                                                             !! per ground area
3828                                                                             !! @tex $(gC m^{-2})$ @endtex
3829    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)        :: lai             !! Leaf area index
3830                                                                             !! @tex $(m^2 m{-2})$ @endtex
3831    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages),INTENT(in) :: leaf_age     !! Age of different leaf classes per PFT (days)
3832    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages),INTENT(in) :: leaf_frac    !! Fraction of leaves in leaf age class per PFT
3833                                                                             !! (unitless; 1)     
3834
3835    !! 0.2 Modified variables
3836    REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(inout) :: assim_param   !! min+max+opt temperatures (K) & vmax for
3837                                                                             !! photosynthesis 
3838   
3839    !! 0.3 Output variables
3840
3841    REAL(r_std),DIMENSION(kjpindex), INTENT (out)         :: deadleaf_cover  !! Fraction of soil covered by dead leaves
3842                                                                             !! (unitless)
3843
3844
3845!gmjc
3846    ! N fertilization factor on Vcmax and SLA
3847    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)       :: N_limfert
3848!end gmjc
3849    ! 0.4 Local variables
3850   
3851    REAL(r_std),PARAMETER                                 :: dt_0 = zero     !! Dummy time step, must be zero
3852    REAL(r_std),DIMENSION(kjpindex,nvm)                   :: vcmax           !! Dummy vcmax
3853                                                                             !! @tex $(\mu mol m^{-2} s^{-1})$ @endtex
3854    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages)         :: leaf_age_tmp    !! Temporary variable
3855    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages)         :: leaf_frac_tmp   !! Temporary variable
3856                                                                             !! (unitless; 1)     
3857    INTEGER(i_std)                                        :: j               !! Index (untiless)
3858   
3859!_ ================================================================================================================================   
3860
3861
3862    ! Calculate assim_param if it was not found in the restart file
3863    IF (ALL(assim_param(:,:,:)==val_exp)) THEN
3864       ! Use temporary leaf_age_tmp and leaf_frac_tmp to preserve the input variables from being modified by the subroutine vmax.
3865       leaf_age_tmp(:,:,:)=leaf_age(:,:,:)
3866       leaf_frac_tmp(:,:,:)=leaf_frac(:,:,:)
3867
3868       !! 1.1 Calculate a temporary vcmax (stomate_vmax.f90)
3869       CALL vmax (kjpindex, dt_0, leaf_age_tmp, leaf_frac_tmp, vcmax, N_limfert)
3870
3871       !! 1.2 transform into nvm vegetation types
3872       assim_param(:,:,ivcmax) = zero
3873       DO j = 2, nvm
3874          assim_param(:,j,ivcmax)=vcmax(:,j)
3875       ENDDO
3876    END IF
3877   
3878    !! 2. Dead leaf cover (stomate_litter.f90)
3879    CALL deadleaf (kjpindex, veget_cov_max, dead_leaves, deadleaf_cover, sla_calc)     
3880   
3881  END SUBROUTINE stomate_var_init
3882
3883
3884!! ================================================================================================================================
3885!! SUBROUTINE   : stomate_accu
3886!!
3887!>\BRIEF        Accumulate a variable for the time period specified by
3888!! ::dt_sechiba or calculate the mean value over ::dt_sechiba.
3889!!
3890!! DESCRIPTION : None
3891!!
3892!! RECENT CHANGE(S) : None
3893!!
3894!! MAIN OUTPUT VARIABLE(S): accumulated or mean variable ::field_out::
3895!!
3896!! REFERENCE(S) : None
3897!!
3898!! FLOWCHART    : None
3899!! \n
3900!_ ================================================================================================================================
3901 
3902  SUBROUTINE stomate_accu (npts, n_dim2, ldmean, field_in, field_out)
3903   
3904  !! 0. Variable and parameter declaration
3905
3906    !! 0.1 Input variables
3907    INTEGER(i_std),INTENT(in)                        :: npts      !! Domain size (unitless)
3908    INTEGER(i_std),INTENT(in)                        :: n_dim2    !! 2nd dimension (1 or nvm)
3909    LOGICAL,INTENT(in)                               :: ldmean    !! Flag to calculate the mean over
3910    REAL(r_std),DIMENSION(npts,n_dim2),INTENT(in)    :: field_in  !! Field that needs to be accumulated
3911   
3912    !! 0.2 Modified variables
3913    REAL(r_std),DIMENSION(npts,n_dim2),INTENT(inout) :: field_out !! Accumulated or mean field
3914
3915!_ ================================================================================================================================
3916
3917  !! 1. Accumulate field
3918
3919    field_out(:,:) = field_out(:,:)+field_in(:,:)*dt_sechiba
3920   
3921  !! 2. Mean fields
3922
3923    IF (ldmean) THEN
3924       field_out(:,:) = field_out(:,:)/dt_stomate
3925    ENDIF
3926
3927  END SUBROUTINE stomate_accu
3928
3929
3930!! ================================================================================================================================
3931!! SUBROUTINE   : init_forcing
3932!!
3933!>\BRIEF        Allocate memory for the variables containing the forcing data.
3934!! The maximum size of the allocated memory is specified in run definition file
3935!! (::max_totsize) and needs to be a compromise between charging the memory and
3936!! accessing disks to get the forcing data.
3937!!
3938!! DESCRIPTION : None
3939!!
3940!! RECENT CHANGE(S) : None
3941!!
3942!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
3943!! variables. However, the routine allocates memory for later use.
3944!!
3945!! REFERENCE(S) : None
3946!!
3947!! FLOWCHART    : None
3948!! \n
3949!_ ================================================================================================================================
3950 
3951  SUBROUTINE init_forcing (kjpindex,nsfm,nsft_loc)
3952   
3953  !! 0. Variable and parameter declaration
3954
3955    !! 0.1 Input variables
3956    INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
3957    INTEGER(i_std),INTENT(in) :: nsfm     !! Number of time steps that can be stored in memory (unitless)
3958    INTEGER(i_std),INTENT(in) :: nsft_loc !! Number of time steps in a year (unitless)
3959
3960   !! 0.2 Output variables
3961
3962   !! 0.3 Modified variables
3963
3964   !! 0.4 Local variables
3965
3966    LOGICAL                   :: l_error  !! Check errors in netcdf call
3967    INTEGER(i_std)            :: ier      !! Check errors in netcdf call
3968!_ ================================================================================================================================
3969   
3970  !! 1. Allocate memory
3971
3972    ! Note ::nvm is number of PFTs and ::nbdl is number of soil layers
3973    l_error = .FALSE.
3974    ALLOCATE(clay_fm(kjpindex,nsfm),stat=ier)
3975    l_error = l_error .OR. (ier /= 0)
3976    IF (l_error) THEN
3977       WRITE(numout,*) 'Problem with memory allocation: forcing variables clay_fm ',kjpindex,nsfm
3978       STOP 'init_forcing'
3979    ENDIF
3980    ALLOCATE(humrel_daily_fm(kjpindex,nvm,nsfm),stat=ier)
3981    l_error = l_error .OR. (ier /= 0)
3982    IF (l_error) THEN
3983       WRITE(numout,*) 'Problem with memory allocation: forcing variables humrel_daily_fm ',kjpindex,nvm,nsfm
3984       STOP 'init_forcing'
3985    ENDIF
3986    ALLOCATE(litterhum_daily_fm(kjpindex,nsfm),stat=ier)
3987    l_error = l_error .OR. (ier /= 0)
3988    IF (l_error) THEN
3989       WRITE(numout,*) 'Problem with memory allocation: forcing variables litterhum_daily_fm ',kjpindex,nsfm
3990       STOP 'init_forcing'
3991    ENDIF
3992    ALLOCATE(t2m_daily_fm(kjpindex,nsfm),stat=ier)
3993    l_error = l_error .OR. (ier /= 0)
3994    IF (l_error) THEN
3995       WRITE(numout,*) 'Problem with memory allocation: forcing variables t2m_daily_fm ',kjpindex,nsfm
3996       STOP 'init_forcing'
3997    ENDIF
3998    ALLOCATE(t2m_min_daily_fm(kjpindex,nsfm),stat=ier)
3999    l_error = l_error .OR. (ier /= 0)
4000    IF (l_error) THEN
4001       WRITE(numout,*) 'Problem with memory allocation: forcing variables t2m_min_daily_fm ',kjpindex,nsfm
4002       STOP 'init_forcing'
4003    ENDIF
4004    ALLOCATE(tsurf_daily_fm(kjpindex,nsfm),stat=ier)
4005    l_error = l_error .OR. (ier /= 0)
4006    IF (l_error) THEN
4007       WRITE(numout,*) 'Problem with memory allocation: forcing variables tsurf_daily_fm ',kjpindex,nsfm
4008       STOP 'init_forcing'
4009    ENDIF
4010    ALLOCATE(tsoil_daily_fm(kjpindex,nbdl,nsfm),stat=ier)
4011    l_error = l_error .OR. (ier /= 0)
4012    IF (l_error) THEN
4013       WRITE(numout,*) 'Problem with memory allocation: forcing variables tsoil_daily_fm ',kjpindex,nbdl,nsfm
4014       STOP 'init_forcing'
4015    ENDIF
4016    ALLOCATE(soilhum_daily_fm(kjpindex,nbdl,nsfm),stat=ier)
4017    l_error = l_error .OR. (ier /= 0)
4018    IF (l_error) THEN
4019       WRITE(numout,*) 'Problem with memory allocation: forcing variables soilhum_daily_fm ',kjpindex,nbdl,nsfm
4020       STOP 'init_forcing'
4021    ENDIF
4022    ALLOCATE(precip_fm(kjpindex,nsfm),stat=ier)
4023    l_error = l_error .OR. (ier /= 0)
4024    IF (l_error) THEN
4025       WRITE(numout,*) 'Problem with memory allocation: forcing variables precip_fm ',kjpindex,nsfm
4026       STOP 'init_forcing'
4027    ENDIF
4028    ALLOCATE(gpp_daily_fm(kjpindex,nvm,nsfm),stat=ier)
4029    l_error = l_error .OR. (ier /= 0)
4030    IF (l_error) THEN
4031       WRITE(numout,*) 'Problem with memory allocation: forcing variables gpp_daily_fm ',kjpindex,nvm,nsfm
4032       STOP 'init_forcing'
4033    ENDIF
4034    ALLOCATE(veget_fm(kjpindex,nvm,nsfm),stat=ier)
4035    l_error = l_error .OR. (ier /= 0)
4036    IF (l_error) THEN
4037       WRITE(numout,*) 'Problem with memory allocation: forcing variables veget_fm ',kjpindex,nvm,nsfm
4038       STOP 'init_forcing'
4039    ENDIF
4040    ALLOCATE(veget_max_fm(kjpindex,nvm,nsfm),stat=ier)
4041    l_error = l_error .OR. (ier /= 0)
4042    IF (l_error) THEN
4043       WRITE(numout,*) 'Problem with memory allocation: forcing variables veget_max_fm ',kjpindex,nvm,nsfm
4044       STOP 'init_forcing'
4045    ENDIF
4046    ALLOCATE(lai_fm(kjpindex,nvm,nsfm),stat=ier)
4047    l_error = l_error .OR. (ier /= 0)
4048    IF (l_error) THEN
4049       WRITE(numout,*) 'Problem with memory allocation: forcing variables lai_fm ',kjpindex,nvm,nsfm
4050       STOP 'init_forcing'
4051    ENDIF
4052    ALLOCATE(isf(nsfm),stat=ier)
4053    l_error = l_error .OR. (ier /= 0)
4054    IF (l_error) THEN
4055       WRITE(numout,*) 'Problem with memory allocation: forcing variables isf ',nsfm
4056       STOP 'init_forcing'
4057    ENDIF
4058    ALLOCATE(nf_written(nsft_loc),stat=ier)
4059    l_error = l_error .OR. (ier /= 0)
4060    IF (l_error) THEN
4061       WRITE(numout,*) 'Problem with memory allocation: forcing variables nf_written ',nsft_loc
4062       STOP 'init_forcing'
4063    ENDIF
4064    ALLOCATE(nf_cumul(nsft_loc),stat=ier)
4065    l_error = l_error .OR. (ier /= 0)
4066    IF (l_error) THEN
4067       WRITE(numout,*) 'Problem with memory allocation: forcing variables nf_cumul ',nsft_loc
4068       STOP 'init_forcing'
4069    ENDIF
4070   
4071  !! 2. Allocate memory for the root processor only (parallel computing)
4072
4073    ! Where, ::nbp_glo is the number of global continental points
4074    IF (is_root_prc) THEN
4075       ALLOCATE(clay_fm_g(nbp_glo,nsfm),stat=ier)
4076       l_error = l_error .OR. (ier /= 0)
4077       IF (l_error) THEN
4078          WRITE(numout,*) 'Problem with memory allocation: forcing variables clay_fm_g ',nbp_glo,nsfm
4079          STOP 'init_forcing'
4080       ENDIF
4081       ALLOCATE(humrel_daily_fm_g(nbp_glo,nvm,nsfm),stat=ier)
4082       l_error = l_error .OR. (ier /= 0)
4083       IF (l_error) THEN
4084          WRITE(numout,*) 'Problem with memory allocation: forcing variables humrel_daily_fm_g ',nbp_glo,nvm,nsfm
4085          STOP 'init_forcing'
4086       ENDIF
4087       ALLOCATE(litterhum_daily_fm_g(nbp_glo,nsfm),stat=ier)
4088       l_error = l_error .OR. (ier /= 0)
4089       IF (l_error) THEN
4090          WRITE(numout,*) 'Problem with memory allocation: forcing variables litterhum_daily_fm_g ',nbp_glo,nsfm
4091          STOP 'init_forcing'
4092       ENDIF
4093       ALLOCATE(t2m_daily_fm_g(nbp_glo,nsfm),stat=ier)
4094       l_error = l_error .OR. (ier /= 0)
4095       IF (l_error) THEN
4096          WRITE(numout,*) 'Problem with memory allocation: forcing variables t2m_daily_fm_g ',nbp_glo,nsfm
4097          STOP 'init_forcing'
4098       ENDIF
4099       ALLOCATE(t2m_min_daily_fm_g(nbp_glo,nsfm),stat=ier)
4100       l_error = l_error .OR. (ier /= 0)
4101       IF (l_error) THEN
4102          WRITE(numout,*) 'Problem with memory allocation: forcing variables t2m_min_daily_fm_g ',nbp_glo,nsfm
4103          STOP 'init_forcing'
4104       ENDIF
4105       ALLOCATE(tsurf_daily_fm_g(nbp_glo,nsfm),stat=ier)
4106       l_error = l_error .OR. (ier /= 0)
4107       IF (l_error) THEN
4108          WRITE(numout,*) 'Problem with memory allocation: forcing variables tsurf_daily_fm_g ',nbp_glo,nsfm
4109          STOP 'init_forcing'
4110       ENDIF
4111       ALLOCATE(tsoil_daily_fm_g(nbp_glo,nbdl,nsfm),stat=ier)
4112       l_error = l_error .OR. (ier /= 0)
4113       IF (l_error) THEN
4114          WRITE(numout,*) 'Problem with memory allocation: forcing variables tsoil_daily_fm_g ',nbp_glo,nbdl,nsfm
4115          STOP 'init_forcing'
4116       ENDIF
4117       ALLOCATE(soilhum_daily_fm_g(nbp_glo,nbdl,nsfm),stat=ier)
4118       l_error = l_error .OR. (ier /= 0)
4119       IF (l_error) THEN
4120          WRITE(numout,*) 'Problem with memory allocation: forcing variables soilhum_daily_fm_g ',nbp_glo,nbdl,nsfm
4121          STOP 'init_forcing'
4122       ENDIF
4123       ALLOCATE(precip_fm_g(nbp_glo,nsfm),stat=ier)
4124       l_error = l_error .OR. (ier /= 0)
4125       IF (l_error) THEN
4126          WRITE(numout,*) 'Problem with memory allocation: forcing variables precip_fm_g ',nbp_glo,nsfm
4127          STOP 'init_forcing'
4128       ENDIF
4129       ALLOCATE(gpp_daily_fm_g(nbp_glo,nvm,nsfm),stat=ier)
4130       l_error = l_error .OR. (ier /= 0)
4131       IF (l_error) THEN
4132          WRITE(numout,*) 'Problem with memory allocation: forcing variables gpp_daily_fm_g ',nbp_glo,nvm,nsfm
4133          STOP 'init_forcing'
4134       ENDIF
4135       ALLOCATE(veget_fm_g(nbp_glo,nvm,nsfm),stat=ier)
4136       l_error = l_error .OR. (ier /= 0)
4137       IF (l_error) THEN
4138          WRITE(numout,*) 'Problem with memory allocation: forcing variables veget_fm_g ',nbp_glo,nvm,nsfm
4139          STOP 'init_forcing'
4140       ENDIF
4141       ALLOCATE(veget_max_fm_g(nbp_glo,nvm,nsfm),stat=ier)
4142       l_error = l_error .OR. (ier /= 0)
4143       IF (l_error) THEN
4144          WRITE(numout,*) 'Problem with memory allocation: forcing variables veget_max_fm_g ',nbp_glo,nvm,nsfm
4145          STOP 'init_forcing'
4146       ENDIF
4147       ALLOCATE(lai_fm_g(nbp_glo,nvm,nsfm),stat=ier)
4148       l_error = l_error .OR. (ier /= 0)
4149       IF (l_error) THEN
4150          WRITE(numout,*) 'Problem with memory allocation: forcing variables lai_fm_g ',nbp_glo,nvm,nsfm
4151          STOP 'init_forcing'
4152       ENDIF
4153    ELSE
4154       ! Allocate memory for co-processors
4155       ALLOCATE(clay_fm_g(0,nsfm),stat=ier)
4156       ALLOCATE(humrel_daily_fm_g(0,nvm,nsfm),stat=ier)
4157       ALLOCATE(litterhum_daily_fm_g(0,nsfm),stat=ier)
4158       ALLOCATE(t2m_daily_fm_g(0,nsfm),stat=ier)
4159       ALLOCATE(t2m_min_daily_fm_g(0,nsfm),stat=ier)
4160       ALLOCATE(tsurf_daily_fm_g(0,nsfm),stat=ier)
4161       ALLOCATE(tsoil_daily_fm_g(0,nbdl,nsfm),stat=ier)
4162       ALLOCATE(soilhum_daily_fm_g(0,nbdl,nsfm),stat=ier)
4163       ALLOCATE(precip_fm_g(0,nsfm),stat=ier)
4164       ALLOCATE(gpp_daily_fm_g(0,nvm,nsfm),stat=ier)
4165       ALLOCATE(veget_fm_g(0,nvm,nsfm),stat=ier)
4166       ALLOCATE(veget_max_fm_g(0,nvm,nsfm),stat=ier)
4167       ALLOCATE(lai_fm_g(0,nvm,nsfm),stat=ier)
4168    ENDIF ! is_root_proc
4169   
4170    IF (l_error) THEN
4171       WRITE(numout,*) 'Problem with memory allocation: forcing variables'
4172       STOP 'init_forcing'
4173    ENDIF
4174
4175  !! 3. Initilaize variables
4176
4177    CALL forcing_zero
4178   
4179  END SUBROUTINE init_forcing
4180
4181
4182!! ================================================================================================================================
4183!! SUBROUTINE   : forcing_zero
4184!!
4185!>\BRIEF        Initialize variables containing the forcing data; variables are
4186!! set to zero.
4187!!
4188!! DESCRIPTION  : None
4189!!
4190!! RECENT CHANGE(S) : None
4191!!
4192!! MAIN OUTPUT VARIABLE(S): None
4193!!
4194!! REFERENCES   : None
4195!!
4196!! FLOWCHART    : None
4197!! \n
4198!_ ================================================================================================================================
4199 
4200  SUBROUTINE forcing_zero
4201   
4202    clay_fm(:,:) = zero
4203    humrel_daily_fm(:,:,:) = zero
4204    litterhum_daily_fm(:,:) = zero
4205    t2m_daily_fm(:,:) = zero
4206    t2m_min_daily_fm(:,:) = zero
4207    tsurf_daily_fm(:,:) = zero
4208    tsoil_daily_fm(:,:,:) = zero
4209    soilhum_daily_fm(:,:,:) = zero
4210    precip_fm(:,:) = zero
4211    gpp_daily_fm(:,:,:) = zero
4212    veget_fm(:,:,:) = zero
4213    veget_max_fm(:,:,:) = zero
4214    lai_fm(:,:,:) = zero
4215   
4216  END SUBROUTINE forcing_zero
4217
4218
4219!! ================================================================================================================================
4220!! SUBROUTINE   : forcing_write
4221!!
4222!>\BRIEF        Appends data values to a netCDF file containing the forcing
4223!! variables of the general processes in stomate.
4224!!
4225!! DESCRIPTION  : None
4226!!
4227!! RECENT CHANGE(S) : None
4228!!
4229!! MAIN OUTPUT VARIABLE(S): netCDF file
4230!!
4231!! REFERENCES   : None
4232!!
4233!! FLOWCHART    : None
4234!! \n
4235!_ ================================================================================================================================
4236 
4237  SUBROUTINE forcing_write(forcing_id,ibeg,iend)
4238   
4239  !! 0. Variable and parameter declaration
4240
4241    !! 0.1 Input variables
4242
4243    INTEGER(i_std),INTENT(in)      :: forcing_id  !! File identifer of forcing file, assigned when netcdf is created
4244    INTEGER(i_std),INTENT(in)      :: ibeg, iend  !! First and last time step to be written
4245
4246    !! 0.2 Output variables
4247
4248    !! 0.3 Modified variables
4249
4250    !! 0.4 Local variables
4251
4252    INTEGER(i_std)                 :: ii          !! Index of isf where isf is the number of time steps that can be
4253                                                  !! stored in memory
4254    INTEGER(i_std)                 :: iblocks     !! Index of block that is written
4255    INTEGER(i_std)                 :: nblocks     !! Number of blocks that needs to be written
4256    INTEGER(i_std)                 :: ier         !! Check errors in netcdf call
4257    INTEGER(i_std),DIMENSION(0:2)  :: ifirst      !! First block in memory - changes with iblocks
4258    INTEGER(i_std),DIMENSION(0:2)  :: ilast       !! Last block in memory - changes with iblocks
4259    INTEGER(i_std),PARAMETER       :: ndm = 10    !! Maximum number of dimensions
4260    INTEGER(i_std),DIMENSION(ndm)  :: start       !! First block to write
4261    INTEGER(i_std)                 :: ndim        !! Dimensions of forcing to be added to the netCDF
4262    INTEGER(i_std),DIMENSION(ndm)  :: count_force !! Number of elements in each dimension 
4263    INTEGER(i_std)                 :: vid         !! Variable identifer of netCDF
4264!_ ================================================================================================================================
4265   
4266  !! 1. Determine number of blocks of forcing variables that are stored in memory
4267
4268    nblocks = 0
4269    ifirst(:) = 1
4270    ilast(:) = 1
4271    DO ii = ibeg, iend
4272       IF (     (nblocks /= 0) &
4273            &      .AND.(isf(ii) == isf(ilast(nblocks))+1)) THEN
4274          ! Last block found
4275          ilast(nblocks) = ii
4276       ELSE
4277          ! First block found
4278          nblocks = nblocks+1
4279          IF (nblocks > 2)  STOP 'Problem in forcing_write'
4280          ifirst(nblocks) = ii
4281          ilast(nblocks) = ii
4282       ENDIF
4283    ENDDO
4284
4285  !! 2. Gather distributed variables (parallel computing)
4286
4287    CALL gather(clay_fm,clay_fm_g)
4288    CALL gather(humrel_daily_fm,humrel_daily_fm_g)
4289    CALL gather(litterhum_daily_fm,litterhum_daily_fm_g)
4290    CALL gather(t2m_daily_fm,t2m_daily_fm_g)
4291    CALL gather(t2m_min_daily_fm,t2m_min_daily_fm_g)
4292    CALL gather(tsurf_daily_fm,tsurf_daily_fm_g)
4293    CALL gather(tsoil_daily_fm,tsoil_daily_fm_g)
4294    CALL gather(soilhum_daily_fm,soilhum_daily_fm_g)
4295    CALL gather(precip_fm,precip_fm_g)
4296    CALL gather(gpp_daily_fm,gpp_daily_fm_g)
4297    CALL gather(veget_fm,veget_fm_g)
4298    CALL gather(veget_max_fm,veget_max_fm_g)
4299    CALL gather(lai_fm,lai_fm_g)
4300 
4301 !! 3. Append data to netCDF file
4302   
4303    IF (is_root_prc) THEN
4304       ! The netCDF file has been created earlier in this module, a file ID is available
4305       ! and variables and dimensions have already been defined
4306       DO iblocks = 1, nblocks
4307          IF (ifirst(iblocks) /= ilast(iblocks)) THEN
4308             ndim = 2
4309             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4310             count_force(1:ndim) = SHAPE(clay_fm_g)
4311             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4312             ier = NF90_INQ_VARID (forcing_id,'clay',vid)
4313             ier = NF90_PUT_VAR (forcing_id,vid, &
4314                  &              clay_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4315                  & start=start(1:ndim), count=count_force(1:ndim))
4316             ndim = 3;
4317             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4318             count_force(1:ndim) = SHAPE(humrel_daily_fm_g)
4319             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4320             ier = NF90_INQ_VARID (forcing_id,'humrel',vid)
4321             ier = NF90_PUT_VAR (forcing_id, vid, &
4322                  &            humrel_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4323                  &            start=start(1:ndim), count=count_force(1:ndim))
4324             ndim = 2;
4325             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4326             count_force(1:ndim) = SHAPE(litterhum_daily_fm_g)
4327             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4328             ier = NF90_INQ_VARID (forcing_id,'litterhum',vid)
4329             ier = NF90_PUT_VAR (forcing_id, vid, &
4330                  &            litterhum_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4331                  & start=start(1:ndim), count=count_force(1:ndim))
4332             ndim = 2;
4333             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4334             count_force(1:ndim) = SHAPE(t2m_daily_fm_g)
4335             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4336             ier = NF90_INQ_VARID (forcing_id,'t2m',vid)
4337             ier = NF90_PUT_VAR (forcing_id, vid, &
4338                  &            t2m_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4339                  & start=start(1:ndim), count=count_force(1:ndim))
4340             ndim = 2;
4341             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4342             count_force(1:ndim) = SHAPE(t2m_min_daily_fm_g)
4343             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4344             ier = NF90_INQ_VARID (forcing_id,'t2m_min',vid)
4345             ier = NF90_PUT_VAR (forcing_id, vid, &
4346                  &            t2m_min_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4347                  & start=start(1:ndim), count=count_force(1:ndim))
4348             ndim = 2;
4349             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4350             count_force(1:ndim) = SHAPE(tsurf_daily_fm_g)
4351             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4352             ier = NF90_INQ_VARID (forcing_id,'tsurf',vid)
4353             ier = NF90_PUT_VAR (forcing_id, vid, &
4354                  &            tsurf_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4355                  & start=start(1:ndim), count=count_force(1:ndim))
4356             ndim = 3;
4357             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4358             count_force(1:ndim) = SHAPE(tsoil_daily_fm_g)
4359             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4360             ier = NF90_INQ_VARID (forcing_id,'tsoil',vid)
4361             ier = NF90_PUT_VAR (forcing_id, vid, &
4362                  &            tsoil_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4363                  & start=start(1:ndim), count=count_force(1:ndim))
4364             ndim = 3;
4365             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4366             count_force(1:ndim) = SHAPE(soilhum_daily_fm_g)
4367             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4368             ier = NF90_INQ_VARID (forcing_id,'soilhum',vid)
4369             ier = NF90_PUT_VAR (forcing_id, vid, &
4370                  &            soilhum_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4371                  & start=start(1:ndim), count=count_force(1:ndim))
4372             ndim = 2;
4373             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4374             count_force(1:ndim) = SHAPE(precip_fm_g)
4375             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4376             ier = NF90_INQ_VARID (forcing_id,'precip',vid)
4377             ier = NF90_PUT_VAR (forcing_id, vid, &
4378                  &            precip_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4379                  & start=start(1:ndim), count=count_force(1:ndim))
4380             ndim = 3;
4381             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4382             count_force(1:ndim) = SHAPE(gpp_daily_fm_g)
4383             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4384             ier = NF90_INQ_VARID (forcing_id,'gpp',vid)
4385             ier = NF90_PUT_VAR (forcing_id, vid, &
4386                  &            gpp_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4387                  &            start=start(1:ndim), count=count_force(1:ndim))
4388             ndim = 3;
4389             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4390             count_force(1:ndim) = SHAPE(veget_fm_g)
4391             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4392             ier = NF90_INQ_VARID (forcing_id,'veget',vid)
4393             ier = NF90_PUT_VAR (forcing_id, vid, &
4394                  &            veget_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4395                  &            start=start(1:ndim), count=count_force(1:ndim))
4396             ndim = 3;
4397             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4398             count_force(1:ndim) = SHAPE(veget_max_fm_g)
4399             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4400             ier = NF90_INQ_VARID (forcing_id,'veget_max',vid)
4401             ier = NF90_PUT_VAR (forcing_id, vid, &
4402                  &            veget_max_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4403                  &            start=start(1:ndim), count=count_force(1:ndim))
4404             ndim = 3;
4405             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4406             count_force(1:ndim) = SHAPE(lai_fm_g)
4407             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4408             ier = NF90_INQ_VARID (forcing_id,'lai',vid)
4409             ier = NF90_PUT_VAR (forcing_id, vid, &
4410                  &            lai_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4411                  &            start=start(1:ndim), count=count_force(1:ndim))
4412          ENDIF
4413       ENDDO
4414    ENDIF
4415   
4416  !! 4. Adjust flag of forcing file
4417    nf_written(isf(:)) = .TRUE.
4418
4419  END SUBROUTINE forcing_write
4420
4421 
4422!! ================================================================================================================================
4423!! SUBROUTINE   : stomate_forcing_read
4424!!
4425!>\BRIEF        Read forcing file.
4426!!
4427!! DESCRIPTION  : None
4428!!
4429!! RECENT CHANGE(S) : None
4430!!
4431!! MAIN OUTPUT VARIABLE(S): None
4432!!
4433!! REFERENCES   : None
4434!!
4435!! FLOWCHART    : None
4436!! \n
4437!_ ================================================================================================================================
4438 
4439  SUBROUTINE stomate_forcing_read(forcing_id,nsfm)
4440   
4441  !! 0. Variable and parameter declaration
4442
4443    !! 0.1 Input variables
4444
4445    INTEGER(i_std),INTENT(in)  :: forcing_id           !! File identifer of forcing file, assigned when netcdf is created
4446    INTEGER(i_std),INTENT(in)  :: nsfm                 !! Number of time steps stored in memory       
4447   
4448    !! 0.2 Output variables
4449
4450    !! 0.3 Modified variables
4451
4452    !! 0.4 Local variables
4453
4454    INTEGER(i_std)                 :: ii                !! Index of isf where isf is the number of time steps that can be stored in
4455                                                        !! memory
4456    INTEGER(i_std)                 :: iblocks           !! Index of block that is written
4457    INTEGER(i_std)                 :: nblocks           !! Number of blocks that needs to be written
4458    INTEGER(i_std)                 :: ier               !! Check error of netcdf call
4459    INTEGER(i_std),DIMENSION(0:2)  :: ifirst            !! First block in memory - changes with iblocks
4460    INTEGER(i_std),DIMENSION(0:2)  :: ilast             !! Last block in memory - changes with iblocks
4461    INTEGER(i_std),PARAMETER       :: ndm = 10          !! Maximum number of dimensions
4462    INTEGER(i_std),DIMENSION(ndm)  :: start             !! First block to write
4463    INTEGER(i_std)                 :: ndim              !! Dimensions of forcing to be added to the netCDF
4464    INTEGER(i_std),DIMENSION(ndm)  :: count_force       !! Number of elements in each dimension
4465    INTEGER(i_std)                 :: vid               !! Variable identifer of netCDF
4466    LOGICAL, PARAMETER             :: check=.FALSE.     !! Flag for debugging
4467    LOGICAL                        :: a_er=.FALSE.      !! Error catching from netcdf file
4468!_ ================================================================================================================================
4469
4470    IF (check) WRITE(numout,*) "stomate_forcing_read "
4471   
4472  !! 1. Set to zero if the corresponding forcing state
4473
4474    ! has not yet been written into the file 
4475    DO ii = 1, nsfm
4476       IF (.NOT.nf_written(isf(ii))) THEN
4477          clay_fm(:,ii) = zero
4478          humrel_daily_fm(:,:,ii) = zero
4479          litterhum_daily_fm(:,ii) = zero
4480          t2m_daily_fm(:,ii) = zero
4481          t2m_min_daily_fm(:,ii) = zero
4482          tsurf_daily_fm(:,ii) = zero
4483          tsoil_daily_fm(:,:,ii) = zero
4484          soilhum_daily_fm(:,:,ii) = zero
4485          precip_fm(:,ii) = zero
4486          gpp_daily_fm(:,:,ii) = zero
4487          veget_fm(:,:,ii) = zero
4488          veget_max_fm(:,:,ii) = zero
4489          lai_fm(:,:,ii) = zero
4490       ENDIF
4491    ENDDO
4492   
4493  !! 2. determine blocks of forcing states that are contiguous in memory
4494
4495    nblocks = 0
4496    ifirst(:) = 1
4497    ilast(:) = 1
4498   
4499    DO ii = 1, nsfm
4500       IF (nf_written(isf(ii))) THEN
4501          IF (     (nblocks /= 0) &
4502               &        .AND.(isf(ii) == isf(ilast(nblocks))+1)) THEN
4503
4504             ! element is contiguous with last element found
4505             ilast(nblocks) = ii
4506          ELSE
4507
4508             ! found first element of new block
4509             nblocks = nblocks+1
4510             IF (nblocks > 2)  STOP 'Problem in stomate_forcing_read'
4511             
4512             ifirst(nblocks) = ii
4513             ilast(nblocks) = ii
4514          ENDIF
4515       ENDIF
4516    ENDDO
4517    IF (check) WRITE(numout,*) "stomate_forcing_read nblocks, ifirst, ilast",nblocks, ifirst, ilast
4518   
4519  !! 3. Read variable values
4520
4521    IF (is_root_prc) THEN
4522       DO iblocks = 1, nblocks
4523          IF (check) WRITE(numout,*) "stomate_forcing_read iblocks, ifirst(iblocks), ilast(iblocks)",iblocks, &
4524               ifirst(iblocks), ilast(iblocks)
4525          IF (ifirst(iblocks) /= ilast(iblocks)) THEN
4526             a_er=.FALSE.
4527             ndim = 2;
4528             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4529             count_force(1:ndim) = SHAPE(clay_fm_g)
4530             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4531             ier = NF90_INQ_VARID (forcing_id,'clay',vid)
4532             a_er = a_er.OR.(ier /= 0)
4533             ier = NF90_GET_VAR (forcing_id, vid, &
4534                  &            clay_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4535                  &            start=start(1:ndim), count=count_force(1:ndim))
4536             a_er = a_er.OR.(ier /= 0)
4537
4538             ndim = 3;
4539             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4540             count_force(1:ndim) = SHAPE(humrel_daily_fm_g)
4541             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4542             ier = NF90_INQ_VARID (forcing_id,'humrel',vid)
4543             a_er = a_er.OR.(ier /= 0)
4544             ier = NF90_GET_VAR (forcing_id, vid, &
4545                  &            humrel_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4546                  &            start=start(1:ndim), count=count_force(1:ndim))
4547             a_er = a_er.OR.(ier /= 0)
4548
4549             ndim = 2;
4550             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4551             count_force(1:ndim) = SHAPE(litterhum_daily_fm_g)
4552             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4553             ier = NF90_INQ_VARID (forcing_id,'litterhum',vid)
4554             a_er = a_er.OR.(ier /= 0)
4555             ier = NF90_GET_VAR (forcing_id, vid, &
4556                  &              litterhum_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4557                  &            start=start(1:ndim), count=count_force(1:ndim))
4558             a_er = a_er.OR.(ier /= 0)
4559
4560             ndim = 2;
4561             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4562             count_force(1:ndim) = SHAPE(t2m_daily_fm_g)
4563             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4564             ier = NF90_INQ_VARID (forcing_id,'t2m',vid)
4565             a_er = a_er.OR.(ier /= 0)
4566             ier = NF90_GET_VAR (forcing_id, vid, &
4567                  &              t2m_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4568                  &            start=start(1:ndim), count=count_force(1:ndim))
4569             a_er = a_er.OR.(ier /= 0)
4570
4571             ndim = 2;
4572             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4573             count_force(1:ndim) = SHAPE(t2m_min_daily_fm_g)
4574             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4575             ier = NF90_INQ_VARID (forcing_id,'t2m_min',vid)
4576             a_er = a_er.OR.(ier /= 0)
4577             ier = NF90_GET_VAR (forcing_id, vid, &
4578                  &              t2m_min_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4579                  &            start=start(1:ndim), count=count_force(1:ndim))
4580             a_er = a_er.OR.(ier /= 0)
4581
4582             ndim = 2;
4583             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4584             count_force(1:ndim) = SHAPE(tsurf_daily_fm_g)
4585             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4586             ier = NF90_INQ_VARID (forcing_id,'tsurf',vid)
4587             a_er = a_er.OR.(ier /= 0)
4588             ier = NF90_GET_VAR (forcing_id, vid, &
4589                  &              tsurf_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4590                  &            start=start(1:ndim), count=count_force(1:ndim))
4591             a_er = a_er.OR.(ier /= 0)
4592
4593             ndim = 3;
4594             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4595             count_force(1:ndim) = SHAPE(tsoil_daily_fm_g)
4596             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4597             ier = NF90_INQ_VARID (forcing_id,'tsoil',vid)
4598             a_er = a_er.OR.(ier /= 0)
4599             ier = NF90_GET_VAR (forcing_id, vid, &
4600                  &              tsoil_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4601                  &            start=start(1:ndim), count=count_force(1:ndim))
4602             a_er = a_er.OR.(ier /= 0)
4603
4604             ndim = 3;
4605             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4606             count_force(1:ndim) = SHAPE(soilhum_daily_fm_g)
4607             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4608             ier = NF90_INQ_VARID (forcing_id,'soilhum',vid)
4609             a_er = a_er.OR.(ier /= 0)
4610             ier = NF90_GET_VAR (forcing_id, vid, &
4611                  &              soilhum_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4612                  &            start=start(1:ndim), count=count_force(1:ndim))
4613             a_er = a_er.OR.(ier /= 0)
4614
4615             ndim = 2;
4616             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4617             count_force(1:ndim) = SHAPE(precip_fm_g)
4618             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4619             ier = NF90_INQ_VARID (forcing_id,'precip',vid)
4620             a_er = a_er.OR.(ier /= 0)
4621             ier = NF90_GET_VAR (forcing_id, vid, &
4622                  &              precip_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
4623                  &            start=start(1:ndim), count=count_force(1:ndim))
4624             a_er = a_er.OR.(ier /= 0)
4625
4626             ndim = 3;
4627             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4628             count_force(1:ndim) = SHAPE(gpp_daily_fm_g)
4629             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4630             ier = NF90_INQ_VARID (forcing_id,'gpp',vid)
4631             a_er = a_er.OR.(ier /= 0)
4632             ier = NF90_GET_VAR (forcing_id, vid, &
4633                  &            gpp_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4634                  &            start=start(1:ndim), count=count_force(1:ndim))
4635             a_er = a_er.OR.(ier /= 0)
4636
4637             ndim = 3;
4638             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4639             count_force(1:ndim) = SHAPE(veget_fm_g)
4640             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4641             ier = NF90_INQ_VARID (forcing_id,'veget',vid)
4642             a_er = a_er.OR.(ier /= 0)
4643             ier = NF90_GET_VAR (forcing_id, vid, &
4644                  &            veget_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4645                  &            start=start(1:ndim), count=count_force(1:ndim))
4646             a_er = a_er.OR.(ier /= 0)
4647
4648             ndim = 3;
4649             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4650             count_force(1:ndim) = SHAPE(veget_max_fm_g)
4651             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4652             ier = NF90_INQ_VARID (forcing_id,'veget_max',vid)
4653             a_er = a_er.OR.(ier /= 0)
4654             ier = NF90_GET_VAR (forcing_id, vid, &
4655                  &            veget_max_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4656                  &            start=start(1:ndim), count=count_force(1:ndim))
4657             a_er = a_er.OR.(ier /= 0)
4658
4659             ndim = 3;
4660             start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
4661             count_force(1:ndim) = SHAPE(lai_fm_g)
4662             count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
4663             ier = NF90_INQ_VARID (forcing_id,'lai',vid)
4664             a_er = a_er.OR.(ier /= 0)
4665             ier = NF90_GET_VAR (forcing_id, vid, &
4666                  &            lai_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
4667                  &            start=start(1:ndim), count=count_force(1:ndim))
4668             a_er = a_er.OR.(ier /= 0)
4669             IF (a_er) THEN
4670                CALL ipslerr_p (3,'stomate_forcing_read', &
4671                     &        'PROBLEM when read forcing file', &
4672                     &        '','')
4673             ENDIF
4674
4675          ENDIF ! (ifirst(iblocks) /= ilast(iblocks))
4676       ENDDO ! iblocks
4677    ENDIF ! is_root_prc
4678
4679  !! 4. Distribute the variable over several processors
4680
4681    CALL scatter(clay_fm_g,clay_fm)
4682    CALL scatter(humrel_daily_fm_g,humrel_daily_fm)
4683    CALL scatter(litterhum_daily_fm_g,litterhum_daily_fm)
4684    CALL scatter(t2m_daily_fm_g,t2m_daily_fm)
4685    CALL scatter(t2m_min_daily_fm_g,t2m_min_daily_fm)
4686    CALL scatter(tsurf_daily_fm_g,tsurf_daily_fm)
4687    CALL scatter(tsoil_daily_fm_g,tsoil_daily_fm)
4688    CALL scatter(soilhum_daily_fm_g,soilhum_daily_fm)
4689    CALL scatter(precip_fm_g,precip_fm)
4690    CALL scatter(gpp_daily_fm_g,gpp_daily_fm)
4691    CALL scatter(veget_fm_g,veget_fm)
4692    CALL scatter(veget_max_fm_g,veget_max_fm)
4693    CALL scatter(lai_fm_g,lai_fm)
4694 
4695  END SUBROUTINE stomate_forcing_read
4696
4697
4698!! ================================================================================================================================
4699!! SUBROUTINE   : setlai
4700!!
4701!>\BRIEF        Routine to force the lai in STOMATE. The code in this routine
4702!! simply CALCULATES lai and is therefore not functional. The routine should be
4703!! rewritten if one wants to force lai.
4704!!
4705!! DESCRIPTION  : None
4706!!
4707!! RECENT CHANGE(S) : None
4708!!
4709!! MAIN OUTPUT VARIABLE(S): ::lai
4710!!
4711!! REFERENCE(S) : None
4712!!
4713!! FLOWCHART : None
4714!! \n
4715!_ ================================================================================================================================
4716 
4717  SUBROUTINE setlai(npts,lai)
4718
4719  !! 0 Variable and parameter declaration
4720 
4721    !! 0.1 Input variables
4722
4723    INTEGER(i_std),INTENT(in)                    :: npts !! Domain size - number of pixels (unitless)
4724   
4725    !! 0.2 Output variables
4726
4727    REAL(r_std),DIMENSION(npts,nvm),INTENT(out)  :: lai  !! PFT leaf area index @tex $(m^{2} m^{-2})$ @endtex
4728
4729    !! 0.3 Modified variables
4730
4731    !! 0.4 Local variables
4732
4733    INTEGER(i_std)                               :: j    !! index (unitless)
4734!_ ================================================================================================================================
4735   
4736    !! 1. Set lai for bare soil to zero
4737
4738    lai(:,ibare_sechiba) = zero
4739
4740    !! 2. Multiply foliage biomass by sla to calculate lai for all PFTs and pixels
4741
4742    DO j=2,nvm
4743       lai(:,j) = biomass(:,j,ileaf,icarbon)*sla_calc(:,j)
4744    ENDDO
4745   
4746  END SUBROUTINE setlai
4747
4748END MODULE stomate
Note: See TracBrowser for help on using the repository browser.