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

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

Integrated correction of diagnostic variables rhSoil and rhLitter as done in the trunk rev [6362]. See ticket #630

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