source: branches/publications/ORCHIDEE_v2_r5968/src_stomate/stomate.f90 @ 8593

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

Modification on co2_to_bm and respiration maintanence. See ticket #419

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