source: branches/publications/ORCHIDEE-LEAK-r5919/src_stomate/stomate_litter.f90

Last change on this file was 5919, checked in by ronny.lauerwald, 5 years ago

ORCHILEAK, version used for trends and biases in NEE and NBP in the Amazon basin

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