source: branches/publications/ORCHIDEE_SOM_13C_r4859/src_stomate/stomate_litter.f90 @ 7346

Last change on this file since 7346 was 4858, checked in by marta.camino, 7 years ago

13C fractionation and root enrichment included and depth dependent diffusion coef

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