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

Last change on this file since 6151 was 5573, checked in by josefine.ghattas, 6 years ago

Do not use t2m/q2m coming from the atmospheric model anymore and instead use temp_air and qair from first atmpospheric model layer. t2m/q2m were previously used only in diffuco_trans_co2 and in stomate_initialize. Also removed diagnostic variables t2m and q2m. See ticket #372

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