source: branches/publications/ORCHIDEE-Clateral/src_stomate/stomate_litter.f90 @ 7346

Last change on this file since 7346 was 7191, checked in by haicheng.zhang, 3 years ago

Implementing latral transfers of sediment and POC from land to ocean through river in ORCHILEAK

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 88.5 KB
Line 
1! =================================================================================================================================
2! MODULE       : stomate_litter
3!
4! CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
5!
6! LICENCE      : IPSL (2006)
7! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF       Update litter and lignine content after litter fall and
10!! calculating litter decomposition.     
11!!
12!!\n DESCRIPTION: None
13!!
14!! RECENT CHANGE(S): None
15!!
16!! REFERENCE(S) : None
17!!
18!! SVN          :
19!! $HeadURL$
20!! $Date$
21!! $Revision$
22!! \n
23!_ ================================================================================================================================
24
25MODULE stomate_litter
26
27  ! modules used:
28
29  USE ioipsl_para
30  USE stomate_data
31  USE constantes
32  USE constantes_soil
33  USE pft_parameters
34
35  IMPLICIT NONE
36
37  ! private & public routines
38
39  PRIVATE
40  PUBLIC littercalc,littercalc_clear, deadleaf
41
42
43
44  LOGICAL, SAVE                        :: firstcall = .TRUE.       !! first call
45!$OMP THREADPRIVATE(firstcall)
46
47CONTAINS
48
49!! ================================================================================================================================
50!! SUBROUTINE   : littercalc_calc
51!!
52!!\BRIEF        Set the flag ::firstcall to .TRUE. and as such activate section
53!! 1.1 of the subroutine littercalc (see below).
54!!
55!! DESCRIPTION  : None
56!!
57!! RECENT CHANGE(S) : None
58!!
59!! MAIN OUTPUT VARIABLE(S) : None
60!!
61!! REFERENCE(S) : None
62!!
63!! FLOWCHART : None
64!! \n
65!_ ================================================================================================================================
66
67  SUBROUTINE littercalc_clear
68    firstcall =.TRUE.
69  END SUBROUTINE littercalc_clear
70
71
72!! ================================================================================================================================
73!! SUBROUTINE   : littercalc
74!!
75!!\BRIEF        Calculation of the litter decomposition and therefore of the
76!! heterotrophic respiration from litter following Parton et al. (1987).
77!!
78!! DESCRIPTION  : The littercal routine splits the litter in 4 pools:
79!! aboveground metaboblic, aboveground structural, belowground metabolic and
80!! belowground structural. the fraction (F) of plant material going to metabolic
81!! and structural is defined following Parton et al. (1987). The below ground
82!! litter is discretized over 11 layers (the same used by the CWRR hydrology
83!! scheme).
84!! \latexonly
85!! \input{littercalc1.tex}
86!! \endlatexonly
87!! \n
88!! where L is the lignin content of the plant carbon pools considered and CN
89!! its CN ratio. L and CN are fixed parameters for each plant carbon pools,
90!! therefore it is the ratio between each plant carbon pool within a PFT, which
91!! controlled the part of the total litter, that will be considered as
92!! recalcitrant (i.e. structural litter) or labile (i.e. metabolic litter).\n 
93!!
94!! The routine calculates the fraction of aboveground litter which is metabolic
95!! or structural (the litterpart variable) which is then used in lpj_fire.f90.\n
96!!
97!! In the section 2, the routine calculate the new plant material entering the
98!! litter pools by phenological death of plants organs (corresponding to the
99!! variable turnover) and by fire, herbivory and others non phenological causes
100!! (variable bm_to_litter). This calculation is first done for each PFT and then
101!! the values calculated for each PFT are added up. Following the same approach
102!! the lignin content of the total structural litter is calculated and will be
103!! then used as a factor control of the decomposition of the structural litter
104!! (lignin_struc) in the section 5.1.2. A test is performed to avoid that we add
105!! more lignin than structural litter. Finally, the variable litterpart is
106!! updated.\n
107!!
108!! In the section 3 and 4 the temperature and the moisture controlling the
109!! decomposition are calculated for above and belowground. For aboveground
110!! litter, air temperature and litter moisture are calculated in sechiba and used
111!! directly. For belowground, soil temperature and moisture are also calculated
112!! in sechiba but are modulated as a function of the soil depth. The modulation
113!! is a multiplying factor exponentially distributed between 0 (in depth) and 1
114!! in surface.\n 
115!!
116!! Then, in the section 5 and 6, the routine calculates the structural litter decomposition
117!! (C) following first order kinetics following Parton et al. (1987).
118!! \latexonly
119!! \input{littercalc2.tex}
120!! \endlatexonly
121!! \n
122!! with k the decomposition rate of the structural litter.
123!! k corresponds to
124!! \latexonly
125!! \input{littercalc3.tex}
126!! \endlatexonly
127!! \n
128!! with littertau the turnover rate, T a function of the temperature and M a function of
129!! the moisture described below.\n
130!! 
131!! Then, the fraction of dead leaves (DL) composed by aboveground structural litter is
132!! calculated as following
133!! \latexonly
134!! \input{littercalc4.tex}
135!! \endlatexonly
136!! \n
137!! with k the decomposition rate of the structural litter previously
138!! described.\n
139!!
140!! In the section 5.1 and 6.1, the fraction of decomposed structural litter
141!! incorporated to the soil (Input) and its associated heterotrophic respiration are
142!! calculated. For structural litter, the C decomposed could go in the active
143!! soil carbon pool or in the slow carbon, as described in
144!! stomate_soilcarbon.f90.\n
145!! \latexonly
146!! \input{littercalc5.tex}
147!! \endlatexonly
148!! \n
149!! with f a parameter describing the fraction of structural litter incorporated
150!! into the considered soil carbon pool, C the amount of litter decomposed and L
151!! the amount of lignin in the litter. The litter decomposed which is not
152!! incorporated into the soil is respired.\n
153!!
154!! In the section 5.2 and 6.2, the fraction of decomposed metabolic litter
155!! incorporated to the soil and its associated heterotrophic respiration are
156!! calculated with the same approaches presented for 5.1 but no control factor
157!! depending on the lignin content are used.\n
158!!
159!! In the section 7 the dead leaf cover is calculated through a call to the
160!! deadleaf subroutine presented below.\n
161!!
162!! MAIN OUTPUT VARIABLES: ::deadleaf_cover, ::resp_hetero_litter,
163!! ::control_temp, ::control_moist
164!!
165!! REFERENCES:
166!! - Parton, WJ, Schimel, DS, Cole, CV, and Ojima, DS. 1987. Analysis
167!! of factors controlling soil organic matter levels in Great Plains
168!! grasslands. Soil Science Society of America journal (USA)
169!! (51):1173-1179.
170!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
171!! Environmental Modelling & Software (2011), doi:10.1016|j.envsoft.2011.05.016
172!!
173!! FLOWCHART    :
174!! \latexonly
175!! \includegraphics(scale=0.5){littercalcflow.jpg}
176!! \endlatexonly
177!! \n
178!_ ================================================================================================================================
179
180  SUBROUTINE littercalc (npts, dt, &
181       turnover, bm_to_litter, Cinp_manure_solid, &
182       veget_max, tsurf, tsoil, soilhum, litterhum, rprof, &
183       litterpart, litter_above, litter_below, dead_leaves, &
184       lignin_struc_above, lignin_struc_below, &
185       deadleaf_cover, resp_hetero_litter, resp_hetero_flood, &
186       control_temp_above, control_temp_soil, &
187       control_moist_above, control_moist_soil, &
188       litter_mc,soilcarbon_input, floodcarbon_input, soil_mc, soiltile, &
189       clay, bulk_dens, soil_ph, poor_soils, carbon, flood_frac)
190
191    !! 0. Variable and parameter declaration
192   
193    !! 0.1 Input variables
194
195    INTEGER(i_std), INTENT(in)                                  :: npts                !! Domain size - number of grid pixels
196    REAL(r_std), INTENT(in)                                     :: dt                  !! Time step of the simulations for stomate
197                                                                                       !! @tex $(dtradia one\_day^{-1})$ @endtex
198    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: turnover          !! Turnover rates of plant biomass
199                                                                                       !! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
200    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: bm_to_litter      !! Conversion of biomass to litter
201                                                                                       !! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
202    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: veget_max          !! PFT "Maximal" coverage fraction of a PFT
203                                                                                       !! defined in the input vegetation map
204                                                                                       !! @tex $(m^2 m^{-2})$ @endtex
205        REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: Cinp_manure_solid  !! Solid manure-C input (metabolic litter-C, gC.m^{-2} dt{-1})
206                                                                                           !! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex                                                                                                                 
207    REAL(r_std), DIMENSION(npts), INTENT(in)                     :: tsurf              !! Temperature (K) at the surface
208    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)                :: tsoil              !! Soil temperature (K)
209    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)                :: soilhum            !! Daily soil humidity of each soil layer
210                                                                                       !! (unitless)
211    REAL(r_std), DIMENSION(npts), INTENT(in)                     :: litterhum          !! Daily litter humidity (unitless)
212    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                 :: rprof              !! PFT root depth as calculated in stomate.f90 from parameter
213                                                                                       !! humcste which is root profile for different PFTs
214                                                                                       !! in slowproc.f90 (m)
215    REAL(r_std),DIMENSION (npts,nstm), INTENT(in)                :: litter_mc          !! soil moisture content \f($m^3 \times m^3$)\f
216    REAL(r_std),DIMENSION (npts,nslm,nstm), INTENT(in)           :: soil_mc            !! soil moisture content \f($m^3 \times m^3$)\f
217    REAL(r_std),DIMENSION (npts,nstm), INTENT (in)               :: soiltile           !! Fraction of each soil tile (0-1, unitless)
218    REAL(r_std), DIMENSION(npts), INTENT(in)                     :: clay               !! Clay fraction (unitless, 0-1)
219    REAL(r_std), DIMENSION(npts), INTENT(inout)                  :: bulk_dens          !! Variable local of bulk density for the moment must change in the futur (kg m-3)
220
221    REAL(r_std), DIMENSION(npts), INTENT(in)                     :: soil_ph            !! soil pH (pH unit, 0-14)
222    REAL(r_std), DIMENSION(npts), INTENT(in)                     :: poor_soils         !! Fraction of poor soils (0-1)
223    REAL(r_std), DIMENSION(npts,ncarb,nvm,nslmd), INTENT(in)     :: carbon             !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f
224    REAL(r_std), DIMENSION(npts), INTENT(in)                     :: flood_frac         !! flooded fraction, needed to calculate heterotrophic respiration input to floodplain
225    !! 0.2 Output variables
226   
227    REAL(r_std), DIMENSION(npts), INTENT(out)                    :: deadleaf_cover     !! Fraction of soil covered by dead leaves
228                                                                                       !! over all PFTs (0-1, unitless)
229    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                :: resp_hetero_litter !! Litter heterotrophic respiration. The unit
230                                                                                       !! is given by m^2 of ground. 
231                                                                                       !! @tex $(gC dtradia one\_day^{-1}) m^{-2})$ @endtex
232    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                :: resp_hetero_flood  !! Litter heterotrophic respiration in flooded areas.
233                                                                                       !! The unit is given by m^2 of ground.
234                                                                                       !! @tex $(gC dtradia one\_day^{-1}) m^{-2})$ @endtex
235    REAL(r_std), DIMENSION(npts,nslmd,npool*2), INTENT(out)       :: control_temp_soil  !! Temperature control of heterotrophic
236                                                                                       !! respiration belowground,(0-1, unitless)
237    REAL(r_std), DIMENSION(npts,nslmd,nvm), INTENT(out)           :: control_moist_soil !! Moisture control of heterotrophic
238                                                                                       !! respiration aboveground(0.25-1,unitless)
239    REAL(r_std), DIMENSION(npts,nlitt), INTENT(out)              :: control_temp_above !! Temperature control of heterotrophic 
240                                                                                       !! respiration, above (0-1, 
241                                                                                       !! unitless)
242    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                :: control_moist_above!! Moisture control of heterotrophic 
243                                                                                       !! respiration aboveground(0.25-1, unitless)
244
245    !! 0.3 Modified variables
246   
247    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)       :: litterpart         !! Fraction of litter above the ground
248                                                                                      !! belonging to the different PFTs (0-1,
249                                                                                      !! unitless)
250    REAL(r_std), DIMENSION(npts,nlitt,nvm,nelements), INTENT(inout) :: litter_above   !! Metabolic and structural litter,above ground
251                                                                                      !! The unit is given by m^2 of
252                                                                                      !! ground @tex $(gC m^{-2})$ @endtex
253    REAL(r_std), DIMENSION(npts,nlitt,nvm,nslmd,nelements), INTENT(inout) ::litter_below  !! Metabolic and structural litter, below ground
254                                                                                      !! The unit is given by m^2 of
255                                                                                      !! ground @tex $(gC m^{-2} $ @endtex
256    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)       :: dead_leaves        !! Dead leaves per ground unit area, per PFT,
257                                                                                      !! metabolic and structural in
258                                                                                      !! @tex $(gC m^{-2})$ @endtex
259    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: lignin_struc_above  !! Ratio Lignin content in structural litter,
260                                                                                      !! above ground, (0-1, unitless)
261    REAL(r_std), DIMENSION(npts,nvm,nslmd), INTENT(inout)       :: lignin_struc_below  !! Ratio Lignin content in structural litter,
262                                                                                      !! below ground, (0-1, unitless)
263    REAL(r_std), DIMENSION(npts,nvm,nslmd,npool,nelements), INTENT(inout) :: soilcarbon_input !! Dissolved Organic Carbon input from litter decomposition
264                                                                                      !! The unit is given by m^2
265                                                                                      !! @tex $(gC m^{-2})$ @endtex
266    REAL(r_std), DIMENSION(npts,nvm,npool,nelements), INTENT(inout) :: floodcarbon_input !! Dissolved Organic Carbon input from litter
267                                                                                      !! decomposition in flooded areas
268                                                                                      !! The unit is given by m^2
269                                                                                      !! @tex $(gC m^{-2})$ @endtex
270    !! 0.4 Local variables
271
272    REAL(r_std), DIMENSION(npts,nvm)                            :: control_flood_above!! Moisture control of heterotrophic
273                                                                                      !! respiration aboveground(0.25-1, unitless)
274!$OMP THREADPRIVATE(control_flood_above)
275    REAL(r_std), SAVE, DIMENSION(nparts,nlitt)                  :: litterfrac         !! The fraction of leaves, wood, etc. that
276                                                                                      !! goes into metabolic and structural
277                                                                                      !! litterpools (0-1, unitless)
278    REAL(r_std), DIMENSION(npts,nparts,nlitt)                   :: litterfrac_pxl     !! The fraction of leaves, wood, etc. that
279                                                                                      !! goes into metabolic and structural
280                                                                                      !! litterpools (0-1, unitless)
281!$OMP THREADPRIVATE(litterfrac)
282    REAL(r_std),DIMENSION(0:nslm)                               :: z_soil             !! Soil levels (m)
283!$OMP THREADPRIVATE(z_soil)
284                                                                                      !! profiles (unitless)
285    REAL(r_std), SAVE, DIMENSION(nlitt)                         :: litter_tau         !! Turnover time in litter pools (days)
286!$OMP THREADPRIVATE(litter_tau)
287    REAL(r_std), SAVE, DIMENSION(npool)                         :: pool_tau         !! Turnover time in litter and soil carbon pools (days)
288!$OMP THREADPRIVATE(pool_tau)
289    REAL(r_std), SAVE, DIMENSION(npool)                         :: DOC_tau          !! Residence time of DOC (days)   
290!$OMP THREADPRIVATE(DOC_tau)
291    REAL(r_std), SAVE, DIMENSION(nlitt,ncarb,nlevs)             :: frac_soil          !! Fraction of litter that goes into soil
292                                                                                      !! (litter -> carbon, above and below). The
293                                                                                      !! remaining part goes to the atmosphere
294!$OMP THREADPRIVATE(frac_soil)
295    REAL(r_std), DIMENSION(npts)                                :: tsoil_decomp       !! Temperature used for decompostition in
296                                                                                      !! soil (K)
297    REAL(r_std), DIMENSION(npts)                                :: soilhum_decomp     !! Humidity used for decompostition in soil
298                                                                                      !! (unitless)
299 !   REAL(r_std), SAVE, DIMENSION(npts)                          :: z_lit              !! Thickness of the above ground litter layer (mm)
300!!$OMP THREADPRIVATE(z_lit)
301
302    REAL(r_std), DIMENSION(npts)                                :: fd                 !! Fraction of structural or metabolic litter
303                                                                                      !! decomposed (unitless)
304    REAL(r_std), DIMENSION(npts,nelements)                      :: qd                 !! Quantity of structural or metabolic litter
305                                                                                      !! decomposed @tex $(gC m^{-2})$ @endtex
306    REAL(r_std), DIMENSION(npts,nelements)                      :: qd_flood           !! Quantity of structural or metabolic litter
307                                                                                      !! decomposed in flooded areas @tex $(gC m^{-2})$ @endtex
308    REAL(r_std), DIMENSION(npts,nvm)                            :: old_struc_above    !! Old structural litter, above ground
309                                                                                      !! @tex $(gC m^{-2})$ @endtex
310    REAL(r_std), DIMENSION(npts,nvm,nslmd)                      :: old_struc_below    !! Old structural litter, below ground
311                                                                                      !! @tex $(gC m^{-2})$ @endtex
312    REAL(r_std), DIMENSION(npts,nvm,nlitt, nelements)           :: litter_inc_PFT_above  !! Increase of litter, per PFT, metabolic and
313                                                                                      !! structural, above ground. The
314                                                                                      !! unit is given by m^2 of ground. 
315                                                                                      !! @tex $(gC m^{-2})$ @endtex
316    REAL(r_std), DIMENSION(npts,nvm,nlitt,nslm,nelements)      :: litter_inc_PFT_below     !! Increase of litter, per PFT, metabolic and
317                                                                                      !! structural, above ground. The
318                                                                                      !!unit is given by m^2 of ground.
319                                                                                      !! @tex $(gC m^{-2})$ @endtex
320    REAL(r_std), DIMENSION(npts,nlitt,nvm,nelements)            :: litter_inc_above   !! Increase of metabolic and structural
321                                                                                      !! litter, above ground. The unit
322                                                                                      !! is given by m^2 of ground.
323                                                                                      !! @tex $(gC m^{-2})$ @endtex
324    REAL(r_std), DIMENSION(npts,nlitt,nvm,nslm,nelements)       :: litter_inc_below   !! Increase of metabolic and structural
325                                                                                      !! litter below ground. The unit
326                                                                                      !! is given by m^2 of ground
327                                                                                      !! @tex $(gC m^{-2})$ @endtex
328
329    REAL(r_std), DIMENSION(npts,nvm)                            :: lignin_struc_inc_above !! Lignin increase in structural litter,
330                                                                                      !! above ground. The unit is given
331                                                                                      !! by m^2 of ground.
332                                                                                      !! @tex $(gC m^{-2})$ @endtex
333    REAL(r_std), DIMENSION(npts,nvm,nslm)                       :: lignin_struc_inc_below   !! Lignin increase in structural litter,
334                                                                                      !! below ground. The unit is given
335                                                                                      !! by m^2 of ground
336                                                                                      !! @tex $(gC m^{-2})$ @endtex
337    REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements)            :: litter_pft         !! Metabolic and structural litter above the
338                                                                                      !! ground per PFT
339    REAL(r_std), DIMENSION(npts)                                :: zdiff_min          !! Intermediate field for looking for minimum
340                                                                                      !! of what? this is not used in the code.
341                                                                                      !! [??CHECK] could we delete it?
342    CHARACTER(LEN=10), DIMENSION(nlitt)                         :: litter_str         !! Messages to write output information about
343                                                                                      !! the litter
344    CHARACTER(LEN=22), DIMENSION(nparts)                        :: part_str           !! Messages to write output information about
345                                                                                      !! the plant
346    CHARACTER(LEN=7), DIMENSION(ncarb)                          :: carbon_str         !! Messages to write output information about
347                                                                                      !! the soil carbon
348    CHARACTER(LEN=5), DIMENSION(nlevs)                          :: level_str          !! Messages to write output information about
349                                                                                      !! the level (aboveground or belowground litter)
350    INTEGER(i_std)                                              :: i,j,k,l,ig,m          !! Indices (unitless)
351    REAL(r_std), DIMENSION(npts)                                :: Dif_coef           !! Diffusion coeficient used for bioturbation (m2 days-1)
352!!$OMP THREADPRIVATE(Dif_coef)
353
354    REAL(r_std), DIMENSION(npts,nlitt,nvm,nslm,nelements)       :: litter_flux        !! Belowground litter carbon flux within pools towards the soil column\f$(gC m^{2} dt^{-1})$\f
355    REAL(r_std), DIMENSION(npts,nlitt,nvm,nslm,nelements)       :: litter_flux_old    !! Belowground litter carbon flux within pools towards the soil column\f$(gC m^{2} dt^{-1})$\f
356                                                                                      !! used for storage
357    REAL(r_std), DIMENSION(npts,nlitt,nvm,nelements)            :: litter_above_flux  !! Above litter carbon flux within pools towards the soil column\f$(gC m^{2} dt^{-1})$\f
358    REAL(r_std), DIMENSION(npts,nlitt,nvm,nelements)            :: litter_above_flux_old  !! Above litter carbon flux within pools towards the soil column\f$(gC m^{2} dt^{-1})$\f
359                                                                                      !! used for storage
360    REAL(r_std), DIMENSION(npts)                                :: rpc                !! Scaling factor for integrating vertical soil
361                                                                                      !! profiles (unitless)
362    REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements)     :: check_intern       !! Contains the components of the internal
363                                                                                      !! mass balance chech for this routine
364                                                                                      !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
365    REAL(r_std), DIMENSION(npts,nvm,nelements)             :: closure_intern     !! Check closure of internal mass balance
366
367    REAL(r_std), DIMENSION(npts,nvm,nelements)             :: pool_start         !! Start and end pool of this routine
368                                                                                      !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
369    REAL(r_std), DIMENSION(npts,nvm,nelements)             :: pool_end           !! End pool of this routine not corrected
370                                                                                      !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
371    REAL(r_std), DIMENSION(npts,nvm,nelements)             :: pool_end_after           !! Start and end pool of this routine
372                                                                                      !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
373    REAL(r_std), DIMENSION(npts,nslm)                      :: moist_soil         !! moisture in soil for one pixel (m3/m3)
374
375   ! REAL(r_std), DIMENSION(npts,nvm,nslm)                       :: resp_hetero_litter_layer !! Litter heterotrophic respiration per layer. The unit
376                                                                                      !! is given by m^2 of ground by z. 
377                                                                                      !! @tex $(gC dtradia one\_day^{-1}) m^{-2} z^{-1})$ @endtex
378    REAL(r_std), DIMENSION(npts,nlitt,nvm,nslmd,nelements)  ::litter_below_old         !! Metabolic and structural litter, below ground
379                                                                                      !! The unit is given by m^2 of
380                                                                                      !! ground @tex $(gC m^{-2} $ @endtex used for storage
381    REAL(r_std), DIMENSION(npts,nlitt,nvm,nslmd,nelements)  ::litter_below_old_buffer  !! Metabolic and structural litter, below ground
382                                                                                      !! The unit is given by m^2 of
383                                                                                      !! ground @tex $(gC m^{-2} $ @endtex used for storage
384    REAL(r_std), DIMENSION(npts,nlitt,nvm,nelements)       :: litter_above_old        !! Metabolic and structural litter,above ground
385                                                                                      !! The unit is given by m^2 of
386                                                                                      !! ground @tex $(gC m^{-2})$ @endtex used for storage
387    REAL(r_std), DIMENSION(npts)                           :: one_array 
388!_ ================================================================================================================================
389
390    IF (printlev>=3) WRITE(numout,*) 'Entering littercalc', nslmd
391
392  !! 1. Initialisations of the different fields during the first call of the routine
393    !! 1.0 Calculation of soil moisture
394    moist_soil(:,:) = zero
395      DO l = 1,nslm
396          DO i = 1, nstm
397             moist_soil(:,l)= moist_soil(:,l) + soil_mc(:,l,i)*soiltile(:,i)
398          ENDDO
399      ENDDO
400
401    WHERE (bulk_dens(:) .LT. 500)
402                bulk_dens(:) = bulk_dens(:)*1e3
403    ENDWHERE
404
405    z_soil(0) = zero
406    z_soil(1:nslm) = zlt(1:nslm)
407
408    !! 1.1 Initialize check for mass balance closure
409    !  The mass balance is calculated at the end of this routine
410    !  in section 14. Initial biomass and harvest pool and all other
411    !  relevant pools were set to zero before calling this routine.
412    pool_start(:,:,:) = zero
413    IF (ld_doc) THEN
414       DO m = 2, nvm
415          DO i = 1,nlitt
416             DO l = 1, nslmd
417                pool_start(:,m,icarbon) = pool_start(:,m,icarbon) + &
418                     ( litter_below(:,i,m,l,icarbon) ) * veget_max(:,m)
419             ENDDO
420          ENDDO
421       ENDDO
422
423       DO m = 2, nvm
424          DO i = 1,nlitt
425             pool_start(:,m,icarbon) = pool_start(:,m,icarbon) + &
426                  ( litter_above(:,i,m,icarbon) ) * veget_max(:,m)
427          ENDDO
428       ENDDO
429
430       DO m = 2, nvm
431          DO i = 1, nparts
432             pool_start(:,m,icarbon) = pool_start(:,m,icarbon) + &
433                  (turnover(:,m,i,icarbon) + bm_to_litter(:,m,i,icarbon)) * veget_max(:,m)
434          ENDDO
435       ENDDO
436    ENDIF
437    IF ( firstcall ) THEN
438
439       !! 1.2.1 litter fractions:
440       !!   what fraction of leaves, wood, etc. goes into metabolic and structural litterpools
441       DO k = 1, nparts
442
443          litterfrac(k,imetabolic) = metabolic_ref_frac - metabolic_LN_ratio * LC(k) * CN(k)
444          litterfrac(k,istructural) = un - litterfrac(k,imetabolic)
445   
446       ENDDO
447
448       !! 1.2.2 residence times in litter pools (days)
449       litter_tau(imetabolic) = tau_metabolic * one_year      ! .5 years
450       litter_tau(istructural) = tau_struct * one_year     ! 3 years
451
452       pool_tau(imetabo) = tau_metabolic 
453       pool_tau(istrabo) = tau_struct 
454       pool_tau(imetbel) = tau_metabolic 
455       pool_tau(istrbel) = tau_struct 
456       pool_tau(iact) = carbon_tau_iactive 
457       pool_tau(islo) = carbon_tau_islow 
458       pool_tau(ipas) = carbon_tau_ipassive 
459       DOC_tau(imetabo) = DOC_tau_labile / one_year
460       DOC_tau(istrabo) = DOC_tau_labile / one_year
461       DOC_tau(imetbel) = DOC_tau_labile / one_year
462       DOC_tau(istrbel) = DOC_tau_labile / one_year
463       DOC_tau(iact) = DOC_tau_labile / one_year
464       DOC_tau(islo) = DOC_tau_stable / one_year
465       DOC_tau(ipas) = DOC_tau_stable / one_year
466
467       !! 1.2.3 decomposition flux fraction that goes into soil
468       !       (litter -> carbon, above and below)
469       !       1-frac_soil goes into atmosphere
470       frac_soil(:,:,:) = zero
471
472       ! structural litter: lignin fraction goes into slow pool + respiration,
473       !                    rest into active pool + respiration
474       frac_soil(istructural,iactive,iabove) = frac_soil_struct_aa
475       frac_soil(istructural,iactive,ibelow) = frac_soil_struct_ab
476       frac_soil(istructural,islow,iabove) = frac_soil_struct_sa
477       frac_soil(istructural,islow,ibelow) = frac_soil_struct_sb
478
479       ! metabolic litter: all goes into active pool + respiration.
480       !   Nothing into slow or passive pool.
481       frac_soil(imetabolic,iactive,iabove) = frac_soil_metab_aa
482       frac_soil(imetabolic,iactive,ibelow) = frac_soil_metab_ab
483       
484 
485       !! 1.4 messages
486       litter_str(imetabolic) = 'metabolic'
487       litter_str(istructural) = 'structural'
488
489       carbon_str(iactive) = 'active'
490       carbon_str(islow) = 'slow'
491       carbon_str(ipassive) = 'passive'
492
493       level_str(iabove) = 'above'
494       level_str(ibelow) = 'below'
495
496       part_str(ileaf) = 'leaves'
497       part_str(isapabove) = 'sap above ground'
498       part_str(isapbelow) = 'sap below ground'
499       part_str(iheartabove) = 'heartwood above ground'
500       part_str(iheartbelow) = 'heartwood below ground'
501       part_str(iroot) = 'roots'
502       part_str(ifruit) = 'fruits'
503       part_str(icarbres) = 'carbohydrate reserve'
504
505       WRITE(numout,*) 'litter:'
506
507       WRITE(numout,*) '   > C/N ratios: '
508       DO k = 1, nparts
509          WRITE(numout,*) '       ', part_str(k), ': ',CN(k)
510       ENDDO
511
512       WRITE(numout,*) '   > Lignine/C ratios: '
513       DO k = 1, nparts
514          WRITE(numout,*) '       ', part_str(k), ': ',LC(k)
515       ENDDO
516
517       WRITE(numout,*) '   > fraction of compartment that goes into litter: '
518       DO k = 1, nparts
519          DO m = 1, nlitt
520             WRITE(numout,*) '       ', part_str(k), '-> ',litter_str(m), ':',litterfrac(k,m)
521          ENDDO
522       ENDDO
523
524       WRITE(numout,*) '   > minimal carbon residence time in litter pools (d):'
525       DO m = 1, nlitt
526          WRITE(numout,*) '       ',litter_str(m),':',litter_tau(m)
527       ENDDO
528
529       WRITE(numout,*) '   > litter decomposition flux fraction that really goes '
530       WRITE(numout,*) '     into carbon pools (rest into the atmosphere):'
531       DO m = 1, nlitt
532          DO l = 1, nlevs
533             DO k = 1, ncarb
534                WRITE(numout,*) '       ',litter_str(m),' ',level_str(l),' -> ',&
535                     carbon_str(k),':', frac_soil(m,k,l)
536             ENDDO
537          ENDDO
538       ENDDO
539
540       firstcall = .FALSE.
541
542    ENDIF
543
544    !! 1.3 Initialization
545
546    DO k = 1, nparts
547
548        litterfrac_pxl(:,k,imetabolic) = MAX(zero, MIN(un, metabolic_ref_frac - metabolic_LN_ratio * LC(k) * CN(k)))! * (un + poor_soils(:))
549        litterfrac_pxl(:,k,istructural) = MAX(zero, MIN(un, un - litterfrac_pxl(:,k,imetabolic)))
550
551    ENDDO
552           
553           
554    litter_above_flux(:,:,:,:) = zero
555    litter_flux(:,:,:,:,:) = zero
556    litter_flux_old(:,:,:,:,:) = zero
557    old_struc_above(:,:) = zero
558    old_struc_below(:,:,:) = zero
559    litter_pft(:,:,:,:) = zero
560    Dif_coef(:) = (Dif/one_year)*dt 
561
562   
563    !! 1.3 litter above the ground per PFT.
564    DO j = 2, nvm ! Loop over # PFTs
565       DO k = 1, nlitt !Loop over litter pool
566          litter_pft(:,j,k,icarbon) = litterpart(:,j,k) * litter_above(:,k,j,icarbon)
567       ENDDO
568    ENDDO
569
570   
571    !! 1.4 set output to zero
572    deadleaf_cover(:) = zero
573    resp_hetero_litter(:,:) = zero
574    resp_hetero_flood(:,:) = zero
575    soilcarbon_input(:,:,:,:,:)= zero
576    floodcarbon_input(:,:,:,:)= zero
577    !resp_hetero_litter_layer(:,:,:) = zero
578
579  !! 2. Add biomass to different litterpools (per m^2 of ground)
580   
581    !! 2.1 first, save old structural litter (needed for lignin fractions).
582    !     above
583       DO m = 2,nvm !Loop over PFTs
584
585          old_struc_above(:,m) = litter_above(:,istructural,m,icarbon)
586
587       ENDDO
588
589    !below
590    DO l = 1, nslmd !Loop over soil levels (below ground)
591       DO m = 2,nvm !Loop over PFTs
592
593          old_struc_below(:,m,l) = litter_below(:,istructural,m,l,icarbon)
594
595       ENDDO
596    ENDDO
597   
598    !! 2.2 update litter, dead leaves, and lignin content in structural litter
599    litter_inc_above(:,:,:,:) = zero
600    lignin_struc_inc_above(:,:) = zero
601    litter_inc_below(:,:,:,:,:) = zero
602    lignin_struc_inc_below(:,:,:) = zero
603    litter_inc_PFT_above(:,:,:,:) = zero
604    litter_inc_PFT_below(:,:,:,:,:) = zero
605
606
607    DO j = 2,nvm !Loop over PFTs
608
609       !! 2.2.1 litter
610       DO k = 1, nlitt    !Loop over litter pools (metabolic and structural)
611
612          !! 2.2.2 calculate litter increase (per m^2 of ground).
613          !       Only a given fracion of fruit turnover is directly coverted into litter.
614          !       Litter increase for each PFT, structural and metabolic, above/below. It must
615          !       be noted that the distribution of above ground litter is done within the soil layers
616          !       as a proportion of the layer size. But for the below ground layer it follows the root profile
617
618          ! For the first litter layer
619          litter_inc_PFT_above(:,j,k,icarbon) = &
620               litterfrac_pxl(:,ileaf,k) * bm_to_litter(:,j,ileaf,icarbon) + &
621               litterfrac_pxl(:,isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
622               litterfrac_pxl(:,iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
623               litterfrac_pxl(:,ifruit,k) * bm_to_litter(:,j,ifruit,icarbon) + &
624               litterfrac_pxl(:,icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
625               litterfrac_pxl(:,ileaf,k) * turnover(:,j,ileaf,icarbon) + &
626               litterfrac_pxl(:,isapabove,k) * turnover(:,j,isapabove,icarbon) + &
627               litterfrac_pxl(:,iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
628               litterfrac_pxl(:,ifruit,k) * turnover(:,j,ifruit,icarbon) + &
629               litterfrac_pxl(:,icarbres,k) * turnover(:,j,icarbres,icarbon) 
630
631             ! litter increase, met/struct, above
632             litter_inc_above(:,k,j,icarbon) = litter_inc_above(:,k,j,icarbon) + litter_inc_PFT_above(:,j,k,icarbon)
633          !! The function used for root profile is coming from stomate_npp.f90
634          rpc(:) = un /( ( un - EXP( -z_soil(nslm) / rprof(:,j) )) ) 
635
636          DO l = 1, nslm
637
638             litter_inc_PFT_below(:,j,k,l,icarbon) = &
639                  litterfrac_pxl(:,isapbelow,k) * bm_to_litter(:,j,isapbelow,icarbon) *  &
640               rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
641                  litterfrac_pxl(:,iheartbelow,k) * bm_to_litter(:,j,iheartbelow,icarbon) * &
642               rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
643                  litterfrac_pxl(:,iroot,k) * bm_to_litter(:,j,iroot,icarbon) * &
644               rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
645                  litterfrac_pxl(:,isapbelow,k) * turnover(:,j,isapbelow,icarbon) * &
646               rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
647                  litterfrac_pxl(:,iheartbelow,k) * turnover(:,j,iheartbelow,icarbon) * &
648               rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
649                  litterfrac_pxl(:,iroot,k) * turnover(:,j,iroot,icarbon) * &
650               rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) )
651
652             ! litter increase, met/struct, below
653             litter_inc_below(:,k,j,l,icarbon) = litter_inc_below(:,k,j,l,icarbon) + litter_inc_PFT_below(:,j,k,l,icarbon)
654          ENDDO
655          !! 2.2.3 dead leaves, for soil cover.
656          dead_leaves(:,j,k) = &
657               dead_leaves(:,j,k) + &
658               litterfrac_pxl(:,ileaf,k) * ( bm_to_litter(:,j,ileaf,icarbon) + turnover(:,j,ileaf,icarbon) )
659
660          !! 2.2.4 lignin increase in structural litter
661          IF ( k .EQ. istructural ) THEN
662
663             lignin_struc_inc_above(:,j) = &
664                  lignin_struc_inc_above(:,j) + &
665                  LC(ileaf) * bm_to_litter(:,j,ileaf,icarbon) + &
666                  LC(isapabove) * bm_to_litter(:,j,isapabove,icarbon) + &
667                  LC(iheartabove) * bm_to_litter(:,j,iheartabove,icarbon) + &
668                  LC(ifruit) * bm_to_litter(:,j,ifruit,icarbon) + &
669                  LC(icarbres) * bm_to_litter(:,j,icarbres,icarbon) + &
670                  LC(ileaf) * turnover(:,j,ileaf,icarbon) + &
671                  LC(isapabove) * turnover(:,j,isapabove,icarbon) + &
672                  LC(iheartabove) * turnover(:,j,iheartabove,icarbon) + &
673                  LC(ifruit) * turnover(:,j,ifruit,icarbon) + &
674                  LC(icarbres) * turnover(:,j,icarbres,icarbon) 
675
676            !! The function used for root profile is coming from stomate_npp.f90
677            DO l = 1, nslm
678
679               lignin_struc_inc_below(:,j,l) = &
680                  lignin_struc_inc_below(:,j,l) + &
681                  LC(isapbelow) * bm_to_litter(:,j,isapbelow,icarbon) * &
682                  rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
683                  LC(iheartbelow) * bm_to_litter(:,j,iheartbelow,icarbon) * &
684                  rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
685                  LC(iroot) * bm_to_litter(:,j,iroot,icarbon) * &
686                  rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
687                  LC(isapbelow)*turnover(:,j,isapbelow,icarbon) * &
688                  rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
689                  LC(iheartbelow)*turnover(:,j,iheartbelow,icarbon) * &
690                  rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) ) + &
691                  LC(iroot)*turnover(:,j,iroot,icarbon) * &
692                  rpc(:) * ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) )
693            ENDDO
694
695          ENDIF
696       ENDDO
697           !! Manure-C input was added to the metabolic & structural litter-C (H. Zhang)
698           litter_inc_above(:,imetabolic,j,icarbon) = litter_inc_above(:,imetabolic,j,icarbon) +Cinp_manure_solid(:,j)*f_met_solmanure
699           litter_inc_above(:,istructural,j,icarbon) = litter_inc_above(:,istructural,j,icarbon) +Cinp_manure_solid(:,j)*(un-f_met_solmanure)                                                                                                                                                             
700     ENDDO
701
702    !! 2.2.5 add new litter (struct/met, above/below)
703        !WRITE(numout,*) 'SUMstomlitt1.0_litter',MINVAL(litter_above),MINVAL(litter_below),MAXVAL(litter_above),MAXVAL(litter_below)
704        !WRITE(numout,*) 'litterin',MINVAL(litter_inc_above),MINVAL(litter_inc_below),MAXVAL(litter_inc_above),MAXVAL(litter_inc_below)
705       
706    litter_above(:,:,:,:) = litter_above(:,:,:,:) + litter_inc_above(:,:,:,:)
707    litter_below(:,:,:,1:nslm,:) = litter_below(:,:,:,1:nslm,:) + litter_inc_below(:,:,:,1:nslm,:)
708                 
709    !WRITE(numout,*) 'SUMstomlitt1.1_litter',SUM(litter_above),SUM(litter_below)
710        !WRITE(numout,*) 'stomlitt_ZHC1',SUM(lignin_struc_below)
711        !WRITE(numout,*) 'litter',MINVAL(litter_above),MINVAL(litter_below),MAXVAL(litter_above),MAXVAL(litter_below)
712        !WRITE(numout,*) 'litterin',MINVAL(litter_inc_above),MINVAL(litter_inc_below),MAXVAL(litter_inc_above),MAXVAL(litter_inc_below)
713    !! 2.2.6 for security: can't add more lignin than structural litter (above/below)
714       DO m = 2,nvm !Lopp over PFTs
715
716          lignin_struc_inc_above(:,m) = &
717               MIN( lignin_struc_inc_above(:,m), litter_inc_above(:,istructural,m,icarbon) )
718
719       ENDDO
720
721    DO l = 1, nslm !Loop over soil levels (below ground)
722       DO m = 2,nvm !Lopp over PFTs
723
724          lignin_struc_inc_below(:,m,l) = &
725               MIN( lignin_struc_inc_below(:,m,l), litter_inc_below(:,istructural,m,l,icarbon) )
726
727       ENDDO
728    ENDDO
729
730    !! 2.2.7.1 new lignin content: add old lignin and lignin increase, divide by
731    !!       total structural litter (above)
732
733       DO m = 2,nvm !Loop over PFTs
734          WHERE( litter_above(:,istructural,m,icarbon) .GT. min_stomate )
735
736       !MM : Soenke modif
737       ! Best vectorization ?
738!!$       lignin_struc(:,:,:) = &
739!!$            ( lignin_struc(:,:,:)*old_struc(:,:,:) + lignin_struc_inc(:,:,:) ) / &
740!!$            litter(:,istructural,:,:,icarbon)
741
742             lignin_struc_above(:,m) = lignin_struc_above(:,m) * old_struc_above(:,m)
743             lignin_struc_above(:,m) = lignin_struc_above(:,m) + lignin_struc_inc_above(:,m)
744             lignin_struc_above(:,m) = lignin_struc_above(:,m) / litter_above(:,istructural,m,icarbon)
745          ELSEWHERE
746             lignin_struc_above(:,m) = zero
747          ENDWHERE
748       ENDDO
749
750    !! 2.2.7.2 new lignin content: add old lignin and lignin increase, divide by
751    !!       total structural litter (below)
752    DO l = 1, nslm !Loop over soil levels (below ground)
753       DO m = 2,nvm !Loop over PFTs
754          WHERE( litter_below(:,istructural,m,l,icarbon) .GT. min_stomate )
755
756       !MM : Soenke modif
757       ! Best vectorization ?
758!!$       lignin_struc(:,:,:) = &
759!!$            ( lignin_struc(:,:,:)*old_struc(:,:,:) + lignin_struc_inc(:,:,:)) / &
760!!$            litter(:,istructural,:,:,icarbon)
761
762             lignin_struc_below(:,m,l) = lignin_struc_below(:,m,l) * old_struc_below(:,m,l)
763             lignin_struc_below(:,m,l) = lignin_struc_below(:,m,l) + lignin_struc_inc_below(:,m,l)
764             lignin_struc_below(:,m,l) = lignin_struc_below(:,m,l) /litter_below(:,istructural,m,l,icarbon)
765          ELSEWHERE
766             lignin_struc_below(:,m,l) = zero
767          ENDWHERE
768       ENDDO
769    ENDDO
770   
771    !! 2.3 new litter fraction per PFT (for structural and metabolic litter, above
772    !!       the ground).
773    DO j = 2,nvm !Loop over PFTs
774      DO l =1,nlitt 
775         WHERE ( litter_above(:,l,j,icarbon) .GT. min_stomate )
776
777          litterpart(:,j,l) = &
778               ( litter_pft(:,j,l,icarbon) + litter_inc_PFT_above(:,j,l,icarbon) ) / litter_above(:,l,j,icarbon)
779
780         ELSEWHERE
781
782          litterpart(:,j,l) = zero
783
784         ENDWHERE
785       ENDDO
786    ENDDO
787  !! 3. Temperature control on decay: Factor between 0 and 1
788
789    !! 3.1 above: surface temperature
790    DO j=1,nlitt
791       control_temp_above(:,j) = control_temp_func (npts, tsurf,litter_tau(j)*(un+poor_soils(:))) 
792    ENDDO         
793    !! 3.2 below: convolution of temperature and decomposer profiles
794    !!            (exponential decomposer profile supposed)
795   
796    !! 3.2.1 integrate over the nslm levels
797    tsoil_decomp(:) = zero
798
799    DO l = 1, nslm
800       DO j = 1,npool
801
802          tsoil_decomp(:) = tsoil(:,l) 
803
804          control_temp_soil(:,l,j) = control_temp_func (npts,tsoil_decomp,pool_tau(j)*(un+poor_soils(:)))
805
806       ENDDO
807    ENDDO
808
809    DO l = 1, nslm
810       DO j = 1,npool
811
812          tsoil_decomp(:) = tsoil(:,l)
813
814          control_temp_soil(:,l,npool+j) = control_temp_func(npts,tsoil_decomp,DOC_tau(j)*(un+poor_soils(:)))
815
816       ENDDO
817    ENDDO
818       
819    control_temp_soil(:,nslmd,:)=control_temp_soil(:,nslm,:)
820        WHERE (ISNAN(control_temp_soil(:,:,:)))
821            control_temp_soil(:,:,:)=un
822        ENDWHERE
823       
824  !! 4. Moisture control. Factor between 0 and 1
825
826    !! 4.1 above the ground: litter humidity
827
828    one_array(:) = un
829    soilhum_decomp(:) = zero
830    DO m= 1,nvm
831       !We used the sum over the 4 first layer as it is done to
832       !compute litterhum in hydrol.f90
833         soilhum_decomp(:) =soil_mc(:,1,pref_soil_veg(m))*(z_soil(1)/z_soil(4)) + &
834                            soil_mc(:,2,pref_soil_veg(m))*((z_soil(2)-z_soil(1))/z_soil(4)) + &
835                            soil_mc(:,3,pref_soil_veg(m))*((z_soil(3)-z_soil(2))/z_soil(4)) + &
836                            soil_mc(:,4,pref_soil_veg(m))*((z_soil(4)-z_soil(3))/z_soil(4))
837         !
838       IF (moist_func_Moyano) THEN
839        control_moist_above(:,m) = control_moist_func_moyano (npts, soilhum_decomp,bulk_dens, &
840                                                    clay, carbon, veget_max)
841       ELSE
842         control_moist_above(:,m) = control_moist_func (npts, soilhum_decomp)
843       ENDIF
844    ENDDO
845    !! 4.2 below: soil humidity for each soil layers
846
847    !! 4.2.1 integrate over the nslm levels
848    soilhum_decomp(:) = zero
849
850    DO l = 1, nslm !Loop over soil levels
851      DO m =1,nvm
852         soilhum_decomp(:) =soil_mc(:,l,pref_soil_veg(m))
853         IF (moist_func_Moyano) THEN
854            control_moist_soil(:,l,m) = control_moist_func_moyano (npts, soilhum_decomp,bulk_dens, &
855                 clay, carbon, veget_max)
856         ELSE
857            control_moist_soil(:,l,m) = control_moist_func (npts, soilhum_decomp)
858         ENDIF
859      ENDDO
860    ENDDO
861       
862        control_moist_soil(:,nslmd,:) = control_moist_soil(:,nslm,:)
863        WHERE (ISNAN(control_moist_soil(:,:,:)))
864            control_moist_soil(:,:,:)=un
865        ENDWHERE
866       
867    !! 5. fluxes from above ground litter to carbon pools and respiration
868    DO m = 2,nvm !Loop over PFTs
869
870      !! 5.1 structural litter: goes into active and slow carbon pools + respiration   
871
872      !! 5.1.1 total quantity of above ground structural litter which is decomposed
873
874      WHERE (soil_mc(:,1,pref_soil_veg(m)) .GT. zero)
875          fd(:) = dt/(litter_tau(istructural)*(un+poor_soils(:))) * &
876             control_temp_above(:,istructural) * control_moist_above(:,m) * exp( -litter_struct_coef * lignin_struc_above(:,m) )
877      ENDWHERE
878      DO k = 1,nelements
879         qd(:,k) = litter_above(:,istructural,m,k) * fd(:) * (un - flood_frac(:))
880         qd_flood(:,k) = litter_above(:,istructural,m,k) * fd(:) * flood_frac(:)/trois
881      ENDDO
882      IF (.NOT. lat_exp_doc) THEN
883         qd = qd + qd_flood
884         qd_flood = zero
885      ELSE
886         !Do nothing
887      ENDIF
888
889      !! 5.1.2 decompose same fraction of structural part of dead leaves. Not exact
890      !!       as lignine content is not the same as that of the total structural litter.
891      dead_leaves(:,m,istructural) = dead_leaves(:,m,istructural) * (un - fd(:)*(un - flood_frac(:)) - fd(:)*flood_frac(:)/trois)
892
893      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
894      ! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
895      ! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
896      ! time step and avoid some constructions that could create bug during future developments.
897      resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
898           ( 1. - frac_soil(istructural,iactive,iabove) ) * qd(:,icarbon) * &
899           ( 1. - lignin_struc_above(:,m) ) / dt
900      resp_hetero_flood(:,m) = resp_hetero_flood(:,m) + &
901           ( 1. - frac_soil(istructural,iactive,iabove) ) * qd_flood(:,icarbon) * &
902           ( 1. - lignin_struc_above(:,m) ) / dt
903
904      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
905      ! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
906      ! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
907      ! time step and avoid some constructions that could create bug during future developments.
908      resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
909           ( 1. - frac_soil(istructural,islow,iabove) ) * qd(:,icarbon) * lignin_struc_above(:,m) / dt
910      resp_hetero_flood(:,m) = resp_hetero_flood(:,m) + &
911           ( 1. - frac_soil(istructural,islow,iabove) ) * qd_flood(:,icarbon) * lignin_struc_above(:,m) / dt
912     
913      ! resp_hetero_litter_layer(:,m,1) = resp_hetero_litter_layer(:,m,1) + &
914      !      ( 1. - frac_soil(istructural,islow,iabove) ) * qd(:,icarbon) *lignin_struc_above(:,m) / dt
915        !  DO ig=1,npts
916        !       IF (qd(ig,k).LT.zero) THEN
917        !       WRITE(numout,*) 'STOMATElitt_ZHC1'
918        !       WRITE(numout,*) 'qd: ', qd(ig,k),MINVAL(soilcarbon_input),MINVAL(floodcarbon_input)
919        !       ENDIF
920        !       IF (fd(ig).LT.zero) THEN
921        !       WRITE(numout,*) 'STOMATElitt_ZHC1'
922        !       WRITE(numout,*) 'fd(ig): ', fd(ig),MINVAL(soilcarbon_input),MINVAL(floodcarbon_input)
923        !       ENDIF
924        !  ENDDO
925      !! 5.1.3 Calculation of the soilcarbon_input in gC m^-2 of water
926
927      DO k = 1, nelements
928         soilcarbon_input(:,m,1,istrabo,k) = (frac_soil(istructural,iactive,iabove) * ( 1. - lignin_struc_above(:,m) ) * qd(:,k) + &
929               frac_soil(istructural,islow,iabove) *  qd(:,k) * lignin_struc_above(:,m) ) /dt 
930         floodcarbon_input(:,m,iact,k) = frac_soil(istructural,iactive,iabove) * ( 1. - lignin_struc_above(:,m) ) * qd_flood(:,k) / dt
931         floodcarbon_input(:,m,islo,k) = frac_soil(istructural,islow,iabove) * qd_flood(:,k) * lignin_struc_above(:,m) /dt
932      ENDDO
933                                         
934      litter_above(:,istructural,m,:) = litter_above(:,istructural,m,:) - qd(:,:) - qd_flood(:,:)
935         
936      !! 5.2 above ground metabolic litter goes into active carbon pool + respiration
937         
938      !! 5.2.1 total quantity of aboveground metabolic litter that is decomposed
939      WHERE (soil_mc(:,1,pref_soil_veg(m)) .GT. zero)
940          fd(:) = dt/(litter_tau(imetabolic)*(un+poor_soils(:))) * control_temp_above(:,imetabolic) * control_moist_above(:,m)
941      ELSEWHERE
942          fd(:) = zero
943      ENDWHERE
944
945      DO k = 1,nelements         
946         qd(:,k) = litter_above(:,imetabolic,m,k) * fd(:) * (un - flood_frac(:))
947         qd_flood(:,k) = litter_above(:,imetabolic,m,k) * fd(:) * flood_frac(:)/trois
948      ENDDO
949          !
950      IF (.NOT. lat_exp_doc) THEN
951         qd = qd + qd_flood
952         qd_flood = zero
953      ELSE
954         !Do nothing
955      ENDIF
956         
957    !  DO ig=1,npts
958        !       IF (qd(ig,k).LT.zero) THEN
959        !       WRITE(numout,*) 'STOMATEcarb_ZHC2'
960        !       WRITE(numout,*) 'qd: ', qd(ig,k),MINVAL(soilcarbon_input),MINVAL(floodcarbon_input)
961        !       ENDIF
962        !       IF (fd(ig).LT.zero) THEN
963        !       WRITE(numout,*) 'STOMATEcarb_ZHC2'
964        !       WRITE(numout,*) 'fd(ig): ', fd(ig),MINVAL(soilcarbon_input),MINVAL(floodcarbon_input)
965        !       ENDIF
966        !  ENDDO
967      !! 5.2.2 decompose same fraction of metabolic part of dead leaves.
968      dead_leaves(:,m,imetabolic) = dead_leaves(:,m,imetabolic) * (un - fd(:)*(un - flood_frac(:)) - fd(:)*flood_frac(:)/trois)
969
970      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
971      ! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
972      ! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
973      ! time step and avoid some constructions that could create bug during future developments.
974      resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
975           ( 1. - frac_soil(imetabolic,iactive,iabove) ) * qd(:,icarbon) / dt
976      resp_hetero_flood(:,m) = resp_hetero_flood(:,m) + &
977           ( 1. - frac_soil(imetabolic,iactive,iabove) ) * qd_flood(:,icarbon) / dt
978
979      !! 5.2.3 Calculation of the soilcarbon_input in gC m^-2 of water
980      ! For above ground litter we assume that the soilcarbon_input coming from AB
981      ! litter decomposition is directly incorporated into the 1st soil
982      ! layers
983
984      DO k = 1, nelements
985         soilcarbon_input(:,m,1,imetabo,k) = frac_soil(imetabolic,iactive,iabove) *  qd(:,k) / dt             
986         floodcarbon_input(:,m,iact,k) = floodcarbon_input(:,m,iact,k) + frac_soil(imetabolic,iactive,iabove) *  qd_flood(:,k) / dt
987      ENDDO
988      litter_above(:,imetabolic,m,:) = litter_above(:,imetabolic,m,:) - qd(:,:) - qd_flood(:,:)
989
990    ENDDO
991
992!! 6. fluxes from below ground litter to carbon pools and respiration
993   DO l = 1,nslmd
994       DO m = 2,nvm !Loop over PFTs
995          !! 6.1 structural litter: goes into active and slow carbon pools respiration
996          !! 6.1.1 total quantity of below ground structural litter which is decomposed
997                 
998                  WHERE (soil_mc(:,l,pref_soil_veg(m)) .GT. zero)
999             fd(:) = dt/(litter_tau(istructural)*(un+poor_soils(:))) * &
1000                    control_temp_soil(:,l,istrbel) * control_moist_soil(:,l,m) * exp(-litter_struct_coef * lignin_struc_below(:,m,l) )
1001          ELSEWHERE
1002             fd(:) = zero
1003          ENDWHERE
1004                 
1005                  !WRITE(numout,*) 'SUMstomlitt2_control_temp_soil',SUM(control_temp_soil(:,l,istrbel))
1006                 
1007          DO k = 1,nelements
1008             qd(:,k) = litter_below(:,istructural,m,l,k) * fd(:) * (un - flood_frac(:))
1009             qd_flood(:,k) = litter_below(:,istructural,m,l,k) * fd(:) * flood_frac(:)/trois
1010          ENDDO
1011          !
1012          IF (.NOT. lat_exp_doc) THEN
1013             qd = qd + qd_flood
1014             qd_flood = zero
1015          ELSE
1016             !Do nothing
1017          ENDIF
1018
1019        !  DO ig=1,npts
1020                !       IF (qd(ig,k).LT.zero) THEN
1021                !               WRITE(numout,*) 'STOMATEcarb_ZHC3'
1022                !               WRITE(numout,*) 'qd: ', qd(ig,k),MINVAL(soilcarbon_input),MINVAL(floodcarbon_input)
1023                !       ENDIF
1024                !    IF (fd(ig).LT.zero) THEN
1025                !               WRITE(numout,*) 'STOMATEcarb_ZHC3'
1026                !               WRITE(numout,*) 'fd(ig): ', fd(ig),MINVAL(soilcarbon_input),MINVAL(floodcarbon_input)
1027                !    ENDIF
1028            !  ENDDO
1029      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
1030      ! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
1031      ! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
1032      ! time step and avoid some constructions that could create bug during future developments.
1033          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
1034               ( 1. - frac_soil(istructural,iactive,ibelow) ) * qd(:,icarbon) * &
1035               ( 1. - lignin_struc_below(:,m,l) ) / dt
1036          resp_hetero_flood(:,m) = resp_hetero_flood(:,m) + &
1037               ( 1. - frac_soil(istructural,iactive,ibelow) ) * qd_flood(:,icarbon) * &
1038               ( 1. - lignin_struc_below(:,m,l) ) / dt
1039
1040     
1041      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
1042      ! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
1043      ! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
1044      ! time step and avoid some constructions that could create bug during future developments.
1045          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
1046               ( 1. - frac_soil(istructural,islow,ibelow) ) * qd(:,icarbon) * lignin_struc_below(:,m,l) / dt
1047          resp_hetero_flood(:,m) = resp_hetero_flood(:,m) + &
1048               ( 1. - frac_soil(istructural,islow,ibelow) ) * qd_flood(:,icarbon) * lignin_struc_below(:,m,l) / dt
1049          !! 6.1 Calculation of the soilcarbon_input in gC m^-2
1050
1051          DO k = 1, nelements
1052             IF (l .GT. sro_bottom) THEN
1053                  soilcarbon_input(:,m,l,istrbel,k) = (frac_soil(istructural,iactive,ibelow) * ( 1. - lignin_struc_below(:,m,l) ) &
1054                       * (qd(:,k)+qd_flood(:,k)) + & 
1055                         (frac_soil(istructural,islow,ibelow) * (qd(:,k)+qd_flood(:,k)) * lignin_struc_below(:,m,l) )) / dt
1056             ELSE
1057                soilcarbon_input(:,m,l,istrbel,k) = (frac_soil(istructural,iactive,ibelow) * ( 1. - lignin_struc_below(:,m,l) ) &
1058                       * (qd(:,k)) + &
1059                         (frac_soil(istructural,islow,ibelow) * (qd(:,k)) * lignin_struc_below(:,m,l) )) / dt
1060                floodcarbon_input(:,m,iact,k) = floodcarbon_input(:,m,iact,k) + &
1061                     frac_soil(istructural,iactive,ibelow) * ( 1. - lignin_struc_below(:,m,l) ) * qd_flood(:,k) /dt
1062                floodcarbon_input(:,m,islo,k) = floodcarbon_input(:,m,islo,k) + (frac_soil(istructural,islow,ibelow) &
1063                     * qd_flood(:,k) * lignin_struc_below(:,m,l)) / dt
1064             ENDIF
1065          ENDDO
1066                 
1067          litter_below(:,istructural,m,l,:) = litter_below(:,istructural,m,l,:) - qd(:,:) - qd_flood(:,:)
1068                 
1069          !! 6.2 below ground metabolic litter goes into active carbon pool + respiration
1070
1071          !! 6.2.1 total quantity of belowground metabolic litter that is decomposed
1072          WHERE (soil_mc(:,l,pref_soil_veg(m)) .GT. zero)
1073             fd(:) = dt/(litter_tau(imetabolic)*(un+poor_soils(:))) * control_temp_soil(:,l,imetbel) * control_moist_soil(:,l,m)
1074          ELSEWHERE
1075             fd(:) = zero
1076          ENDWHERE
1077
1078          DO k = 1,nelements
1079             qd(:,k) = litter_below(:,imetabolic,m,l,k) * fd(:) * (un-flood_frac(:))
1080             qd_flood(:,k) = litter_below(:,imetabolic,m,l,k) * fd(:) * flood_frac(:)/trois
1081          ENDDO
1082          !
1083          IF (.NOT. lat_exp_doc) THEN
1084             qd = qd + qd_flood
1085             qd_flood = zero
1086          ELSE
1087             !Do nothing
1088          ENDIF
1089         
1090          DO ig=1,npts
1091                        IF (qd(ig,k).LT.zero) THEN
1092                                WRITE(numout,*) 'STOMATEcarb_ZHC4'
1093                                WRITE(numout,*) 'qd: ', qd(ig,k),MINVAL(soilcarbon_input),MINVAL(floodcarbon_input)
1094                        ENDIF
1095                    IF (fd(ig).LT.zero) THEN
1096                                WRITE(numout,*) 'STOMATEcarb_ZHC4'
1097                                WRITE(numout,*) 'fd(ig): ', fd(ig),MINVAL(soilcarbon_input),MINVAL(floodcarbon_input)
1098                    ENDIF
1099              ENDDO
1100      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
1101      ! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
1102      ! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
1103      ! time step and avoid some constructions that could create bug during future developments.
1104          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
1105               ( 1. - frac_soil(imetabolic,iactive,ibelow) ) * qd(:,icarbon) /dt
1106          resp_hetero_flood(:,m) = resp_hetero_flood(:,m) + &
1107               ( 1. - frac_soil(imetabolic,iactive,ibelow) ) * qd_flood(:,icarbon) /dt
1108         
1109          !! 6.2.3 Calculation of the soilcarbon_input in gC m^-2 of water
1110
1111          DO k = 1, nelements
1112             IF (l .GT. sro_bottom) THEN
1113                soilcarbon_input(:,m,l,imetbel,k) = frac_soil(imetabolic,iactive,ibelow) * (qd(:,k)+qd_flood(:,k))/ dt 
1114             ELSE
1115                soilcarbon_input(:,m,l,imetbel,k) = frac_soil(imetabolic,iactive,ibelow) * qd(:,k)/ dt
1116                floodcarbon_input(:,m,iact,k) = floodcarbon_input(:,m,iact,k) + frac_soil(imetabolic,iactive,ibelow) * qd_flood(:,k)/ dt
1117             ENDIF
1118          ENDDO
1119          !
1120          litter_below(:,imetabolic,m,l,:) = litter_below(:,imetabolic,m,l,:) - qd(:,:) - qd_flood(:,:)
1121          !
1122       ENDDO
1123   ENDDO
1124
1125
1126   DO m =2,nvm
1127
1128      litter_below_old(:,:,m,:,:)=zero
1129      litter_below_old_buffer(:,:,m,:,:)=zero
1130      litter_below_old(:,:,m,:,:)=litter_below(:,:,m,:,:)
1131      litter_above_old(:,:,m,:)=zero
1132      litter_above_old(:,:,m,:)=litter_above(:,:,m,:)
1133      DO k = 1, nelements ! Loop over elements
1134         DO j = 1,nlitt ! Loop over litter pools
1135
1136            litter_flux(:,j,m,1,k) = Dif_coef(:)*dt*ABS(litter_below_old(:,j,m,1,k)/(z_soil(1))-litter_below_old(:,j,m,2,k)/(z_soil(2)-z_soil(1)))/(z_soil(2)-z_soil(1))
1137            litter_above_flux(:,j,m,k) = Dif_coef(:)*dt*ABS(litter_above(:,j,m,k)/(z_litter*1e-3)-litter_below(:,j,m,1,k)/(z_litter*1e-3-z_soil(1)))
1138            ! BE CAREFUL: Here very important assumption we assume the above ground litter layers to be 10mm thick by default.
1139
1140            DO l= 2, nslm-1 ! Loop over soil layers
1141               litter_flux(:,j,m,l,k) = Dif_coef(:)*dt*ABS(litter_below_old(:,j,m,l,k)/(z_soil(l)-z_soil(l-1))-litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l)))/(z_soil(l+1)-z_soil(l))
1142            ENDDO
1143            litter_flux_old(:,j,m,:,k) = litter_flux(:,j,m,:,k)
1144            litter_above_flux_old(:,j,m,k) = litter_above_flux(:,j,m,k)
1145
1146            !Below we checked if in case that, in a given layer, you have diffusion in the above and below litters, both fluxes are not higher than the stocks of the given layer.
1147            WHERE (litter_above_old(:,j,m,k)/(z_litter*1e-3) .LT. litter_below_old(:,j,m,1,k)/(z_soil(1)) .AND. &
1148                 litter_below_old(:,j,m,2,k)/(z_soil(2)-z_soil(1)) .LT. litter_below_old(:,j,m,1,k)/(z_soil(1)) .AND. &
1149                 litter_above_flux_old(:,j,m,k)+litter_flux_old(:,j,m,1,k) .GT. litter_below_old(:,j,m,1,k))
1150               litter_above_flux(:,j,m,k) = litter_below_old(:,j,m,1,k)*(litter_above_flux_old(:,j,m,k)/(litter_above_flux_old(:,j,m,k)+litter_flux_old(:,j,m,1,k)))
1151               litter_flux(:,j,m,1,k) = litter_below_old(:,j,m,1,k)*(litter_flux_old(:,j,m,1,k)/(litter_above_flux_old(:,j,m,k)+litter_flux_old(:,j,m,1,k)))
1152            ELSEWHERE (litter_above_old(:,j,m,k)/(z_litter*1e-3) .GE. litter_below_old(:,j,m,1,k)/(z_soil(1)) .AND. &
1153                 litter_below_old(:,j,m,2,k)/(z_soil(2)-z_soil(1)) .LT. litter_below_old(:,j,m,1,k)/(z_soil(1)) .AND. &
1154                 litter_below_old(:,j,m,1,k) + litter_above_flux_old(:,j,m,k) - litter_flux_old(:,j,m,1,k) .LE. min_stomate)
1155               litter_above_flux(:,j,m,k) = litter_above_flux_old(:,j,m,k)
1156               litter_flux(:,j,m,1,k) = litter_below_old(:,j,m,1,k) + litter_above_flux_old(:,j,m,k)
1157            ELSEWHERE (litter_above_old(:,j,m,k)/(z_litter*1e-3) .LT. litter_below_old(:,j,m,1,k)/(z_soil(1)) .AND. &
1158                 litter_below_old(:,j,m,2,k)/(z_soil(2)-z_soil(1)) .GE. litter_below_old(:,j,m,1,k)/(z_soil(1)) .AND. &
1159                 litter_below_old(:,j,m,1,k) - litter_above_flux_old(:,j,m,k) + litter_flux_old(:,j,m,1,k) .LE. min_stomate)
1160               litter_above_flux(:,j,m,k) = litter_below_old(:,j,m,1,k) + litter_flux_old(:,j,m,1,k) 
1161               litter_flux(:,j,m,1,k) = litter_flux_old(:,j,m,1,k)
1162            ELSEWHERE
1163               litter_above_flux(:,j,m,k) = litter_above_flux_old(:,j,m,k)
1164               litter_flux(:,j,m,1,k) = litter_flux_old(:,j,m,1,k)
1165            ENDWHERE
1166
1167
1168            DO l =1, nslm-2
1169               WHERE (litter_below_old(:,j,m,l,k)/(z_soil(l)-z_soil(l-1)) .LT. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l)) .AND. &
1170                    litter_below_old(:,j,m,l+2,k)/(z_soil(l+2)-z_soil(l+1)) .LT. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l)) .AND. &
1171                    litter_flux_old(:,j,m,l,k)+litter_flux_old(:,j,m,l+1,k) .GT. litter_below_old(:,j,m,l+1,k))
1172                  litter_flux(:,j,m,l,k) = litter_below_old(:,j,m,l+1,k)*(litter_flux_old(:,j,m,l,k)/(litter_flux_old(:,j,m,l,k)+litter_flux_old(:,j,m,l+1,k)))
1173                  litter_flux(:,j,m,l+1,k) = litter_below_old(:,j,m,l+1,k)*(litter_flux_old(:,j,m,l+1,k)/(litter_flux_old(:,j,m,l,k)+litter_flux_old(:,j,m,l+1,k)))
1174               ELSEWHERE (litter_below_old(:,j,m,l,k)/(z_soil(l)-z_soil(l-1)) .GE. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l)) .AND. &
1175                    litter_below_old(:,j,m,l+2,k)/(z_soil(l+2)-z_soil(l+1)) .LT. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l)) .AND. &
1176                    litter_below_old(:,j,m,l+1,k) + litter_flux_old(:,j,m,l,k) - litter_flux_old(:,j,m,l+1,k) .LE. min_stomate)
1177                  litter_flux(:,j,m,l,k) = litter_flux_old(:,j,m,l,k)
1178                  litter_flux(:,j,m,l+1,k) = litter_below_old(:,j,m,l+1,k) + litter_flux_old(:,j,m,l,k)
1179               ELSEWHERE (litter_below_old(:,j,m,l,k)/(z_soil(l)-z_soil(l-1)) .LT. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l)) .AND. &
1180                    litter_below_old(:,j,m,l+2,k)/(z_soil(l+2)-z_soil(l+1)) .GE. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l)) .AND. &
1181                    litter_below_old(:,j,m,l+1,k) - litter_flux_old(:,j,m,l,k) + litter_flux_old(:,j,m,l+1,k) .LE. min_stomate)
1182                  litter_flux(:,j,m,l,k) = litter_below_old(:,j,m,l+1,k) + litter_flux_old(:,j,m,l+1,k)
1183                  litter_flux(:,j,m,l+1,k) = litter_flux_old(:,j,m,l+1,k)
1184               ELSEWHERE
1185                  litter_flux(:,j,m,l,k) = litter_flux_old(:,j,m,l,k)
1186                  litter_flux(:,j,m,l+1,k) = litter_flux_old(:,j,m,l+1,k)
1187               ENDWHERE
1188            ENDDO
1189
1190            WHERE ((litter_above_old(:,j,m,k)/(z_litter*1e-3) .LT. litter_below_old(:,j,m,1,k)/(z_soil(1))) .AND. &
1191                 ((litter_above_flux(:,j,m,k) - litter_below_old(:,j,m,1,k)) .GE. min_stomate))
1192               litter_above(:,j,m,k) = litter_above_old(:,j,m,k) + litter_below_old(:,j,m,1,k) 
1193               litter_below_old_buffer(:,j,m,1,k) = -litter_below_old(:,j,m,1,k) 
1194               litter_below(:,j,m,1,k) = zero
1195            ELSEWHERE ((litter_above_old(:,j,m,k)/(z_litter*1e-3) .LT. litter_below_old(:,j,m,1,k)/(z_soil(1))) .AND. &
1196                 ((litter_below_old(:,j,m,1,k) - litter_above_flux(:,j,m,k)) .GT. min_stomate))
1197               litter_above(:,j,m,k) = litter_above_old(:,j,m,k) + litter_above_flux(:,j,m,k) 
1198               litter_below_old_buffer(:,j,m,1,k) =  - litter_above_flux(:,j,m,k)
1199               litter_below(:,j,m,1,k) = litter_below_old(:,j,m,1,k) - litter_above_flux(:,j,m,k)
1200            ELSEWHERE ((litter_above_old(:,j,m,k)/(z_litter*1e-3) .GT. litter_below_old(:,j,m,1,k)/(z_soil(1))) .AND. &
1201                 ((litter_above_flux(:,j,m,k) - litter_above_old(:,j,m,k)) .GE. min_stomate))
1202               litter_above(:,j,m,k) = zero 
1203               litter_below_old_buffer(:,j,m,1,k) = litter_above_old(:,j,m,k)
1204               litter_below(:,j,m,1,k) = litter_below_old(:,j,m,1,k) + litter_above_old(:,j,m,k)
1205            ELSEWHERE ((litter_above_old(:,j,m,k)/(z_litter*1e-3) .GT. litter_below_old(:,j,m,1,k)/(z_soil(1))) .AND. &
1206                 ((litter_above_old(:,j,m,k) - litter_above_flux(:,j,m,k)) .GT. min_stomate))
1207               litter_above(:,j,m,k) = litter_above_old(:,j,m,k) - litter_above_flux(:,j,m,k)
1208               litter_below_old_buffer(:,j,m,1,k) = litter_above_flux(:,j,m,k)
1209               litter_below(:,j,m,1,k) = litter_below_old(:,j,m,1,k) + litter_above_flux(:,j,m,k)
1210            ELSEWHERE
1211               litter_above(:,j,m,k) = litter_above_old(:,j,m,k)
1212               litter_below(:,j,m,1,k) = litter_below_old(:,j,m,1,k)
1213            ENDWHERE
1214                 
1215            DO l= 1, nslm-1 ! Loop over soil layers
1216               WHERE ((litter_below_old(:,j,m,l,k)/(z_soil(l)-z_soil(l-1)) .LT. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l))) .AND. &
1217                    ((litter_flux(:,j,m,l,k) - litter_below_old(:,j,m,l+1,k)) .GE. min_stomate))
1218                  litter_below(:,j,m,l,k) = litter_below_old(:,j,m,l,k) + litter_below_old(:,j,m,l+1,k) + litter_below_old_buffer(:,j,m,l,k)
1219                  litter_below_old_buffer(:,j,m,l+1,k) = -litter_below_old(:,j,m,l+1,k)
1220                  litter_below(:,j,m,l+1,k) = zero
1221               ELSEWHERE ((litter_below_old(:,j,m,l,k)/(z_soil(l)-z_soil(l-1)) .LT. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l))) .AND. &
1222                    ((litter_below_old(:,j,m,l+1,k) - litter_flux(:,j,m,l,k)) .GT. min_stomate))
1223                  litter_below(:,j,m,l,k) = litter_below_old(:,j,m,l,k) + litter_flux(:,j,m,l,k) + litter_below_old_buffer(:,j,m,l,k)
1224                  litter_below_old_buffer(:,j,m,l+1,k) =  - litter_flux(:,j,m,l,k)
1225                  litter_below(:,j,m,l+1,k) = litter_below_old(:,j,m,l+1,k) - litter_flux(:,j,m,l,k)
1226               ELSEWHERE ((litter_below_old(:,j,m,l,k)/(z_soil(l)-z_soil(l-1)) .GT. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l))) .AND. &
1227                    ((litter_flux(:,j,m,l,k) - litter_below_old(:,j,m,l,k)) .GE. min_stomate))
1228                  litter_below(:,j,m,l,k) = zero + litter_below_old_buffer(:,j,m,l,k)
1229                  litter_below_old_buffer(:,j,m,l+1,k) = litter_below_old(:,j,m,l,k)
1230                  litter_below(:,j,m,l+1,k) = litter_below_old(:,j,m,l+1,k) + litter_below_old(:,j,m,l,k)
1231               ELSEWHERE ((litter_below_old(:,j,m,l,k)/(z_soil(l)-z_soil(l-1)) .GT. litter_below_old(:,j,m,l+1,k)/(z_soil(l+1)-z_soil(l))) .AND. &
1232                    ((litter_below_old(:,j,m,l,k) - litter_flux(:,j,m,l,k)) .GT. min_stomate))
1233                  litter_below(:,j,m,l,k) = litter_below_old(:,j,m,l,k) - litter_flux(:,j,m,l,k) + litter_below_old_buffer(:,j,m,l,k)
1234                  litter_below_old_buffer(:,j,m,l+1,k) = litter_flux(:,j,m,l,k)
1235                  litter_below(:,j,m,l+1,k) = litter_below_old(:,j,m,l+1,k) + litter_flux(:,j,m,l,k)
1236               ELSEWHERE
1237                  litter_below(:,j,m,l,k) = litter_below_old(:,j,m,l,k)
1238                  litter_below(:,j,m,l+1,k) = litter_below_old(:,j,m,l+1,k)
1239               ENDWHERE
1240            ENDDO ! End loop over soil layers
1241         ENDDO ! End loop over litter pools
1242      ENDDO ! End loop over elements
1243   ENDDO ! End loop over PFT
1244
1245  !! 8. calculate fraction of total soil covered by dead leaves
1246
1247   CALL deadleaf (npts, veget_max, dead_leaves, deadleaf_cover)
1248
1249   !! 9. Check mass balance closure
1250
1251   !! 9.1 Calculate components of the mass balance
1252   pool_end_after(:,:,:) = zero
1253   pool_end(:,:,:) = zero
1254   IF (ld_doc) THEN
1255      DO m = 2, nvm
1256         DO i = 1,nlitt
1257            DO l = 1, nslmd
1258               pool_end(:,m,icarbon) = pool_end(:,m,icarbon) + &
1259                    ( litter_below(:,i,m,l,icarbon) ) * veget_max(:,m)
1260            ENDDO
1261         ENDDO
1262      ENDDO
1263
1264      DO m = 2, nvm
1265         DO i = 1,nlitt
1266            pool_end(:,m,icarbon) = pool_end(:,m,icarbon) + &
1267                 ( litter_above(:,i,m,icarbon) ) * veget_max(:,m)
1268         ENDDO
1269      ENDDO
1270
1271
1272    !! 9.2 Calculate mass balance
1273    !  Note that soilcarbon is transfered to other pools but that the pool
1274    !  was not updated. We should not account for it in ::pool_end
1275
1276      check_intern(:,:,:,:) = zero
1277      check_intern(:,:,iatm2land,icarbon) = zero 
1278
1279   
1280      check_intern(:,:,iland2atm,icarbon) = -un * (resp_hetero_litter(:,:)+resp_hetero_flood(:,:)) * &
1281           veget_max(:,:) * dt
1282   
1283
1284      DO m = 1,nvm
1285         DO l = 1, nslmd
1286            check_intern(:,m,ilat2out,icarbon) = check_intern(:,m,ilat2out,icarbon) -un * &
1287                 (soilcarbon_input(:,m,l,imetbel,icarbon)+soilcarbon_input(:,m,l,istrbel,icarbon)) * &
1288                 veget_max(:,m) * dt
1289         ENDDO
1290      ENDDO
1291      DO m = 1,nvm
1292         check_intern(:,m,ilat2out,icarbon) = check_intern(:,m,ilat2out,icarbon) -un * &
1293              (soilcarbon_input(:,m,1,imetabo,icarbon)+soilcarbon_input(:,m,1,istrabo,icarbon)    &
1294              +floodcarbon_input(:,m,iact,icarbon)+floodcarbon_input(:,m,islo,icarbon)) * &
1295              veget_max(:,m) * dt
1296      ENDDO
1297
1298      check_intern(:,:,ilat2in,icarbon) = zero
1299   
1300      check_intern(:,:,ipoolchange,icarbon) = -un * (pool_end(:,:,icarbon) - &
1301           pool_start(:,:,icarbon))
1302   
1303      closure_intern = zero
1304      DO i = 1,nmbcomp
1305         closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + &
1306              check_intern(:,:,i,icarbon)
1307      ENDDO
1308
1309    !! 9.3 Write verdict
1310
1311      IF (SUM(SUM(closure_intern(:,:,icarbon),2),1) .LT. min_stomate .AND. &
1312           SUM(SUM(closure_intern(:,:,icarbon),2),1) .GT. -min_stomate) THEN
1313         WRITE(numout,*) 'Mass balance closure in stomate_litter.f90'
1314      ELSE
1315         WRITE(numout,*) 'Error: mass balance is not closed in stomate_litter.f90'
1316         WRITE(numout,*) '   Difference is, ', SUM(SUM(closure_intern(:,:,icarbon),2),1)
1317         WRITE(numout,*) '   Difference is, ', closure_intern(:,:,icarbon)
1318      ENDIF
1319   ENDIF
1320
1321
1322   IF (printlev>=4) WRITE(numout,*) 'Leaving littercalc'
1323
1324 END SUBROUTINE littercalc
1325
1326
1327!! ==============================================================================================================================\n
1328!! SUBROUTINE   : deadleaf
1329!!
1330!>\BRIEF        This routine calculates the deadleafcover.
1331!!
1332!! DESCRIPTION  : It first calculates the lai corresponding to the dead leaves (LAI) using
1333!! the dead leaves carbon content (DL) the specific leaf area (sla) and the
1334!! maximal coverage fraction of a PFT (vegetmax) using the following equations:
1335!! \latexonly
1336!! \input{deadleaf1.tex}
1337!! \endlatexonly
1338!! \n
1339!! Then, the dead leaf cover (DLC) is calculated as following:\n
1340!! \latexonly
1341!! \input{deadleaf2.tex}
1342!! \endlatexonly
1343!! \n
1344!!
1345!! RECENT CHANGE(S) : None
1346!!
1347!! MAIN OUTPUT VARIABLE: ::deadleaf_cover
1348!!
1349!! REFERENCE(S) : None
1350!!
1351!! FLOWCHART : None
1352!! \n
1353!_ ================================================================================================================================
1354
1355  SUBROUTINE deadleaf (npts, veget_max, dead_leaves, deadleaf_cover)
1356
1357  !! 0. Variable and parameter declaration
1358   
1359    !! 0.1 Input variables
1360
1361    INTEGER(i_std), INTENT(in)                          :: npts           !! Domain size - number of grid pixels (unitless)
1362    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in)  :: dead_leaves    !! Dead leaves per ground unit area, per PFT,
1363                                                                          !! metabolic and structural 
1364                                                                          !! @tex $(gC m^{-2})$ @endtex
1365    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)          :: veget_max      !! PFT "Maximal" coverage fraction of a PFT defined in
1366                                                                          !! the input vegetation map
1367                                                                          !! @tex $(m^2 m^{-2})$ @endtex
1368   
1369    !! 0.2 Output variables
1370   
1371    REAL(r_std), DIMENSION(npts), INTENT(out)           :: deadleaf_cover !! Fraction of soil covered by dead leaves over all PFTs
1372                                                                          !! (0-1, unitless)
1373
1374    !! 0.3 Modified variables
1375
1376    !! 0.4 Local variables
1377
1378    REAL(r_std), DIMENSION(npts)                        :: dead_lai       !! LAI of dead leaves @tex $(m^2 m^{-2})$ @endtex
1379    INTEGER(i_std)                                      :: j              !! Index (unitless)
1380!_ ================================================================================================================================
1381   
1382  !! 1. LAI of dead leaves
1383 
1384    dead_lai(:) = zero
1385
1386    DO j = 2,nvm !Loop over PFTs
1387       dead_lai(:) = dead_lai(:) + ( dead_leaves(:,j,imetabolic) + dead_leaves(:,j,istructural) ) * sla(j) &
1388            * veget_max(:,j)
1389    ENDDO
1390
1391  !! 2. fraction of soil covered by dead leaves
1392
1393    deadleaf_cover(:) = un - exp( - 0.5 * dead_lai(:) )
1394
1395    IF (printlev>=4) WRITE(numout,*) 'Leaving deadleaf'
1396
1397  END SUBROUTINE deadleaf
1398
1399
1400!! ================================================================================================================================
1401!! FUNCTION     : control_moist_func_moyano
1402!!
1403!>\BRIEF        Calculate moisture control for litter and soild C decomposition
1404!!
1405!! DESCRIPTION  : Calculate moisture control factor applied
1406!! to litter decomposition and to soil carbon decomposition in
1407!! stomate_soilcarbon.f90 using the following equation: \n
1408!! \latexonly
1409!! \input{control_moist_func1.tex}
1410!! \endlatexonly
1411!! \n
1412!! with M the moisture control factor and soilmoisutre, the soil moisture
1413!! calculated in sechiba.
1414!! Then, the function is ranged between 0.25 and 1:\n
1415!! \latexonly
1416!! \input{control_moist_func2.tex}
1417!! \endlatexonly
1418!! \n
1419!! RECENT CHANGE(S) : None
1420!!
1421!! RETURN VALUE : ::moistfunc_result
1422!!
1423!! REFERENCE(S) : None
1424!!
1425!! FLOWCHART : None
1426!! \n
1427!_ ================================================================================================================================
1428 
1429  FUNCTION control_moist_func_moyano (npts, moist_in,bulk_dens, clay, carbon,veget_max) RESULT (moistfunc_result)
1430
1431  !! 0. Variable and parameter declaration
1432
1433    !! 0.1 Input variables
1434
1435    INTEGER(i_std), INTENT(in)               :: npts                !! Domain size - number of grid pixel (unitless)
1436    REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in            !! relative humidity (unitless)
1437    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: clay               !! Clay fraction (unitless, 0-1)
1438    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: bulk_dens          !! Variable local of bulk density for the moment must change in the futur (kg m-3)
1439    REAL(r_std), DIMENSION(npts,ncarb,nvm,nslmd), INTENT(in)    :: carbon             !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f
1440    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                  :: veget_max          !! Maximum fraction of vegetation type including
1441                                                                                      !! non-biological fraction (unitless)
1442
1443    !! 0.2 Output variables
1444
1445    REAL(r_std), DIMENSION(npts)          :: moistfunc_result    !! Moisture control factor (0-1, unitless)
1446
1447    !! 0.3 Modified variables
1448
1449    !! 0.4 Local variables
1450   REAL(r_std), DIMENSION(npts)           :: total_soc              !! Total soil organic carbon for a grid pixel (g-C g-1 soil)
1451   INTEGER(i_std)                         :: k,i,j                  !! Indices
1452   INTEGER, PARAMETER                     :: nummoist = 120         !! We used numoist then divived by 100. because to calculate the function
1453                                                                    !! we need to predict PRsr for all the possible soil moisture values
1454   REAL(r_std), DIMENSION(npts,nummoist)  :: PRsr                   !! Proportional response or soil respiration at moist_in
1455   REAL(r_std), DIMENSION(npts,nummoist)  :: SR                     !! Soil respiration at moist_in to facilitate the understanding regarding
1456                                                                    !! the Moyano et al paper we keep the same nomenclature than in the paper
1457                                                                    !! but **BE CAREFUL** it does not correspond to the variable resp_hetero_soil
1458                                                                    !! or resp_hetero_litter in ORCHIDEE
1459   REAL(r_std), DIMENSION(0:nslm)   :: z_soil                 !! Soil levels (m)
1460!$OMP THREADPRIVATE(z_soil)
1461   REAL(r_std), DIMENSION(npts,nummoist)  :: moistfunc              !! An intermediate variable to calculate the Moisture control factor
1462   INTEGER(i_std)                         :: ind                    !! A local variable to know what is the index of the max value for the statistical function
1463   INTEGER(i_std), DIMENSION(npts)        :: ind_i                  !! A local variable to know what is the index of the max value for the
1464                                                                    !! statistical function at each points
1465   REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: moistfunc_below_1 !! A allocatable variable to store the value of moistfunc that are below the max value
1466                                                                       !! to then rescale at is it done in Moyano et al. 2012 BG
1467!$OMP THREADPRIVATE(moistfunc_below_1)
1468   INTEGER(i_std)                         :: ier                    !! Check errors in netcdf call (unitless)
1469   LOGICAL                                :: l_error                !! Check errors in netcdf call
1470
1471!_ ================================================================================================================================
1472
1473   !! 1.3 soil levels
1474   z_soil(0) = zero
1475   z_soil(1:nslm) = zlt(1:nslm)
1476
1477   ! In Moyano et al. 2012 BG, the function obtained is based on soil
1478   ! incubations were litter is generally manually removed, therefore we used
1479   ! only the variable carbon corresponding to the soil organic carbon even
1480   ! if some litter (metabolic litter) is not easy to remove manually since it
1481   ! corresponds to roots exudates. Therefore, we probably underestimate a bit
1482   ! the soil organic carbon compared to what is used in Moyano et al., 2012
1483   ! BG.
1484
1485   PRsr(:,:) = zero
1486   SR(:,:) = zero
1487   moistfunc(:,:) = zero
1488   l_error = .FALSE.
1489
1490   total_soc(:) = zero
1491   DO k = 1,nvm
1492      DO j = 1,nslm
1493         total_soc(:) =  total_soc(:) + &
1494              ((carbon(:,iactive,k,j)*veget_max(:,k) + &
1495              carbon(:,islow,k,j)*veget_max(:,k) + &
1496              carbon(:,ipassive,k,j)*veget_max(:,k))* &
1497              ((z_soil(j)-z_soil(j-1))/z_soil(nslm))) / &
1498              (bulk_dens(:) ) 
1499      ENDDO
1500   ENDDO
1501
1502   DO k = 1,nummoist
1503
1504      PRsr(:,k) = 1.11066-0.83344*(k/100.) + 1.48095*((k/100.)**2) - 1.02959*((k/100.)**3) + 0.07995*clay(:) + 1.27892*total_soc(:)
1505     
1506   ENDDO
1507
1508   SR(:,1)=PRsr(:,1)
1509   DO k=2,nummoist
1510      SR(:,k)= SR(:,k-1)* PRsr(:,k)
1511   ENDDO
1512
1513   DO i = 1,npts
1514      DO k= 1, nummoist
1515         moistfunc(i,k) = SR(i,k)/MAXVAL(SR(i,:))
1516      ENDDO
1517   ENDDO
1518
1519   ind_i(:) = MAXLOC(moistfunc(:,:),DIM=2)
1520   
1521   DO i = 1,npts
1522      ind=ind_i(i)
1523
1524      ALLOCATE(moistfunc_below_1(npts,ind),stat=ier)
1525      l_error = l_error .OR. (ier /= 0)
1526      IF (l_error) THEN
1527         WRITE(numout,*) 'Memory allocation error for moistfunc_below_1. '
1528      ENDIF
1529
1530      DO k= 1,ind
1531         moistfunc_below_1(i,k)= moistfunc(i,k)
1532      ENDDO
1533
1534      DO k= 1,ind           
1535         moistfunc(i,k)=moistfunc(i,k) - MINVAL(moistfunc_below_1(i,:))
1536      ENDDO
1537
1538      DO k= 1,ind
1539         moistfunc_below_1(i,k)= moistfunc(i,k)
1540      ENDDO
1541
1542      DO k=1,ind
1543         moistfunc(i,k)=moistfunc(i,k)/MAXVAL(moistfunc_below_1(i,:))
1544      ENDDO
1545
1546      IF (ALLOCATED(moistfunc_below_1)) DEALLOCATE(moistfunc_below_1)
1547
1548   ENDDO
1549   DO i=1,npts
1550      IF (NINT(moist_in(i)*100.) .GT. zero) THEN
1551         moistfunc_result(i) = moistfunc(i,NINT(moist_in(i)*100.))
1552      ELSE
1553         moistfunc_result(i) = zero
1554      ENDIF
1555   ENDDO
1556
1557 END FUNCTION control_moist_func_moyano
1558
1559!! ================================================================================================================================
1560!! FUNCTION     : control_moist_func
1561!!
1562!>\BRIEF        Calculate moisture control for litter and soild C decomposition
1563!!
1564!! DESCRIPTION  : Calculate moisture control factor applied
1565!! to litter decomposition and to soil carbon decomposition in
1566!! stomate_soilcarbon.f90 using the following equation: \n
1567!! \latexonly
1568!! \input{control_moist_func1.tex}
1569!! \endlatexonly
1570!! \n
1571!! with M the moisture control factor and soilmoisutre, the soil moisture
1572!! calculated in sechiba.
1573!! Then, the function is ranged between 0.25 and 1:\n
1574!! \latexonly
1575!! \input{control_moist_func2.tex}
1576!! \endlatexonly
1577!! \n
1578!! RECENT CHANGE(S) : None
1579!!
1580!! RETURN VALUE : ::moistfunc_result
1581!!
1582!! REFERENCE(S) : None
1583!!
1584!! FLOWCHART : None
1585!! \n
1586!_ ================================================================================================================================
1587
1588  FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)
1589
1590  !! 0. Variable and parameter declaration
1591
1592    !! 0.1 Input variables
1593
1594    INTEGER(i_std), INTENT(in)               :: npts                !! Domain size - number of grid pixel (unitless)
1595    REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in            !! relative humidity (unitless)
1596
1597
1598    !! 0.2 Output variables
1599
1600    REAL(r_std), DIMENSION(npts)          :: moistfunc_result    !! Moisture control factor (0-1, unitless)
1601
1602    !! 0.3 Modified variables
1603
1604
1605!_ ================================================================================================================================
1606
1607
1608          moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
1609          moistfunc_result(:) = MAX( 0.25_r_std, MIN( un, moistfunc_result(:) ) )
1610   
1611  END FUNCTION control_moist_func
1612
1613!! ================================================================================================================================
1614!! FUNCTION     : control_temp_func
1615!!
1616!>\BRIEF        Calculate temperature control for litter and soild C decomposition
1617!!
1618!! DESCRIPTION  : Calculate temperature control factor applied
1619!! to litter decomposition and to soil carbon decomposition in
1620!! stomate_soilcarbon.f90 using the following equation: \n
1621!! \latexonly
1622!! \input{control_temp_func1.tex}
1623!! \endlatexonly
1624!! \n
1625!! with T the temperature control factor, temp the temperature in Kelvin of
1626!! the air (for aboveground litter) or of the soil (for belowground litter
1627!! and soil)
1628!! Then, the function is limited in its maximal range to 1:\n
1629!! \latexonly
1630!! \input{control_temp_func2.tex}
1631!! \endlatexonly
1632!! \n
1633!! RECENT CHANGE(S) : None
1634!!
1635!! RETURN VALUE: ::tempfunc_result
1636!!
1637!! REFERENCE(S) : None
1638!!
1639!! FLOWCHART : None
1640!! \n
1641!_ ================================================================================================================================
1642
1643  FUNCTION control_temp_func (npts, temp_in,tau) RESULT (tempfunc_result)
1644
1645  !! 0. Variable and parameter declaration
1646   
1647    !! 0.1 Input variables
1648    INTEGER(i_std), INTENT(in)                 :: npts            !! Domain size - number of land pixels (unitless)
1649    REAL(r_std), DIMENSION(npts), INTENT(in)   :: temp_in         !! Temperature (K)
1650    REAL(r_std), DIMENSION(npts), INTENT(in)   :: tau             !! residence time of the pool considered (days)
1651    !! 0.2 Output variables
1652    REAL(r_std), DIMENSION(npts)               :: tempfunc_result !! Temperature control factor (0-1, unitless)
1653
1654    !! 0.3 Modified variables
1655
1656    !! 0.4 Local variables
1657    REAL(r_std), DIMENSION(npts)               :: Ea              !! Activation Energy
1658    REAL(r_std), DIMENSION(npts)               :: Q10_soil        !! Q10 values calculated based on Ea
1659!_ ================================================================================================================================
1660
1661    ! This relationship is based on a reanalysis of the data published by LefÚvre et al. 2014 in global change biology
1662    ! We use MAX(cent,(temp_in(:)-ZeroCelsius)*(temp_in(:)-ZeroCelsius)))) because for values between -10 and 10 degC
1663    ! the Q10 values are close to infinity
1664    !    Ea(:)=a_term_Q10_soil*tau + b_term_Q10_soil
1665    !    Q10_soil(:)=exp((dix*Ea(:))/(RR*MAX(cent,(temp_in(:)-ZeroCelsius)*(temp_in(:)-ZeroCelsius))))
1666    !    tempfunc_result(:) = exp( log(Q10_soil(:)) * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
1667       
1668        tempfunc_result(:) = exp( soil_Q10 * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
1669    tempfunc_result(:) = MIN( un, tempfunc_result(:) )
1670    WHERE (temp_in(:) .LT. zero) 
1671        tempfunc_result(:) = zero
1672    ENDWHERE
1673  END FUNCTION control_temp_func
1674
1675END MODULE stomate_litter
Note: See TracBrowser for help on using the repository browser.