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

Last change on this file since 8398 was 4183, checked in by josefine.ghattas, 7 years ago

Adapted to avoid warning messages when compiling with ifort debug options "check all". See ticket #254

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