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

Last change on this file since 8418 was 8418, checked in by bertrand.guenet, 5 months ago

The Moyano function describing the soil moisture effect on OM decomposition is added

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