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

Last change on this file since 7326 was 7326, checked in by josefine.ghattas, 3 years ago

Corrected bug on carbon balance closure. See ticket #785
Integration in branch 2_2 done by P. Cadule

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