source: branches/publications/ORCHIDEE-LEAK-r5919/src_stomate/stomate_soilcarbon.f90 @ 5925

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

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

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 86.2 KB
Line 
1
2! MODULE       : stomate_soilcarbon
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       Calculate soil dynamics largely following the Century model
10!
11!!     
12!!\n DESCRIPTION: 11 layers discretization: Layers 1-11 are soil layers
13!!
14!! RECENT CHANGE(S): None
15!!
16!! SVN          :
17!! $HeadURL$
18!! $Date$
19!! $Revision$
20!! \n
21!_ ================================================================================================================================
22
23MODULE stomate_soilcarbon
24
25  ! modules used:
26  USE xios_orchidee
27  USE ioipsl_para
28  USE stomate_data
29  USE constantes
30  USE constantes_soil
31  USE pft_parameters
32  USE pft_parameters_var
33  USE vertical_soil
34  IMPLICIT NONE
35
36  ! private & public routines
37
38  PRIVATE
39  PUBLIC soilcarbon,soilcarbon_clear
40
41  ! Variables shared by all subroutines in this module
42 
43  LOGICAL, SAVE    :: firstcall = .TRUE.   !! Is this the first call? (true/false)
44!$OMP THREADPRIVATE(firstcall)
45
46CONTAINS
47
48
49!! ================================================================================================================================
50!!  SUBROUTINE   : soilcarbon_clear
51!!
52!>\BRIEF        Set the flag ::firstcall to .TRUE. and as such activate sections 1.1.2 and 1.2 of the subroutine soilcarbon
53!! (see below).
54!!
55!_ ================================================================================================================================
56 
57  SUBROUTINE soilcarbon_clear
58    firstcall=.TRUE.
59  ENDSUBROUTINE soilcarbon_clear
60
61
62!! ================================================================================================================================
63!!  SUBROUTINE   : soilcarbon
64!!
65!>\BRIEF        Computes the soil respiration and carbon stocks, essentially
66!! following Parton et al. (1987).
67!!
68!! DESCRIPTION  : The soil is divided into 3 carbon pools, with different
69!! characteristic turnover times : active (1-5 years), slow (20-40 years)
70!! and passive (200-1500 years).\n
71!! There are three types of carbon transferred in the soil:\n
72!! - carbon input in active and slow pools from litter decomposition,\n
73!! - carbon fluxes between the three pools,\n
74!! - carbon losses from the pools to the atmosphere, i.e., soil respiration.\n
75!!
76!! The subroutine performs the following tasks:\n
77!!
78!! Section 1.\n
79!! The flux fractions (f) between carbon pools are defined based on Parton et
80!! al. (1987). The fractions are constants, except for the flux fraction from
81!! the active pool to the slow pool, which depends on the clay content,\n
82!! \latexonly
83!! \input{soilcarbon_eq1.tex}
84!! \endlatexonly\n
85!! In addition, to each pool is assigned a constant turnover time.\n
86!!
87!! Section 2.\n
88!! The carbon input, calculated in the stomate_litter module, is added to the
89!! carbon stock of the different pools.\n
90!!
91!! Section 3.\n
92!! First, the outgoing carbon flux of each pool is calculated. It is
93!! proportional to the product of the carbon stock and the ratio between the
94!! iteration time step and the residence time:\n
95!! \latexonly
96!! \input{soilcarbon_eq2.tex}
97!! \endlatexonly
98!! ,\n
99!! Note that in the case of crops, the additional multiplicative factor
100!! integrates the faster decomposition due to tillage (following Gervois et
101!! al. (2008)).
102!! In addition, the flux from the active pool depends on the clay content:\n
103!! \latexonly
104!! \input{soilcarbon_eq3.tex}
105!! \endlatexonly
106!! ,\n
107!! Each pool is then cut from the carbon amount corresponding to each outgoing
108!! flux:\n
109!! \latexonly
110!! \input{soilcarbon_eq4.tex}
111!! \endlatexonly\n
112!! Second, the flux fractions lost to the atmosphere is calculated in each pool
113!! by subtracting from 1 the pool-to-pool flux fractions. The soil respiration
114!! is then the summed contribution of all the pools,\n
115!! \latexonly
116!! \input{soilcarbon_eq5.tex}
117!! \endlatexonly\n
118!! Finally, each carbon pool accumulates the contribution of the other pools:
119!! \latexonly
120!! \input{soilcarbon_eq6.tex}
121!! \endlatexonly
122!!
123!! RECENT CHANGE(S): None
124!!
125!! MAIN OUTPUTS VARIABLE(S): carbon, resp_hetero_soil
126!!
127!! REFERENCE(S)   :
128!! - Parton, W.J., D.S. Schimel, C.V. Cole, and D.S. Ojima. 1987. Analysis of
129!! factors controlling soil organic matter levels in Great Plains grasslands.
130!! Soil Sci. Soc. Am. J., 51, 1173-1179.
131!! - Gervois, S., P. Ciais, N. de Noblet-Ducoudre, N. Brisson, N. Vuichard,
132!! and N. Viovy (2008), Carbon and water balance of European croplands
133!! throughout the 20th century, Global Biogeochem. Cycles, 22, GB2022,
134!! doi:10.1029/2007GB003018.
135!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
136!! Environmental Modelling & Software (2011), doi:10.1016|j.envsoft.2011.05.016
137!!
138!! FLOWCHART    :
139!! \latexonly
140!! \includegraphics[scale=0.5]{soilcarbon_flowchart.jpg}
141!! \endlatexonly
142!! \n
143!_ ================================================================================================================================
144
145  SUBROUTINE soilcarbon (npts, dt, clay, &
146       soilcarbon_input, floodcarbon_input, control_temp_soil, control_moist_soil, &
147       carbon, &
148       resp_hetero_soil, resp_flood_soil, &
149       litter_above,litter_below,&
150       soilhum, DOC, moist_soil, DOC_EXP, &
151       lignin_struc_above, lignin_struc_below, &
152       floodout, runoff_per_soil, drainage_per_soil, wat_flux0, wat_flux,&
153       bulk_dens, soil_ph, poor_soils, veget_max, soil_mc, soiltile, &
154       DOC_to_topsoil, DOC_to_subsoil, flood_frac, &
155       precip2ground, precip2canopy, canopy2ground, &
156       dry_dep_canopy, DOC_precip2ground, DOC_precip2canopy, DOC_canopy2ground, &
157       DOC_infil, DOC_noinfil, interception_storage, biomass, fastr)
158
159!! 0. Variable and parameter declaration
160
161    !! 0.1 Input variables
162   
163    INTEGER(i_std), INTENT(in)                                        :: npts                  !! Domain size (unitless)
164    REAL(r_std), INTENT(in)                                           :: dt                    !! Time step \f$(dtradia one_day^{-1})$\f
165    REAL(r_std), DIMENSION(npts), INTENT(inout)                       :: clay                  !! Clay fraction (unitless, 0-1)
166    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)
167
168    REAL(r_std), DIMENSION(npts), INTENT(in)                          :: soil_ph               !! soil pH (pH unit, 0-14)
169    REAL(r_std), DIMENSION(npts), INTENT(in)                          :: poor_soils            !! Fraction of poor soils (0-1)
170    REAL(r_std), DIMENSION(npts,nbdl,npool*2), INTENT(in)             :: control_temp_soil     !! Temperature control of heterotrophic respiration (unitless: 0->1)
171    REAL(r_std), DIMENSION(npts,nbdl,nvm), INTENT(in)                 :: control_moist_soil    !! Moisture control of heterotrophic respiration (unitless: 0.25->1)
172
173    REAL(r_std), DIMENSION(npts,nlitt,nvm,nelements), INTENT(in)      :: litter_above          !! Metabolic and structural litter,above ground
174                                                                                               !! The unit is given by m^2 of ground
175                                                                                               !! @tex $(gC m^{-2})$ @endtex
176    REAL(r_std), DIMENSION(npts,nlitt,nvm,nbdl,nelements), INTENT(in) ::litter_below           !! Metabolic and structural litter, below ground
177                                                                                               !! The unit is given by m^2 of
178                                                                                               !! ground @tex $(gC m^{-2}) $ @endtex
179    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)                     :: soilhum               !! Daily soil humidity of each soil layer
180                                                                                               !! (unitless)
181    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                      :: lignin_struc_above    !! Ratio Lignin content in structural litter,
182                                                                                               !! above ground, (0-1, unitless)
183    REAL(r_std), DIMENSION(npts,nvm,nbdl), INTENT(in)                 :: lignin_struc_below    !! Ratio Lignin content in structural litter,
184                                                                                               !! below ground, (0-1, unitless)
185    REAL(r_std),DIMENSION (npts), INTENT (in)                         :: floodout              !! flux out of floodplains
186    REAL(r_std),DIMENSION (npts,nstm), INTENT (in)                    :: runoff_per_soil       !! Runoff per soil type [mm]
187    REAL(r_std),DIMENSION (npts,nstm), INTENT (in)                    :: drainage_per_soil     !! Drainage per soil type [mm]
188    REAL(r_std),DIMENSION (npts,nstm), INTENT(in)                     :: wat_flux0             !! Water flux in the first soil layers exported for soil C calculations(mm)
189    REAL(r_std),DIMENSION (npts,nslm,nstm), INTENT(in)                :: wat_flux              !! Water flux in the soil layers exported for soil C calculations(mm)
190    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                        :: veget_max             !! PFT "Maximal" coverage fraction of a PFT
191                                                                                               !! defined in the input vegetation map
192                                                                                               !! @tex $(m^2 m^{-2})$ @endtex
193                                                                                               !! The unit is given by m^2 of
194                                                                                               !! water @tex $(gC m{-2} of water)$ @endtex
195    REAL(r_std),DIMENSION (npts,nbdl,nstm), INTENT(in)                :: soil_mc               !! soil moisture content \f($m^3 \times m^3$)\f
196    REAL(r_std),DIMENSION (npts,nstm), INTENT (in)                    :: soiltile              !! Fraction of each soil tile (0-1, unitless)
197
198    REAL(r_std),DIMENSION (npts,nflow), INTENT(in)                    :: DOC_to_topsoil        !! DOC inputs to top of the soil column, from reinfiltration on
199                                                                                               !! floodplains and from irrigation
200                                                                                               !! @tex $(gC m^{-2} day{-1})$ @endtex
201    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                      :: precip2ground         !! Precipitation not intercepted by vegetation
202                                                                                               !! @tex $(mm.dt^{-1})$ @endtex
203    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                      :: precip2canopy         !! Precipitation intercepted by vegetation
204                                                                                               !! @tex $(mm.dt^{-1})$ @endtex
205    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                      :: canopy2ground         !! Water flux from canopy to ground
206                                                                                               !! @tex $(mm.dt^{-1})$ @endtex
207    REAL(r_std),DIMENSION (npts,nflow), INTENT(in)                    :: DOC_to_subsoil        !! DOC inputs to bottom of the soil column, from returnflow
208                                                                                               !! in swamps and lakes
209                                                                                               !! @tex $(gC m^{-2} day{-1})$ @endtex
210    REAL(r_std), DIMENSION(npts), INTENT(in)                          :: flood_frac            !! Flooded fraction of grid box
211                                                                                               !! @tex $(-)$ @endtex
212    REAL(r_std), DIMENSION(npts), INTENT(in)                          :: fastr                 !! Fast reservoir (mm)
213
214    !! 0.2 Output variables
215   
216    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                     :: resp_hetero_soil      !! Soil heterotrophic respiration \f$(gC m^{-2} (dtradia one_day^{-1})^{-1})$\f
217    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                     :: resp_flood_soil       !! Soil heterotrophic respiration when flooded
218                                                                                               !! \f$(gC m^{-2} (dtradia one_\day^{-1})^{-1})$\f
219    REAL(r_std), DIMENSION(npts,nvm,nexp,npool,nelements), INTENT(out):: DOC_EXP               !! Exported DOC through runoff, drainage, flood,
220                                                                                               !! The unit is given by m^2 of
221                                                                                               !! ground \f$(gC m^{-2} day^{-1})$\f
222    REAL(r_std), DIMENSION(npts,nbdl), INTENT(out)                    :: moist_soil            !! moisture in soil for one pixel (m3/m3)
223    !! TF-DOC, maybe put all togetehr into one matrix, note that all Throughfall related fluxes are expressed relative to the total area of the pixel
224    REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(out):: dry_dep_canopy        !! Increase in canopy storage of soluble OC & DOC
225                                                                                    !! @tex $(gC.m^{-2})$ @endtex
226    REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(out):: DOC_precip2canopy     !! Wet deposition of DOC onto canopy
227                                                                                    !! @tex $(gC.m^{-2})$ @endtex
228    REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(out):: DOC_precip2ground     !! Wet deposition of DOC not intecepted by canopy
229                                                                                    !! @tex $(gC.m^{-2})$ @endtex
230    REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(out):: DOC_canopy2ground     !! Wet deposition of DOC not intecepted by canopy
231                                                                                    !! @tex $(gC.m^{-2})$ @endtex
232    REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(out):: DOC_infil             !! Wet deposition of DOC infiltrating into the ground
233                                                                                    !! @tex $(gC.m^{-2})$ @endtex
234    REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(out):: DOC_noinfil           !! Wet deposition of DOC going into surface runoff
235                                                                                    !! @tex $(gC.m^{-2})$ @endtex 
236    !! 0.3 Modified variables
237
238    REAL(r_std), DIMENSION(npts,nvm,nbdl,npool,nelements), INTENT(inout) :: soilcarbon_input   !! Amount of carbon going into the carbon pools
239                                                                                               !! from litter decomposition \f$(gC m^{-2} day^{-1})$\f
240    REAL(r_std), DIMENSION(npts,nvm,npool,nelements), INTENT(inout)    :: floodcarbon_input    !! Amount of carbon going directly into the water column
241                                                                                          !! from litter decomposition \f$(gC m^{-2} day^{-1})$\f   
242    REAL(r_std), DIMENSION(npts,ncarb,nvm,nbdl), INTENT(inout)                :: carbon   !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f
243    REAL(r_std), DIMENSION(npts,nvm,nbdl,ndoc,npool,nelements), INTENT(inout) :: DOC      !! Dissolved Organic Carbon in soil
244                                                                                          !! The unit is given by m^2 of
245                                                                                          !! ground @tex $(gC m{-2} of ground)$ @endtex
246    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass           !! Biomass @tex $(gC m^{-2})$ @endtex
247    REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(inout):: interception_storage      !! Storage of DOC and solube OC in canopy
248                                                                                          !! @tex $(gC.m^{-2})$ @endtex 
249    !! 0.4 Local variables
250    REAL(r_std), SAVE, DIMENSION(ncarb)                           :: carbon_tau       !! Residence time in carbon pools (years)
251!$OMP THREADPRIVATE(carbon_tau)
252    REAL(r_std), DIMENSION(npts,ncarb,ncarb)                      :: frac_carb        !! Flux fractions between carbon pools
253                                                                                      !! (second index=origin, third index=destination)
254                                                                                      !! (unitless, 0-1)
255    REAL(r_std), DIMENSION(npts,ncarb,nelements,nbdl)             :: fluxtot          !! Total flux out of carbon pools \f$(gC m^{2})$\f
256    REAL(r_std), DIMENSION(npts,nvm,nbdl,npool)                   :: fluxtot_DOC      !! Total flux out of dissolved organic carbon pools \f$(gC m^{2}ground)$\f
257    REAL(r_std), DIMENSION(npts,ncarb,nelements,nbdl)             :: fluxtot_flood    !! Total flux out of carbon pools \f$(gC m^{2})$\f
258    REAL(r_std), DIMENSION(npts,nvm,nbdl,npool)                   :: fluxtot_DOC_flood!! Total flux out of dissolved organic carbon pools \f$(gC m^{2}ground\)$\f
259
260!! Add things for extrapolating passive SOC
261    CHARACTER(LEN=30)                                             :: var_name         !! To store variables with xios
262    CHARACTER(LEN=2)                                              :: part_str         !! To store variables with xios
263!!!!
264    REAL(r_std), DIMENSION(npts,npool,ncarb,nelements,nbdl)       :: flux             !! Fluxes between carbon pools \f$(gC m^{2})$\f
265    CHARACTER(LEN=7), DIMENSION(ncarb)                            :: carbon_str       !! Name of the carbon pools for informative outputs (unitless)
266    REAL(r_std), SAVE, DIMENSION(ncarb)                           :: priming_param    !! parmater controlling the sensitvity of priming effect on  each pools to the amount of LOM
267!$OMP THREADPRIVATE(priming_param)
268    REAL(r_std), SAVE, DIMENSION(npool)                           :: DOC_tau          !! Residence time of DOC (days)
269!$OMP THREADPRIVATE(DOC_tau)
270    REAL(r_std), DIMENSION(npts,nvm,nbdl)                         :: litter_tot       !! total litter (gC/m**2 of ground)
271    REAL(r_std), DIMENSION(npts,nvm,nbdl)                         :: LOM              !! total labile organic matter that could induce the decompostion rate
272                                                                                      !! (gC/m**2 of ground)
273    INTEGER(i_std)                                                :: k,kk,m,j,l,i,kkk !! Indices (unitless)
274    REAL(r_std),DIMENSION(0:nbdl)                                 :: z_soil           !! Soil levels (m)
275    REAL(r_std), DIMENSION(npts,nvm, nbdl,npool,nelements)        :: DOC_RE           !! Adsoprtion Desorption flux of DOC (mmol C/kg)
276    REAL(r_std), DIMENSION(npts,nvm, nbdl,npool,nelements)        :: DOC_FLUX         !! DOC flux between layers in f$(gC m^{-2} day^{-1})$\f
277    REAL(r_std), DIMENSION(npts,nvm, nbdl,npool,nelements)        :: DOC_FLUX_DIFF    !! DOC flux by diffusion between layers in f$(gC m^{-2} day^{-1})$\f
278    REAL(r_std), DIMENSION(npts,nvm, nbdl,npool,nelements)        :: DOC_FLUX_DIFF_old !! DOC flux by diffusion between layers in f$(gC m^{-2} day^{-1})$\f for storage
279    REAL(r_std), DIMENSION(npts,nvm,npool,nelements)              :: DOC_RUN          !! DOC lost with runoff in f$(gC m^{-2} day^{-1})$\f
280    REAL(r_std), DIMENSION(npts,nvm,npool,nelements)              :: DOC_FLOOD        !! DOC lost with flooding in f$(gC m^{-2} day^{-1})$\f
281    REAL(r_std), DIMENSION(npts,nvm,npool,nelements)              :: DOC_DRAIN        !! DOC lost with drainage in f$(gC m^{-2} day^{-1})$\f
282    REAL(r_std), DIMENSION(npts,nvm,nbdl,ndoc,npool,nelements)    :: DOC_old          !! Dissolved Organic Carbon in soil
283                                                                                      !! The unit is given by m^2 of
284                                                                                      !! ground @tex $(gC m{-2} of ground)$ @endtex variable used for storage
285    REAL(r_std), DIMENSION(npts,nvm,nbdl,ndoc,npool,nelements)    :: DOC_old2         !! Dissolved Organic Carbon in soil
286                                                                                      !! The unit is given by m^2 of
287                                                                                      !! ground @tex $(gC m{-2} of ground)$ @endtex variable used for storage
288    REAL(r_std), DIMENSION(npts,nvm,nbdl,ndoc,npool,nelements)    :: DOC_old_buffer   !! Dissolved Organic Carbon in soil
289                                                                                      !! The unit is given by m^2 of
290                                                                                      !! ground @tex $(gC m{-2} of ground)$ @endtex variable used for storage
291
292    REAL(r_std), DIMENSION(npts,ncarb,nvm,nbdl)                   :: carbon_old       !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f used for storage
293    REAL(r_std), DIMENSION(npts,ncarb,nvm,nbdl)                   :: carbon_old_buffer!! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f used for storage
294    REAL(r_std), DIMENSION(npts)                                  :: Qmax             !! Maximum adsorption capacity for DOC in the Langmuir equation (mg kg¿1)
295    REAL(r_std), DIMENSION(npts,nbdl)                             :: Kd               !! distribution coefficient between DOC adsorbed and DOC free (unitless [0-1])
296    REAL(r_std), DIMENSION(npts,nbdl)                             :: b                !! desorption term of the Initial Mass (IM) isotherm (mmol C/kg)
297    REAL(r_std), DIMENSION(npts,nbdl)                             :: IM               !! sorption affinity of the Initial Mass (IM) isotherm (unitless)
298
299    REAL(r_std), DIMENSION(npts)                                  :: Dif_coef         !! Diffusion coeficient used for bioturbation (m2 dt-1)
300
301    REAL(r_std), DIMENSION(npts)                                  :: Dif_DOC          !! Diffusion coeficient used for DOC diffusion (m2 dt-1)
302
303    REAL(r_std), DIMENSION(npts)                                  :: CUE_coef         !! Partition coefficient for DOC decomposed between CO2 and POC (unitless, 0-1)
304                                                                                   
305    REAL(r_std), DIMENSION(npts,ncarb,nvm,nbdl)                   :: carbon_flux      !! Soil carbon flux within pools towards the soil column\f$(gC m^{2} dt^{-1})$\f
306    REAL(r_std), DIMENSION(npts,ncarb,nvm,nbdl)                   :: carbon_flux_old  !! Soil carbon flux within pools towards the soil column\f$(gC m^{2} dt^{-1})$\f
307                                                                                      !! used for storage
308    REAL(r_std), DIMENSION(npts,nstm)                             :: soilwater_31mm   !! Soil water content of the first x layers (m)
309    REAL(r_std), DIMENSION(npts)                                  :: soilDOC_31mm     !! Soil DOC content of the first x layers (m)
310    REAL(r_std), DIMENSION(npts)                                  :: soilDOC_corr     !! Soil DOC content of the first x layers (m)
311    REAL(r_std), DIMENSION(npts)                                  :: flux_red         !! Flux reduction in upper meter of the soil
312    REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements)            :: check_intern     !! Contains the components of the internal
313                                                                                      !! mass balance chech for this routine
314                                                                                      !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
315    REAL(r_std), DIMENSION(npts,nvm,nelements)                    :: closure_intern   !! Check closure of internal mass balance
316
317    REAL(r_std), DIMENSION(npts,nvm,nelements)                    :: pool_start       !! Start and end pool of this routine
318                                                                                      !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
319    REAL(r_std), DIMENSION(npts,nvm,nelements)                    :: pool_end         !! Start and end pool of this routine
320                                                                                      !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
321    REAL(r_std), DIMENSION(npts,nvm,nelements)                    :: pool_end_after   !! Start and end pool of this routine after mass balance correction
322                                                                                      !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
323    REAL(r_std), DIMENSION(npts,nvm,nelements)                    :: wet_dep_ground    !! total wet deposition of DOC onto the ground
324    REAL(r_std), DIMENSION(npts,nvm,nelements)                    :: wet_dep_flood     !! total wet deposition of DOC onto the flooded
325    REAL(r_std), DIMENSION(npts)                                  :: bio_frac          !! total wet deposition of DOC onto the flooded
326    REAL(r_std), DIMENSION(npts)                                  :: fastr_corr        !! correction factor based on fast reservoir
327    REAL(r_std)                                                   :: lat_exp           !! factor used to disable soil-water export of DOC
328!_ ================================================================================================================================
329
330    !! bavard is the level of diagnostic information, 0 (none) to 4 (full)
331    IF (printlev>=3) WRITE(numout,*) 'Entering soilcarbon' 
332
333!! 1. Initializations
334!    clay(:) = 1.E-06
335    !! 1.0 Calculation of soil moisture
336    IF (lat_exp_doc) THEN
337       lat_exp = un
338    ELSE
339       lat_exp = zero
340    ENDIF
341
342    moist_soil(:,:) = zero
343      DO l = 1,nbdl
344          DO i = 1, nstm
345             moist_soil(:,l)= moist_soil(:,l) + soil_mc(:,l,i)*soiltile(:,i)
346          ENDDO
347      ENDDO
348      !! 1.1 soil levels
349       z_soil(0) = zero
350       z_soil(1:nbdl) = diaglev(1:nbdl)
351
352       !! DOC flux reduction factor, effects of poor soils
353       flux_red(:) = flux_red_sro * (un - poor_soils(:)) + poor_soils(:)
354
355       !! 1.2 TF-DOC: Initialize wet and dry deposition of DOC
356       DOC_precip2canopy(:,:,icarbon) = zero
357       DOC_precip2ground(:,:,icarbon) = zero
358       dry_dep_canopy(:,:,icarbon) = zero
359       bio_frac(:)=SUM(veget_max(:,2:13),DIM=2)
360
361       IF (ok_TF_DOC) THEN   
362          !! Calculate the DOC depoited onto the ground
363          !! TF Calculate thoughfall DOC
364
365          !! TF.1 Flux of DOC with precipitation onto the ground and onto the canopy (interception storage)
366          DOC_precip2canopy(:,2:13,icarbon) = precip2canopy(:,2:13) * conc_DOC_rain * 1e-3
367          DOC_precip2ground(:,2:13,icarbon) = precip2ground(:,2:13) * conc_DOC_rain * 1e-3
368          !! TF.2 Daily contirubitoon of dry deposition and exsuadtes and feces
369
370          DO j = 2,9! Loop over PFTs composed of trees
371             DO i = 1, npts! Loop over # of pixels
372                IF (veget_max(i,j) .GT. 0) THEN
373                   dry_dep_canopy(i,j,icarbon) = DOC_incr_per_LEAF_M2 * dt * biomass(i,j,ileaf,icarbon) * veget_max(i,j)
374                   !! So far, this is only the parametrisation for PFT=2
375                   !! sa the rate is per day, we have to multiply by dt
376                ELSE
377                   dry_dep_canopy(i,j,icarbon) = zero
378                ENDIF
379             ENDDO ! Loop over # of pixels
380          ENDDO ! # End Loop over # of PFTs
381       ENDIF !(ok_TF_DOC)
382
383    !! 1.4 Initialize check for mass balance closure
384    !  The mass balance is calculated at the end of this routine
385    !  in section 14. Initial biomass and harvest pool and all other
386    !  relevant pools were set to zero before calling this routine.
387    pool_start(:,:,:) = zero
388    IF (ld_doc) THEN
389       !Starting carbon pool
390       DO m = 1, nvm
391          DO i = 1,ncarb
392             DO l = 1, nbdl
393                pool_start(:,m,icarbon) = pool_start(:,m,icarbon) + &
394                     ( carbon(:,i,m,l) ) * veget_max(:,m)
395             ENDDO
396          ENDDO
397          pool_start(:,m,icarbon) = pool_start(:,m,icarbon) + interception_storage(:,m,icarbon)
398       ENDDO
399
400       !Starting DOC pool
401       DO m = 1, nvm
402          DO i = 1,npool
403             DO l = 1, nbdl
404                DO j = 1, ndoc
405                   pool_start(:,m,icarbon) = pool_start(:,m,icarbon) + &
406                        ( DOC(:,m,l,j,i,icarbon)) * veget_max(:,m)
407                ENDDO
408             ENDDO
409          ENDDO
410       ENDDO
411
412
413       !Starting litter below pool
414       DO m = 1,nvm
415          DO l = 1, nbdl
416             pool_start(:,m,icarbon) = pool_start(:,m,icarbon) + &
417                  (soilcarbon_input(:,m,l,imetbel,icarbon)+soilcarbon_input(:,m,l,istrbel,icarbon)+ &
418                  soilcarbon_input(:,m,l,imetabo,icarbon)+soilcarbon_input(:,m,l,istrabo,icarbon)) * &
419                  veget_max(:,m) * dt
420          ENDDO
421       ENDDO
422
423       !Starting litter above pool
424       DO m = 1,nvm
425          WHERE (veget_max(:,m) .GT. zero)
426             pool_start(:,m,icarbon) = pool_start(:,m,icarbon) + &
427                  (DOC_to_topsoil(:,idocl) + DOC_to_subsoil(:,idocl) + DOC_to_topsoil(:,idocr) + DOC_to_subsoil(:,idocr) + &
428                  floodcarbon_input(:,m,iact,icarbon) + floodcarbon_input(:,m,islo,icarbon))*dt*veget_max(:,m) + &
429                  DOC_precip2canopy(:,m,icarbon) + DOC_precip2ground(:,m,icarbon) + dry_dep_canopy(:,m,icarbon)
430          ELSEWHERE
431             pool_start(:,m,icarbon) = pool_start(:,m,icarbon) 
432          ENDWHERE
433       ENDDO
434    ENDIF
435
436    !! 1.2 Get soil "constants"
437    !! 1.2.1 Flux fractions between carbon pools: depend on clay content, recalculated each time
438    ! From active pool: depends on clay content
439    frac_carb(:,iactive,iactive) = zero
440    frac_carb(:,iactive,ipassive) = frac_carb_ap /(un - metabolic_ref_frac + active_to_pass_clay_frac*clay(:))
441    frac_carb(:,iactive,islow) = un - frac_carb(:,iactive,ipassive) 
442    ! 1.2.1.2 from slow pool
443    frac_carb(:,islow,islow) = zero
444    frac_carb(:,islow,iactive) = frac_carb_sa
445    frac_carb(:,islow,ipassive) = un - frac_carb(:,islow,iactive) 
446    ! From passive pool
447    frac_carb(:,ipassive,ipassive) = zero
448    frac_carb(:,ipassive,iactive) = frac_carb_pa
449    frac_carb(:,ipassive,islow) = un - frac_carb(:,ipassive,iactive)
450
451    ! 1.2.1.3 Initialization
452    carbon_flux(:,:,:,:) =  zero
453    carbon_flux_old(:,:,:,:) =  zero
454    IF ( firstcall ) THEN
455
456        !! 1.2.2 Residence times in carbon pools (in year converted in days)
457        carbon_tau(iactive) = carbon_tau_iactive * one_year      !residence times when SOM decomposition is function of litter coming from from optimization with the theoretical model SOM-FOM-DEP
458        carbon_tau(islow) = carbon_tau_islow * one_year          ! See Guenet et al., (In prep)
459        carbon_tau(ipassive) = carbon_tau_ipassive * one_year   
460
461        priming_param(iactive)=priming_param_iactive  !Priming parameter for mineralization obtained from optimization of SOM-FOM-DEP see Guenet et al., (In prep)
462        priming_param(islow)=priming_param_islow
463        priming_param(ipassive)=priming_param_ipassive
464
465        ! The values of DOC_tau are coming from Kalbitz et al 2003 Geoderma
466        ! (already in days)
467        DOC_tau(imetabo) = DOC_tau_labile
468        DOC_tau(istrabo) = DOC_tau_labile 
469        DOC_tau(imetbel) = DOC_tau_labile 
470        DOC_tau(istrbel) = DOC_tau_labile 
471        DOC_tau(iact) = DOC_tau_labile 
472        DOC_tau(islo) = DOC_tau_stable 
473        DOC_tau(ipas) = DOC_tau_stable 
474
475       !! 1.3 soil levels
476        !! 1.4 Messages : display the residence times 
477        carbon_str(iactive) = 'active'
478        carbon_str(islow) = 'slow'
479        carbon_str(ipassive) = 'passive'
480       
481        WRITE(numout,*) 'soilcarbon:'
482       
483        WRITE(numout,*) '   > minimal carbon residence time in carbon pools (d):'
484        DO k = 1, ncarb ! Loop over carbon pools
485          WRITE(numout,*) '(1, ::carbon_str(k)):',carbon_str(k),' : (1, ::carbon_tau(k)):',carbon_tau(k)
486        ENDDO
487       
488        WRITE(numout,*) '   > flux fractions between carbon pools: depend on clay content'
489        firstcall = .FALSE.
490       
491    ENDIF
492
493    !! 1.5 Initialization
494    !! 1.5.1 Initialization of parameters
495        Dif_coef(:) = (Dif/one_year)*dt
496        Dif_DOC(:)=D_DOC*dt        !We divided by one_year 
497                                   !because the unit of the Dif
498                                   !parameter is m2 yr-1 and the time stept of the model
499                                   !is dt (1/48, by default, i.e. the fraction
500                                   !of 30 min expressed in days)
501        CUE_coef(:)=CUE
502
503    !! 1.5.2 Set soil respiration and DOC fluxes to zero
504    resp_hetero_soil(:,:) = zero
505    resp_flood_soil(:,:) = zero
506    DOC_EXP(:,:,:,:,:) = zero
507    DOC_RE(:,:,:,:,:) = zero
508    DOC_FLUX(:,:,:,:,:) = zero
509    DOC_FLUX_DIFF(:,:,:,:,:) = zero
510    DOC_FLUX_DIFF_old(:,:,:,:,:) = zero
511    DOC_RUN(:,:,:,:)= zero
512    DOC_FLOOD(:,:,:,:)= zero
513    DOC_DRAIN(:,:,:,:)= zero
514    fluxtot(:,:,:,:)= zero
515    fluxtot_DOC(:,:,:,:) = zero
516    fluxtot_flood(:,:,:,:) = zero
517    fluxtot_DOC_flood(:,:,:,:) = zero
518    litter_tot(:,:,:) = zero
519    DOC_canopy2ground(:,:,:) = zero
520    wet_dep_flood(:,:,:) = zero
521    wet_dep_ground(:,:,:) = zero
522    LOM(:,:,:)= zero
523
524    ! bulk_density must be in kg m-3 and is generally given in g cm-3
525    WHERE (bulk_dens(:) .LT. 500)
526    bulk_dens(:) = bulk_dens(:)*1e3
527    ENDWHERE
528 
529    !!1.6 Calculating water content of the first 5 layers
530    soilwater_31mm(:,:)= zero
531    soilwater_31mm(:,:)= soilwater_31mm(:,:) + (z_soil(1) * soil_mc(:,1,:))   
532     IF (sro_bottom .GT. 1) THEN
533      DO l=2, sro_bottom
534         soilwater_31mm(:,:)= soilwater_31mm(:,:) + ((z_soil(l)-z_soil(l-1)) * soil_mc(:,l,:))
535      ENDDO
536    ENDIF
537
538    !! 2. Fluxes between carbon reservoirs, and to the atmosphere (respiration) \n
539    !! 2.1 calculation of TF fluxes
540    !! Calculation of TF-Fluxes
541
542       !! TF.3 Calculate the flow to the ground
543
544       !! Calculate max possible DOC flux, assuming a maximum concentration
545       DO j = 2,nvm! Loop over # of PFTs
546          DO i = 1, npts! Loop over # of pixels
547             IF (veget_max(i,j) > zero) THEN
548                DOC_canopy2ground(i,j,icarbon) = min(canopy2ground(i,j)*conc_DOC_max * 1e-3, &
549                     &                                              interception_storage(i,j,icarbon))
550             ELSE
551                DOC_canopy2ground(i,j,icarbon) = zero
552             ENDIF
553          ENDDO ! Loop over # of pixels
554       ENDDO ! # End Loop over # of PFTs
555
556       DO j = 1,nvm! Loop over # of PFTs
557          wet_dep_flood(:,j,icarbon)  = (DOC_precip2ground(:,j,icarbon) + DOC_canopy2ground(:,j,icarbon)) * flood_frac(:)
558          wet_dep_ground(:,j,icarbon) = (DOC_precip2ground(:,j,icarbon) + DOC_canopy2ground(:,j,icarbon)) * (un-flood_frac(:))
559       ENDDO
560
561       IF (.NOT. lat_exp_doc) THEN
562          wet_dep_ground = wet_dep_ground + wet_dep_flood
563          wet_dep_flood = zero
564       ELSE
565          !Do nothing
566       ENDIF
567
568       !! TF.4 Update the canopy storage of soluble OC and DOC
569
570       DO j = 2,nvm! Loop over # of PFTs
571          DO i = 1, npts! Loop over # of pixels
572             IF (veget_max(i,j) > zero) THEN
573                interception_storage(i,j,icarbon) = interception_storage(i,j,icarbon) + dry_dep_canopy(i,j,icarbon) + &
574                     &                          DOC_precip2canopy(i,j,icarbon) - DOC_canopy2ground(i,j,icarbon)
575             ELSE
576                interception_storage(i,j,icarbon) = zero
577             ENDIF
578          ENDDO ! Loop over # of pixels
579       ENDDO ! # End Loop over # of PFTs
580
581       !! TF.5 Calculate the fraction that is infiltrating into the soil
582       DOC_infil = zero
583       DOC_noinfil = zero
584
585    !! 2.2 Input from litter to DOC and from reinflitration in floodplains
586    !! For the variables DOC_to_topsoil and DOC_to_subsoil we multiplied by
587    !veget_max because it is not defined per PFT and therefore we redistribute
588    !it following veget_max. Moreover we assume the all the DOC coming from
589    !floodplains reinfiltration is going to DOC coming from the active pool.
590    DO m = 2, nvm
591       DO l = 1,nbdl
592
593          DOC(:,m,l,ifree,imetabo,icarbon) = DOC(:,m,l,ifree,imetabo,icarbon) + soilcarbon_input(:,m,l,imetabo,icarbon)*dt 
594          DOC(:,m,l,ifree,istrabo,icarbon) = DOC(:,m,l,ifree,istrabo,icarbon) + soilcarbon_input(:,m,l,istrabo,icarbon)*dt
595          DOC(:,m,l,ifree,imetbel,icarbon) = DOC(:,m,l,ifree,imetbel,icarbon) + soilcarbon_input(:,m,l,imetbel,icarbon)*dt
596          DOC(:,m,l,ifree,istrbel,icarbon) = DOC(:,m,l,ifree,istrbel,icarbon) + soilcarbon_input(:,m,l,istrbel,icarbon)*dt
597
598       ENDDO
599          WHERE (veget_max(:,m) .GT. zero)
600          DOC(:,m,1,ifree,iact,icarbon) = DOC(:,m,1,ifree,iact,icarbon) + DOC_to_topsoil(:,idocl)*dt/bio_frac(:) &
601               + wet_dep_ground(:,m,icarbon) * undemi/veget_max(:,m)
602          DOC(:,m,nbdl,ifree,iact,icarbon) = DOC(:,m,nbdl,ifree,iact,icarbon) + DOC_to_subsoil(:,idocl)*dt/bio_frac(:)
603          DOC(:,m,1,ifree,islo,icarbon) = DOC(:,m,1,ifree,islo,icarbon) + DOC_to_topsoil(:,idocr)*dt/bio_frac(:) &
604               + wet_dep_ground(:,m,icarbon) * undemi/veget_max(:,m)
605          DOC(:,m,nbdl,ifree,islo,icarbon) = DOC(:,m,nbdl,ifree,islo,icarbon) + DOC_to_subsoil(:,idocr)*dt/bio_frac(:)
606          ELSEWHERE
607          DOC(:,m,1,ifree,iact,icarbon) = DOC(:,m,1,ifree,iact,icarbon) 
608          DOC(:,m,nbdl,ifree,iact,icarbon) = DOC(:,m,nbdl,ifree,iact,icarbon)
609          ENDWHERE
610    ENDDO
611
612    !! 2.3. Calculate fluxes
613    DO m = 2,nvm
614      !calculate total litter available
615      litter_tot(:,m,1)=litter_above(:,imetabolic,m,icarbon)+litter_above(:,istructural,m,icarbon)
616   
617      DO l = 1, nbdl
618      litter_tot(:,m,l)=litter_tot(:,m,l) + litter_below(:,imetabolic,m,l,icarbon)+litter_below(:,istructural,m,l,icarbon)
619      ENDDO
620   ENDDO
621
622    DO m = 1,nvm ! Loop over # PFTs
623
624
625      !! 2.3.1. Flux out of pools
626
627      DO k = 1, ncarb ! Loop over carbon pools from which the flux comes
628
629        DO l = 1, nbdl !Loop over soil layers
630
631         ! Determine total flux out of pool
632         ! S.L. Piao 2006/05/05 - for crop multiply tillage factor of decomposition
633         ! Not crop
634         
635
636          IF (k == iactive) THEN
637             LOM(:,m,l) = litter_tot(:,m,l) + DOC(:,m,l,ifree,imetabo,icarbon) + DOC(:,m,l,ifree,istrabo,icarbon) + DOC(:,m,l,ifree,imetbel,icarbon) + DOC(:,m,l,ifree,istrbel,icarbon)
638          ELSEIF (k == islow) THEN
639             LOM(:,m,l) = litter_tot(:,m,l) + carbon(:,iactive,m,l) + DOC(:,m,l,ifree,iact,icarbon) + &
640                          DOC(:,m,l,ifree,imetabo,icarbon) + DOC(:,m,l,ifree,istrabo,icarbon) + DOC(:,m,l,ifree,imetbel,icarbon) + DOC(:,m,l,ifree,istrbel,icarbon)
641          ELSEIF (k == ipassive) THEN
642             LOM(:,m,l) = litter_tot(:,m,l) + carbon(:,iactive,m,l) + carbon(:,islow,m,l) + DOC(:,m,l,ifree,iact,icarbon) + DOC(:,m,l,ifree,islo,icarbon) + &
643                          DOC(:,m,l,ifree,imetabo,icarbon) + DOC(:,m,l,ifree,istrabo,icarbon) + DOC(:,m,l,ifree,imetbel,icarbon) + DOC(:,m,l,ifree,istrbel,icarbon)
644          ELSE
645             WRITE(numout,*) "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
646             WRITE(numout,*) "BE CAREFUL SOMETHING STRANGE HAPPENS"
647             WRITE(numout,*) "IN STOMATE_SOILCARBON FOR LOM CALCULATION"
648             WRITE(numout,*) "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
649
650          ENDIF
651             
652       WHERE ((soil_mc(:,l,pref_soil_veg(m)) .GT. zero) .AND. (natural(m)))
653               fluxtot(:,k,icarbon,l) =(dt/(carbon_tau(k)*(un+poor_soils(:))))*(un-exp(-priming_param(k)*LOM(:,m,l))) * carbon(:,k,m,l) * &
654                    control_moist_soil(:,l,m) * control_temp_soil(:,l,k+4) * (un-flood_frac(:))
655               fluxtot_flood(:,k,icarbon,l) =(dt/(carbon_tau(k)*(un+poor_soils(:))))*(un-exp(-priming_param(k)*LOM(:,m,l))) * carbon(:,k,m,l) * &
656                    control_moist_soil(:,l,m) * control_temp_soil(:,l,k+4) * flood_frac(:) * un/trois
657            ! C3 crop
658       ELSEWHERE ((soil_mc(:,l,pref_soil_veg(m)) .GT. zero) .AND.( (.NOT. natural(m)) .AND. (.NOT. is_c4(m)) ))
659               fluxtot(:,k,icarbon,l) =(dt/(carbon_tau(k)*(un+poor_soils(:))))*(un-exp(-priming_param(k)*LOM(:,m,l))) * carbon(:,k,m,l) * &
660                    control_moist_soil(:,l,m) * control_temp_soil(:,l,k+4) * flux_tot_coeff(1) * (un-flood_frac(:))
661               fluxtot_flood(:,k,icarbon,l) =(dt/(carbon_tau(k)*(un+poor_soils(:))))*(un-exp(-priming_param(k)*LOM(:,m,l))) * carbon(:,k,m,l) * &
662                    control_moist_soil(:,l,m) * control_temp_soil(:,l,k+4) * flux_tot_coeff(1) * flood_frac(:) * un/trois
663            ! C4 crop
664       ELSEWHERE ((soil_mc(:,l,pref_soil_veg(m)) .GT. zero) .AND.( (.NOT. natural(m)) .AND. (is_c4(m)) )) 
665               fluxtot(:,k,icarbon,l) =(dt/(carbon_tau(k)*(un+poor_soils(:))))*(un-exp(-priming_param(k)*LOM(:,m,l))) * carbon(:,k,m,l) * &
666                    control_moist_soil(:,l,m) * control_temp_soil(:,l,k+4) * flux_tot_coeff(2) * (un-flood_frac(:))
667               fluxtot_flood(:,k,icarbon,l) =(dt/(carbon_tau(k)*(un+poor_soils(:))))*(un-exp(-priming_param(k)*LOM(:,m,l))) * carbon(:,k,m,l) * &
668                    control_moist_soil(:,l,m) * control_temp_soil(:,l,k+4) * flux_tot_coeff(2) * flood_frac(:) * un/trois
669       ELSEWHERE
670               fluxtot(:,k,icarbon,l) = zero
671               fluxtot_flood(:,k,icarbon,l) = zero
672       ENDWHERE
673
674          ! END - S.L. Piao 2006/05/05 - for crop multiply tillage factor of decomposition
675
676          ! Carbon flux from active pools depends on clay content
677          IF ( k .EQ. iactive ) THEN
678              fluxtot(:,k,icarbon,l) = fluxtot(:,k,icarbon,l) * ( un - flux_tot_coeff(3) * clay(:) )
679              fluxtot_flood(:,k,icarbon,l) = fluxtot_flood(:,k,icarbon,l) * ( un - flux_tot_coeff(3) * clay(:) )
680          ENDIF
681        ENDDO ! End of loop over soil layers
682
683      ENDDO ! End of loop over carbon pools
684
685      IF (.NOT. lat_exp_doc) THEN
686         fluxtot = fluxtot + fluxtot_flood
687         fluxtot_flood = zero
688      ELSE
689         !do nothing
690      ENDIF
691
692      !! 2.3.2 respiration
693      !BE CAREFUL: Here resp_hetero_soil is divided by dt to have a value which corresponds to
694      ! the sechiba time step but then in stomate.f90 resp_hetero_soil is multiplied by dt.
695      ! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
696      ! time step and avoid some constructions that could create bug during future developments.
697      !
698      DO l = 1,nbdl !Loop over soil layers
699         resp_hetero_soil(:,m) = resp_hetero_soil(:,m)+&
700            ( (un-CUE_coef(:)) * fluxtot(:,iactive,icarbon,l) + &
701            (un-CUE_coef(:)) * fluxtot(:,islow,icarbon,l) + &
702            (un-CUE_coef(:))* fluxtot(:,ipassive,icarbon,l)  ) / dt
703         resp_flood_soil(:,m) = resp_flood_soil(:,m)+&
704            ( (un-CUE_coef(:)) * fluxtot_flood(:,iactive,icarbon,l) + &
705            (un-CUE_coef(:)) * fluxtot_flood(:,islow,icarbon,l) + &
706            (un-CUE_coef(:))* fluxtot_flood(:,ipassive,icarbon,l)  ) / dt
707      ENDDO !End loop over soil layers
708
709      !! 2.3.3 DOC calculation
710      !! 2.3.3.1 DOC decomposition
711      ! In the Kalbitz et al., 2003 Geoderma as well
712      ! as in Qualls and Haines 1992 SSAJ
713      ! the equation are defined with exponential but for the time
714      ! step we used we can approximate the values with 1st order kinetics.
715
716     DO j = 1,npool  !Loop over pools
717       DO l = 1,nbdl !Loop over soil layers
718        WHERE (soil_mc(:,l,pref_soil_veg(m)) .GT. zero) 
719                 fluxtot_DOC(:,m,l,j) = (dt/(DOC_tau(j)*(un+poor_soils(:))))* (un - flood_frac(:)) * &
720                                         DOC(:,m,l,ifree,j,icarbon)*control_temp_soil(:,l,npool+j) 
721                 fluxtot_DOC_flood(:,m,l,j) = (dt/(DOC_tau(j)*(un+poor_soils(:))))* flood_frac(:) * un/trois * &
722                                         DOC(:,m,l,ifree,j,icarbon)*control_temp_soil(:,l,npool+j)
723         ELSEWHERE
724                 fluxtot_DOC(:,m,l,j) = zero
725                 fluxtot_DOC_flood(:,m,l,j) = zero
726         ENDWHERE
727
728          resp_hetero_soil(:,m) = resp_hetero_soil(:,m)+&
729             (un-CUE_coef(:))*fluxtot_DOC(:,m,l,j)/dt
730
731          resp_flood_soil(:,m) = resp_flood_soil(:,m)+&
732             (un-CUE_coef(:))*fluxtot_DOC_flood(:,m,l,j)/dt
733
734      ENDDO !End loop over soil layers!
735     ENDDO !End loop over soil pool!
736
737      !! 2.3.4 Pools update
738      !! 2.3.4.1 For POC
739      DO l= 1, nbdl ! Loop over soil layers
740
741       carbon(:,iactive,m,l) = carbon(:,iactive,m,l)+ &
742                               frac_carb(:,ipassive,iactive)*(CUE_coef(:))*fluxtot_DOC(:,m,l,ipas) + &
743                               frac_carb(:,islow,iactive)*(CUE_coef(:))*fluxtot_DOC(:,m,l,islo) + &
744                               (CUE_coef(:))*fluxtot_DOC(:,m,l,imetbel) + &
745                               (CUE_coef(:))*fluxtot_DOC(:,m,l,istrbel) * (un - lignin_struc_below(:,m,l)) + &
746                               (CUE_coef(:))*fluxtot_DOC(:,m,l,imetabo) + &
747                               (CUE_coef(:))*fluxtot_DOC(:,m,l,istrabo) * (un - lignin_struc_above(:,m)) + &
748                               frac_carb(:,ipassive,iactive)*(CUE_coef(:))*fluxtot_DOC_flood(:,m,l,ipas) + &
749                               frac_carb(:,islow,iactive)*(CUE_coef(:))*fluxtot_DOC_flood(:,m,l,islo) + &
750                               (CUE_coef(:))*fluxtot_DOC_flood(:,m,l,imetbel) + &
751                               (CUE_coef(:))*fluxtot_DOC_flood(:,m,l,istrbel) * (un - lignin_struc_below(:,m,l)) + &
752                               (CUE_coef(:))*fluxtot_DOC_flood(:,m,l,imetabo) + &
753                               (CUE_coef(:))*fluxtot_DOC_flood(:,m,l,istrabo) * (un - lignin_struc_above(:,m))
754
755       carbon(:,islow,m,l) = carbon(:,islow,m,l)+ &
756                               frac_carb(:,ipassive,islow)*(CUE_coef(:))*fluxtot_DOC(:,m,l,ipas) + &
757                               frac_carb(:,iactive,islow)*(CUE_coef(:))*fluxtot_DOC(:,m,l,iact) + &
758                               (CUE_coef(:))*fluxtot_DOC(:,m,l,istrbel) *lignin_struc_below(:,m,l) + &
759                               (CUE_coef(:))*fluxtot_DOC(:,m,l,istrabo) *lignin_struc_above(:,m) + &
760                               frac_carb(:,ipassive,islow)*(CUE_coef(:))*fluxtot_DOC_flood(:,m,l,ipas) + &
761                               frac_carb(:,iactive,islow)*(CUE_coef(:))*fluxtot_DOC_flood(:,m,l,iact) + &
762                               (CUE_coef(:))*fluxtot_DOC_flood(:,m,l,istrbel) *lignin_struc_below(:,m,l) + &
763                               (CUE_coef(:))*fluxtot_DOC_flood(:,m,l,istrabo) *lignin_struc_above(:,m)
764
765       carbon(:,ipassive,m,l) = carbon(:,ipassive,m,l)+ &
766                               frac_carb(:,iactive,ipassive)*(CUE_coef(:))*fluxtot_DOC(:,m,l,iact) + &
767                               frac_carb(:,islow,ipassive)*(CUE_coef(:))*fluxtot_DOC(:,m,l,islo) + &
768                               frac_carb(:,iactive,ipassive)*(CUE_coef(:))*fluxtot_DOC_flood(:,m,l,iact) + &
769                               frac_carb(:,islow,ipassive)*(CUE_coef(:))*fluxtot_DOC_flood(:,m,l,islo)
770
771      ENDDO ! End loop over soil layers
772 
773      DO k = 1, ncarb ! Loop over carbon pools from which the flux comes
774        DO l = 1, nbdl !Loop over soil layers
775
776         carbon(:,k,m,l) = carbon(:,k,m,l) - fluxtot(:,k,icarbon,l) - fluxtot_flood(:,k,icarbon,l)
777
778        ENDDO ! End of loop over soil layers
779      ENDDO ! End of loop over carbon pools
780
781      !! 2.3.4.2 For DOC
782      DO l = 1,nbdl !Loop over soil layers
783         DO kk = 1,nelements !Loop over elements
784           IF (l .GT. sro_bottom) THEN
785              !For Active pool
786              DOC(:,m,l,ifree,iact,kk) = DOC(:,m,l,ifree,iact,kk) + (CUE_coef(:)) * fluxtot(:,iactive,kk,l) + &
787                   (CUE_coef(:)) * fluxtot_flood(:,iactive,kk,l)
788             
789              !For Slow pool   
790              DOC(:,m,l,ifree,islo,kk) = DOC(:,m,l,ifree,islo,kk) + (CUE_coef(:)) * fluxtot(:,islow,kk,l) + &
791                   (CUE_coef(:)) * fluxtot_flood(:,islow,kk,l)
792
793              !For Passive pool
794              DOC(:,m,l,ifree,ipas,kk) = DOC(:,m,l,ifree,ipas,kk) + (CUE_coef(:)) * fluxtot(:,ipassive,kk,l) + &
795                   (CUE_coef(:)) * fluxtot_flood(:,ipassive,kk,l)
796           ELSE
797              !For Active pool
798              DOC(:,m,l,ifree,iact,kk) = DOC(:,m,l,ifree,iact,kk) + (CUE_coef(:)) * fluxtot(:,iactive,kk,l)
799              !For Slow pool
800              DOC(:,m,l,ifree,islo,kk) = DOC(:,m,l,ifree,islo,kk) + (CUE_coef(:)) * fluxtot(:,islow,kk,l) 
801              !For Passive pool
802              DOC(:,m,l,ifree,ipas,kk) = DOC(:,m,l,ifree,ipas,kk) + (CUE_coef(:)) * fluxtot(:,ipassive,kk,l)
803           ENDIF
804        ENDDO !End loop over elements
805     ENDDO !End loop over soil layers
806
807     DO j = 1,npool  !Loop over pools
808       DO l = 1,nbdl !Loop over soil layers
809
810          DOC(:,m,l,ifree,j,icarbon) = DOC(:,m,l,ifree,j,icarbon) - fluxtot_DOC(:,m,l,j) - fluxtot_DOC_flood(:,m,l,j)
811
812       ENDDO !End loop over soil layers!
813     ENDDO !End loop over soil pool!
814
815      !! 2.3.5 Adsorption/resorption of DOC using linear Mass isotherm assuming
816      !that DOC from native soil is zero (b parameters in Nodvin et al., 1986,
817      !Soil Science)
818      !For the moment, fixed parameters taken from Kaiser and others,1996, b in mg/g soil
819     Kd(:,:)=kd_ads
820     Kd(:,1) = dix**(-3.1+0.2*log10(clay(:)*cent)+log10(z_soil(1)*cent/deux))
821     DO l = 2,nbdl
822        ! IM, b or Kd calculating following the statistical model obtained in
823        ! Camino Serrano et al. in prep
824        Kd(:,l) = dix**(-3.1+0.2*log10(clay(:)*cent)+log10((z_soil(l)+z_soil(l-1))/deux))
825        ! Kd(:,l) = 0.001225841-0.000212346*soil_ph(:)+0.003739918*clay(:)
826     ENDDO
827     Kd(:,:)=MAX(Kd(:,:),dix**(-3.2))
828     Kd(:,:)=MIN(Kd(:,:),dix**(-1.5))
829
830     ! Calculation of amount adsorbed/desorbed based on
831     ! Ota et al., 2013 JGR
832     DO l = 1,nbdl !Loop over soil layers
833        DO kk = 1,nelements !Loop over elements
834           DO kkk = 1,npool !Loop over soil and litter pool
835              DOC_RE(:,m,l,kkk,kk) = (Kd(:,l)*bulk_dens(:))/(Kd(:,l)*bulk_dens(:)+un)* &
836                   (DOC(:,m,l,ifree,kkk,kk)+DOC(:,m,l,iadsorbed,kkk,kk)) 
837           ENDDO !End loop over soil and litter pool
838        ENDDO !End loop over elements
839     ENDDO !End loop over soil layers
840
841     DOC_old(:,m,:,:,:,:)= zero
842     DOC_old(:,m,:,:,:,:)=DOC(:,m,:,:,:,:)
843     DO l = 1,nbdl !Loop over soil layers
844        DO kk = 1,nelements !Loop over elements
845           DO kkk = 1,npool !Loop over soil and litter pool
846              WHERE (((DOC_RE(:,m,l,kkk,kk)-DOC_old(:,m,l,iadsorbed,kkk,kk)) .GT. zero) .AND. (abs(DOC_RE(:,m,l,kkk,kk)-DOC_old(:,m,l,iadsorbed,kkk,kk)) .GT. DOC_old(:,m,l,ifree,kkk,kk)))
847                 !To avoid adsorbing more DOC than existing (when there is not enough free carbon)
848                 DOC(:,m,l,iadsorbed,kkk,kk)= DOC_old(:,m,l,iadsorbed,kkk,kk)+ DOC_old(:,m,l,ifree,kkk,kk)
849                 DOC(:,m,l,ifree,kkk,kk)= zero           
850              ELSEWHERE (((DOC_RE(:,m,l,kkk,kk)-DOC_old(:,m,l,iadsorbed,kkk,kk)) .LT. zero) .AND. (abs(DOC_RE(:,m,l,kkk,kk)-DOC_old(:,m,l,iadsorbed,kkk,kk)) .GT. DOC_old(:,m,l,iadsorbed,kkk,kk)))
851               !To avoid desorbing more DOC than existing (when there is not enough adsorbed carbon)
852                 DOC(:,m,l,iadsorbed,kkk,kk)= zero 
853                 DOC(:,m,l,ifree,kkk,kk)= DOC_old(:,m,l,iadsorbed,kkk,kk)+ DOC_old(:,m,l,ifree,kkk,kk)       
854              ELSEWHERE ((DOC_RE(:,m,l,kkk,kk)-DOC_old(:,m,l,iadsorbed,kkk,kk)) .NE. zero) 
855                !The normal functioning of the distribution coefficient to reach equilibrium between DOC adsorbed and DOC free   
856                 DOC(:,m,l,ifree,kkk,kk)=DOC_old(:,m,l,ifree,kkk,kk)-(DOC_RE(:,m,l,kkk,kk)-DOC_old(:,m,l,iadsorbed,kkk,kk)) 
857                 DOC(:,m,l,iadsorbed,kkk,kk)= DOC_old(:,m,l,iadsorbed,kkk,kk)+(DOC_RE(:,m,l,kkk,kk)-DOC_old(:,m,l,iadsorbed,kkk,kk))
858              ELSEWHERE
859                 DOC(:,m,l,ifree,kkk,kk)=DOC_old(:,m,l,ifree,kkk,kk)
860                 DOC(:,m,l,iadsorbed,kkk,kk)= DOC_old(:,m,l,iadsorbed,kkk,kk)
861              ENDWHERE
862           ENDDO !End loop over soil and litter pool
863        ENDDO !End loop over elements
864     ENDDO !End loop over soil layers
865
866
867
868     !! 2.3.6 Transport of DOC with the water fluxes following Futter et al.,
869     ! (2007) Water ressources research.
870     
871     !! 2.3.6.1 Change DOC free from gC m-2 of soil to g C m^-3 of water 
872     DO kk = 1, nelements !Loop over elements
873       
874        DO k = 1, npool ! Loop over soil and litter pools
875           
876           WHERE( soil_mc(:,1,pref_soil_veg(m)) .GT. zero)
877
878              DOC(:,m,1,ifree,k,kk) = DOC(:,m,1,ifree,k,kk) / (z_soil(1)*soil_mc(:,1,pref_soil_veg(m)))
879
880           ELSEWHERE
881              DOC(:,m,1,ifree,k,kk) = DOC(:,m,1,ifree,k,kk)
882           ENDWHERE
883
884        ENDDO ! End loop over soil and litter pools
885
886     ENDDO !End loop over elements
887
888     DO l = 2,nbdl !Loop over soil layers
889
890        DO kk = 1, nelements !Loop over elements
891
892           DO k = 1, npool ! Loop over soil and litter pools
893
894              WHERE( soil_mc(:,l,pref_soil_veg(m)) .GT. zero)
895
896                 DOC(:,m,l,ifree,k,kk) = DOC(:,m,l,ifree,k,kk) / ((z_soil(l)-z_soil(l-1))*soil_mc(:,l,pref_soil_veg(m)))
897
898              ELSEWHERE
899                 DOC(:,m,l,ifree,k,kk) = DOC(:,m,l,ifree,k,kk)
900              ENDWHERE
901
902           ENDDO ! End loop over soil and litter pools
903
904        ENDDO !End loop over elements
905
906     ENDDO !End loop over soil layers
907     ! DOC fluxes between layers calculated (we multiplied by
908     ! 1E-3 to convert from kg/m2 to m3/m2)
909     DO l = 1,nbdl-1 !Loop over soil layers
910        DO kk = 1,nelements !Loop over elements
911           DO kkk = 1,npool !Loop over soil and litter pool
912              WHERE ((wat_flux(:,l,pref_soil_veg(m)) .GT. 0.) .AND. (soil_mc(:,l,pref_soil_veg(m)) .GT. zero))
913                 DOC_FLUX(:,m,l,kkk,kk) = DOC(:,m,l,ifree,kkk,kk)*wat_flux(:,l,pref_soil_veg(m))*1E-3*flux_red(:)
914              ELSEWHERE ((wat_flux(:,l,pref_soil_veg(m)) .LT. 0.) .AND. (soil_mc(:,l,pref_soil_veg(m)) .GT. zero))
915                 DOC_FLUX(:,m,l,kkk,kk) = DOC(:,m,l+1,ifree,kkk,kk)*abs(wat_flux(:,l,pref_soil_veg(m)))*1E-3*flux_red(:)
916              ELSEWHERE 
917                 DOC_FLUX(:,m,l,kkk,kk) = zero
918              ENDWHERE
919           ENDDO !End loop over soil and litter pool
920        ENDDO !End loop over elements
921     ENDDO !End loop over soil layers
922
923     !For the last layer
924     DO kk = 1,nelements !Loop over elements
925        DO kkk = 1,npool !Loop over soil and litter pool
926           WHERE ((wat_flux(:,nbdl,pref_soil_veg(m)) .GT. 0.) .AND. (soil_mc(:,nbdl,pref_soil_veg(m)) .GT. zero))
927              DOC_FLUX(:,m,nbdl,kkk,kk) = DOC(:,m,nbdl,ifree,kkk,kk)*wat_flux(:,nbdl,pref_soil_veg(m))*1E-3*flux_red(:)
928           ELSEWHERE 
929              DOC_FLUX(:,m,nbdl,kkk,kk) = zero
930           ENDWHERE
931        ENDDO !End loop over soil and litter pool
932     ENDDO !End loop over elements
933
934     !! 2.3.7 DOC is reconverted from gC m^⁻3 of water in gC m^-2 of ground to exchange it with the other carbon pools
935     !For the first layer
936
937     DO k = 1, npool ! Loop over soil and litter pools
938        DO kk = 1,nelements !Loop over elements
939           WHERE(soil_mc(:,1,pref_soil_veg(m)) .GT. zero)             
940              DOC(:,m,1,ifree,k,kk)=DOC(:,m,1,ifree,k,kk) * soil_mc(:,1,pref_soil_veg(m)) * z_soil(1)
941           ELSEWHERE
942              DOC(:,m,1,ifree,k,kk)=DOC(:,m,1,ifree,k,kk) 
943           ENDWHERE
944        ENDDO !End loop over elements
945     ENDDO ! End loop over soil and litter pools
946
947
948     !For the other layers
949
950     DO l= 2, nbdl ! Loop over soil layers
951        DO k = 1, npool ! Loop over soil and litter pools
952           DO kk = 1,nelements !Loop over elements
953              WHERE(soil_mc(:,l,pref_soil_veg(m)) .GT. zero)
954                 DOC(:,m,l,ifree,k,kk)=DOC(:,m,l,ifree,k,kk) * soil_mc(:,l,pref_soil_veg(m)) * (z_soil(l)-z_soil(l-1))
955              ELSEWHERE
956                 DOC(:,m,l,ifree,k,kk)=DOC(:,m,l,ifree,k,kk) 
957              ENDWHERE
958           ENDDO !End loop over elements
959        ENDDO ! End loop over soil and litter pools
960     ENDDO ! End loop over soil layers
961
962     ! DOC update
963     !For the first layer
964     DOC_old(:,m,:,:,:,:)= zero
965     DOC_old(:,m,:,:,:,:)=DOC(:,m,:,:,:,:)
966
967
968     DO kk = 1,nelements !Loop over elements
969        DO kkk = 1,npool !Loop over soil and litter pool
970             
971           WHERE (wat_flux(:,1,pref_soil_veg(m)) .GT. 0.) 
972              DOC(:,m,1,ifree,kkk,kk) = DOC_old(:,m,1,ifree,kkk,kk) - MIN(DOC_old(:,m,1,ifree,kkk,kk),DOC_FLUX(:,m,1,kkk,kk)) 
973           ELSEWHERE (wat_flux(:,1,pref_soil_veg(m)) .LT. 0.) 
974              DOC(:,m,1,ifree,kkk,kk) = DOC_old(:,m,1,ifree,kkk,kk) + MIN(DOC_old(:,m,2,ifree,kkk,kk),DOC_FLUX(:,m,1,kkk,kk))
975           ELSEWHERE
976              DOC(:,m,1,ifree,kkk,kk) = DOC_old(:,m,1,ifree,kkk,kk)
977           ENDWHERE
978        ENDDO !End loop over soil and litter pool
979     ENDDO !End loop over elements
980
981     ! For the other layers
982     DO l = 2,nbdl-1 !Loop over soil layers
983        DO kk = 1,nelements !Loop over elements
984           DO kkk = 1,npool !Loop over soil and litter pool
985              WHERE ((wat_flux(:,l,pref_soil_veg(m)) .GT. 0.) .AND. (wat_flux(:,l-1,pref_soil_veg(m)) .GT. 0.))
986                 DOC(:,m,l,ifree,kkk,kk) = DOC_old(:,m,l,ifree,kkk,kk) - MIN(DOC_FLUX(:,m,l,kkk,kk),DOC_old(:,m,l,ifree,kkk,kk)) + MIN(DOC_FLUX(:,m,l-1,kkk,kk),DOC_old(:,m,l-1,ifree,kkk,kk))
987              ELSEWHERE ((wat_flux(:,l,pref_soil_veg(m)) .LT. 0.) .AND. (wat_flux(:,l-1,pref_soil_veg(m)) .LT. 0.))
988                 DOC(:,m,l,ifree,kkk,kk) = DOC_old(:,m,l,ifree,kkk,kk) + MIN(DOC_FLUX(:,m,l,kkk,kk),DOC_old(:,m,l+1,ifree,kkk,kk)) - MIN(DOC_FLUX(:,m,l-1,kkk,kk),DOC_old(:,m,l,ifree,kkk,kk))
989              ELSEWHERE ((wat_flux(:,l,pref_soil_veg(m)) .GT. 0.) .AND. (wat_flux(:,l-1,pref_soil_veg(m)) .LT. 0.))
990                 DOC(:,m,l,ifree,kkk,kk) = DOC_old(:,m,l,ifree,kkk,kk) - MIN(DOC_FLUX(:,m,l,kkk,kk),DOC_old(:,m,l,ifree,kkk,kk)) - MIN(DOC_FLUX(:,m,l-1,kkk,kk),DOC_old(:,m,l,ifree,kkk,kk))
991              ELSEWHERE ((wat_flux(:,l,pref_soil_veg(m)) .LT. 0.) .AND. (wat_flux(:,l-1,pref_soil_veg(m)) .GT. 0.))
992                 DOC(:,m,l,ifree,kkk,kk) = DOC_old(:,m,l,ifree,kkk,kk) + MIN(DOC_FLUX(:,m,l,kkk,kk),DOC_old(:,m,l+1,ifree,kkk,kk)) + MIN(DOC_FLUX(:,m,l-1,kkk,kk),DOC_old(:,m,l-1,ifree,kkk,kk))
993              ELSEWHERE 
994                 DOC(:,m,l,ifree,kkk,kk) = DOC_old(:,m,l,ifree,kkk,kk)
995              ENDWHERE
996           ENDDO !End loop over soil and litter pool
997        ENDDO !End loop over elements
998     ENDDO !End loop over soil layers
999
1000     !For the last layer
1001     DO kk = 1,nelements !Loop over elements
1002        DO kkk = 1,npool !Loop over soil and litter pool
1003           WHERE (wat_flux(:,nbdl-1,pref_soil_veg(m)) .GT. 0.)
1004              DOC(:,m,nbdl,ifree,kkk,kk) = DOC_old(:,m,nbdl,ifree,kkk,kk) + MIN(DOC_FLUX(:,m,nbdl-1,kkk,kk),DOC_old(:,m,nbdl-1,ifree,kkk,kk))
1005           ELSEWHERE (wat_flux(:,nbdl-1,pref_soil_veg(m)) .LT. 0.)
1006              DOC(:,m,nbdl,ifree,kkk,kk) = DOC_old(:,m,nbdl,ifree,kkk,kk) - MIN(DOC_FLUX(:,m,nbdl-1,kkk,kk),DOC_old(:,m,nbdl,ifree,kkk,kk))
1007           ELSEWHERE
1008              DOC(:,m,nbdl,ifree,kkk,kk) = DOC_old(:,m,nbdl,ifree,kkk,kk)
1009           ENDWHERE
1010        ENDDO !End loop over soil and litter pool
1011     ENDDO !End loop over elements
1012
1013     !! 2.3.8 DOC transport by diffusion
1014     !! 2.3.8.1 Calculation of the DOC fluxes by diffusion
1015     ! Here we use the DOC concentration in the soil and not in water to
1016     ! stabilize the calcul of
1017     ! water distribution is continuous in the soil column (no dry places)
1018
1019
1020     DOC_old2(:,m,:,:,:,:)=zero
1021     DOC_old_buffer(:,m,:,:,:,:)=zero
1022     DOC_old2(:,m,:,:,:,:)=DOC(:,m,:,:,:,:)
1023
1024     DO kkk = 1, npool ! Loop over soil and litter pools
1025        DO kk = 1, nelements ! Loop over elements
1026               
1027           DOC_FLUX_DIFF(:,m,1,kkk,kk) = Dif_DOC(:)*ABS(DOC_old2(:,m,1,ifree,kkk,kk)/(z_soil(1))-DOC_old2(:,m,2,ifree,kkk,kk)/(z_soil(2)-z_soil(1)))/(z_soil(2)-z_soil(1))
1028
1029           DO l= 2, nbdl-1 ! Loop over soil layers
1030              DOC_FLUX_DIFF(:,m,l,kkk,kk) = Dif_DOC(:)*ABS(DOC_old2(:,m,l,ifree,kkk,kk)/(z_soil(l)-z_soil(l-1))-DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l)))/(z_soil(l+1)-z_soil(l))
1031           ENDDO
1032
1033       !Below we checked if in case that, in a given layer, you have diffusion
1034       !in the above and below litters, both fluxes are not higher than the
1035       !stocks of the given layer.
1036           DO l =1, nbdl-2
1037              WHERE (DOC_old2(:,m,l,ifree,kkk,kk)/(z_soil(l)-z_soil(l-1)) .LT. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l)) .AND. &
1038                   DOC_old2(:,m,l+2,ifree,kkk,kk)/(z_soil(l+2)-z_soil(l+1)) .LT. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l)) .AND. &
1039                   DOC_FLUX_DIFF(:,m,l,kkk,kk) + DOC_FLUX_DIFF(:,m,l+1,kkk,kk) .GT. DOC_old2(:,m,l+1,ifree,kkk,kk))
1040                 DOC_FLUX_DIFF(:,m,l,kkk,kk) = DOC_old2(:,m,l+1,ifree,kkk,kk)*(DOC_FLUX_DIFF(:,m,l,kkk,kk)/(DOC_FLUX_DIFF(:,m,l,kkk,kk)+DOC_FLUX_DIFF(:,m,l+1,kkk,kk)))
1041                 DOC_FLUX_DIFF(:,m,l+1,kkk,kk) = DOC_old2(:,m,l+1,ifree,kkk,kk)*(DOC_FLUX_DIFF(:,m,l+1,kkk,kk)/(DOC_FLUX_DIFF(:,m,l,kkk,kk)+DOC_FLUX_DIFF(:,m,l+1,kkk,kk)))
1042              ELSEWHERE (DOC_old2(:,m,l,ifree,kkk,kk)/(z_soil(l)-z_soil(l-1)) .GE. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l)) .AND. &
1043                   DOC_old2(:,m,l+2,ifree,kkk,kk)/(z_soil(l+2)-z_soil(l+1)) .LT. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l)) .AND. &
1044                   DOC_old2(:,m,l+1,ifree,kkk,kk) + DOC_FLUX_DIFF(:,m,l,kkk,kk) - DOC_FLUX_DIFF(:,m,l+1,kkk,kk) .LE. min_stomate)
1045                 DOC_FLUX_DIFF(:,m,l,kkk,kk) = DOC_FLUX_DIFF(:,m,l,kkk,kk)
1046                 DOC_FLUX_DIFF(:,m,l+1,kkk,kk) = DOC_old2(:,m,l+1,ifree,kkk,kk) + DOC_FLUX_DIFF(:,m,l,kkk,kk)
1047              ELSEWHERE (DOC_old2(:,m,l,ifree,kkk,kk)/(z_soil(l)-z_soil(l-1)) .LT. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l)) .AND. &
1048                   DOC_old2(:,m,l+2,ifree,kkk,kk)/(z_soil(l+2)-z_soil(l+1)) .GE. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l)) .AND. &
1049                   DOC_old2(:,m,l+1,ifree,kkk,kk) - DOC_FLUX_DIFF(:,m,l,kkk,kk) + DOC_FLUX_DIFF(:,m,l+1,kkk,kk) .LE. min_stomate)
1050                 DOC_FLUX_DIFF(:,m,l,kkk,kk) = DOC_old2(:,m,l+1,ifree,kkk,kk) + DOC_FLUX_DIFF(:,m,l+1,kkk,kk)
1051                 DOC_FLUX_DIFF(:,m,l+1,kkk,kk) = DOC_FLUX_DIFF(:,m,l+1,kkk,kk)
1052              ELSEWHERE
1053                 DOC_FLUX_DIFF(:,m,l,kkk,kk) = DOC_FLUX_DIFF(:,m,l,kkk,kk)
1054                 DOC_FLUX_DIFF(:,m,l+1,kkk,kk) = DOC_FLUX_DIFF(:,m,l+1,kkk,kk)
1055              ENDWHERE
1056           ENDDO
1057
1058           DO l= 1, nbdl-1 ! Loop over soil layers
1059              WHERE ((DOC_old2(:,m,l,ifree,kkk,kk)/(z_soil(l)-z_soil(l-1)) .LT. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l))) .AND. &
1060                   ((DOC_FLUX_DIFF(:,m,l,kkk,kk) - DOC_old2(:,m,l+1,ifree,kkk,kk)) .GE. min_stomate))
1061                 DOC(:,m,l,ifree,kkk,kk) = DOC_old2(:,m,l,ifree,kkk,kk) + DOC_old2(:,m,l+1,ifree,kkk,kk) + DOC_old_buffer(:,m,l,ifree,kkk,kk)
1062                 DOC_old_buffer(:,m,l+1,ifree,kkk,kk) = - DOC_old2(:,m,l+1,ifree,kkk,kk) 
1063                 DOC(:,m,l+1,ifree,kkk,kk) = zero
1064              ELSEWHERE ((DOC_old2(:,m,l,ifree,kkk,kk)/(z_soil(l)-z_soil(l-1)) .LT. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l))) .AND. &
1065                   ((DOC_old2(:,m,l+1,ifree,kkk,kk) - DOC_FLUX_DIFF(:,m,l,kkk,kk)) .GT. min_stomate))
1066                 DOC(:,m,l,ifree,kkk,kk) = DOC_old2(:,m,l,ifree,kkk,kk) + DOC_FLUX_DIFF(:,m,l,kkk,kk) + DOC_old_buffer(:,m,l,ifree,kkk,kk)
1067                 DOC_old_buffer(:,m,l+1,ifree,kkk,kk) =  - DOC_FLUX_DIFF(:,m,l,kkk,kk)
1068                 DOC(:,m,l+1,ifree,kkk,kk) = DOC_old2(:,m,l+1,ifree,kkk,kk) - DOC_FLUX_DIFF(:,m,l,kkk,kk)
1069              ELSEWHERE ((DOC_old2(:,m,l,ifree,kkk,kk)/(z_soil(l)-z_soil(l-1)) .GT. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l))) .AND. &
1070                   ((DOC_FLUX_DIFF(:,m,l,kkk,kk) - DOC_old2(:,m,l,ifree,kkk,kk)) .GE. min_stomate))
1071                 DOC(:,m,l,ifree,kkk,kk) = zero + DOC_old_buffer(:,m,l,ifree,kkk,kk)
1072                 DOC_old_buffer(:,m,l+1,ifree,kkk,kk) = DOC_old2(:,m,l,ifree,kkk,kk)
1073                 DOC(:,m,l+1,ifree,kkk,kk) = DOC_old2(:,m,l+1,ifree,kkk,kk) + DOC_old2(:,m,l,ifree,kkk,kk)
1074              ELSEWHERE ((DOC_old2(:,m,l,ifree,kkk,kk)/(z_soil(l)-z_soil(l-1)) .GT. DOC_old2(:,m,l+1,ifree,kkk,kk)/(z_soil(l+1)-z_soil(l))) .AND. &
1075                   ((DOC_old2(:,m,l,ifree,kkk,kk) - DOC_FLUX_DIFF(:,m,l,kkk,kk)) .GT. min_stomate))
1076                 DOC(:,m,l,ifree,kkk,kk) = DOC_old2(:,m,l,ifree,kkk,kk) - DOC_FLUX_DIFF(:,m,l,kkk,kk) + DOC_old_buffer(:,m,l,ifree,kkk,kk)
1077                 DOC_old_buffer(:,m,l+1,ifree,kkk,kk) = DOC_FLUX_DIFF(:,m,l,kkk,kk)
1078                 DOC(:,m,l+1,ifree,kkk,kk) = DOC_old2(:,m,l+1,ifree,kkk,kk) + DOC_FLUX_DIFF(:,m,l,kkk,kk)
1079              ELSEWHERE
1080                 DOC(:,m,l,ifree,kkk,kk) = DOC_old2(:,m,l,ifree,kkk,kk)
1081                 DOC(:,m,l+1,ifree,kkk,kk) = DOC_old2(:,m,l+1,ifree,kkk,kk)
1082              ENDWHERE
1083           ENDDO ! End loop over soil layers
1084        ENDDO ! End loop over elements
1085     ENDDO ! End loop over DOC pools
1086
1087        !! 2.3.9 DOC export
1088        !! Calculate total DOC going out of the system (the factor 1E-3 is
1089        !! applied to convert runoff, floodout, drainage from kg/m2 to m3/m2)
1090        !! We assume that runoff and floodout affects the first layer of soil
1091        !! and drainage only the last layer.
1092        !! Drainage and runoff are given by soil type by hydrol.f90 but not floodout (it is a weighed mean) explaining why
1093        !! the calculation is a bit different between DOC_RUN and DOC_FLOOD. Moreover, floodout is computed as evaporation-precipitation,
1094        !! so it represent an export to the rivers only if floodout is negative explaing why we have to multiply by -un.
1095        !! Be careful we do not * by dt because it is already done in hydrol.f90
1096
1097        !! based on fast reservoir, we calculate a reduction factor. The calculation is based on the assumption that there
1098        !! head waters that is the main source of DOC to inland waters, while areas more remote from the head waters have a neglibile
1099        !! effect on the DOC in streams. We assume that the export from soils to head waters is proportional to the
1100        !! square root of the fast reservoir storage. This is consitnet with empirical finding relating stream width to discharge. As we
1101        !! do not directly consider a stream surface area, we simply assign a reference fast reservoir for which we assume that
1102        !! all top soils export DOC into the head waters (which can include streamlets). Note that very high runoff appears where swamps
1103        !! are pesent.
1104
1105
1106        fastr_corr(:) = (fastr(:)**undemi)/(fastr_ref**undemi)
1107        WHERE (fastr_corr(:) .LT. zero)
1108           fastr_corr(:) = zero
1109        ELSEWHERE
1110           fastr_corr(:) = fastr_corr(:)
1111        ENDWHERE
1112
1113        !! Based on maximum export concentration, we calculate a correction factor
1114         soilDOC_corr(:) = zero
1115         soilDOC_31mm(:) = zero
1116         DO kk = 1,nelements !Loop over elements
1117            DO kkk = 1,npool !Loop over soil and litter pool
1118               soilDOC_31mm(:) = soilDOC_31mm + DOC(:,m,1,ifree,kkk,kk)
1119               IF (sro_bottom .GT. 1) THEN
1120                  DO l=2,sro_bottom
1121                     soilDOC_31mm(:) = soilDOC_31mm(:) + DOC(:,m,l,ifree,kkk,kk)
1122                  ENDDO !l=2,sro_bottom
1123               ENDIF !(sro_bottom .GT. 1)
1124            ENDDO !kkk = 1,npool
1125         ENDDO !kk = 1,nelements
1126         !
1127         WHERE (soilDOC_31mm(:)*flux_red(:)*fastr_corr(:)/soilwater_31mm(:,pref_soil_veg(m)) .GT. DOCexp_max)
1128            soilDOC_corr(:) = DOCexp_max/(soilDOC_31mm(:)*flux_red(:)*fastr_corr(:)/soilwater_31mm(:,pref_soil_veg(m)))
1129         ELSEWHERE
1130            soilDOC_corr(:) = un
1131         ENDWHERE
1132         soilDOC_corr(:) = soilDOC_corr(:) * flux_red(:) * fastr_corr(:)
1133         !
1134         DO kk = 1,nelements !Loop over elements
1135            DO kkk = 1,npool !Loop over soil and litter pool
1136               !DOC lost through runoff for the first layer
1137               WHERE (soilwater_31mm(:,pref_soil_veg(m)) .GT. zero)
1138                  DOC_RUN(:,m,kkk,kk) = DOC_RUN(:,m,kkk,kk) + MIN(DOC(:,m,1,ifree,kkk,kk),lat_exp*soilDOC_corr(:) &
1139                       &              * (runoff_per_soil(:,pref_soil_veg(m))*1E-3/soilwater_31mm(:,pref_soil_veg(m))) &
1140                       &              * DOC(:,m,1,ifree,kkk,kk))
1141               ELSEWHERE
1142                  DOC_RUN(:,m,kkk,kk) = DOC_RUN(:,m,kkk,kk)
1143               ENDWHERE
1144               !
1145               IF (kkk .EQ. imetabo) THEN
1146                  DOC_RUN(:,m,iact,kk) = DOC_RUN(:,m,iact,kk) + DOC_RUN(:,m,imetabo,kk) 
1147                  DOC_RUN(:,m,imetabo,kk) = zero
1148               ELSEIF (kkk .EQ. istrabo) THEN                 
1149                  DOC_RUN(:,m,iact,kk) = DOC_RUN(:,m,iact,kk) + DOC_RUN(:,m,istrabo,kk) * (un - lignin_struc_above(:,m))       
1150                  DOC_RUN(:,m,islo,kk) = DOC_RUN(:,m,islo,kk) + DOC_RUN(:,m,istrabo,kk) * lignin_struc_above(:,m) 
1151                  DOC_RUN(:,m,istrabo,kk) = zero
1152               ELSEIF (kkk .EQ. istrbel) THEN
1153                  DOC_RUN(:,m,iact,kk) = DOC_RUN(:,m,iact,kk) + DOC_RUN(:,m,istrbel,kk) * (un - lignin_struc_below(:,m,1))       
1154                  DOC_RUN(:,m,islo,kk) = DOC_RUN(:,m,islo,kk) + DOC_RUN(:,m,istrbel,kk) * lignin_struc_below(:,m,1) 
1155                  DOC_RUN(:,m,istrbel,kk) = zero 
1156               ELSEIF (kkk .EQ. imetbel) THEN             
1157                  DOC_RUN(:,m,iact,kk) = DOC_RUN(:,m,iact,kk) + DOC_RUN(:,m,imetbel,kk) 
1158                  DOC_RUN(:,m,imetbel,kk) = zero
1159               ELSE
1160                  DOC_RUN(:,m,kkk,kk) = DOC_RUN(:,m,kkk,kk)
1161               ENDIF
1162               !
1163               !And for the layers 2-? layers
1164               IF (sro_bottom .GT. 1) THEN
1165                  DO l=2,sro_bottom
1166                     WHERE (soilwater_31mm(:,pref_soil_veg(m)) .GT. zero)
1167                        DOC_RUN(:,m,kkk,kk) = DOC_RUN(:,m,kkk,kk) + MIN(DOC(:,m,l,ifree,kkk,kk), &
1168                             &      lat_exp * soilDOC_corr(:) * (runoff_per_soil(:,pref_soil_veg(m))*1E-3/soilwater_31mm(:,pref_soil_veg(m))) &
1169                             &    * DOC(:,m,l,ifree,kkk,kk))
1170                     ELSEWHERE
1171                        DOC_RUN(:,m,kkk,kk) = DOC_RUN(:,m,kkk,kk)
1172                     ENDWHERE
1173                     !
1174                     IF (kkk .EQ. imetabo) THEN
1175                        DOC_RUN(:,m,iact,kk) = DOC_RUN(:,m,iact,kk) + DOC_RUN(:,m,imetabo,kk) 
1176                        DOC_RUN(:,m,imetabo,kk) = zero
1177                     ELSEIF (kkk .EQ. istrabo) THEN                 
1178                        DOC_RUN(:,m,iact,kk) = DOC_RUN(:,m,iact,kk) + DOC_RUN(:,m,istrabo,kk) * (un - lignin_struc_above(:,m))       
1179                        DOC_RUN(:,m,islo,kk) = DOC_RUN(:,m,islo,kk) + DOC_RUN(:,m,istrabo,kk) * lignin_struc_above(:,m) 
1180                        DOC_RUN(:,m,istrabo,kk) = zero
1181                     ELSEIF (kkk .EQ. istrbel) THEN
1182                        DOC_RUN(:,m,iact,kk) = DOC_RUN(:,m,iact,kk) + DOC_RUN(:,m,istrbel,kk) * (un - lignin_struc_below(:,m,l))       
1183                        DOC_RUN(:,m,islo,kk) = DOC_RUN(:,m,islo,kk) + DOC_RUN(:,m,istrbel,kk) * lignin_struc_below(:,m,l) 
1184                        DOC_RUN(:,m,istrbel,kk) = zero 
1185                     ELSEIF (kkk .EQ. imetbel) THEN             
1186                        DOC_RUN(:,m,iact,kk) = DOC_RUN(:,m,iact,kk) + DOC_RUN(:,m,imetbel,kk) 
1187                        DOC_RUN(:,m,imetbel,kk) = zero
1188                     ELSE
1189                        DOC_RUN(:,m,kkk,kk) = DOC_RUN(:,m,kkk,kk)
1190                     ENDIF
1191                  ENDDO
1192               ENDIF !sro_bottom .GT. 1
1193                     !
1194               WHERE (soil_mc(:,nbdl,pref_soil_veg(m)) .GT. zero)
1195                  DOC_DRAIN(:,m,kkk,kk) = DOC_DRAIN(:,m,kkk,kk) + MIN(DOC(:,m,nbdl,ifree,kkk,kk),&
1196                       &                                         lat_exp*DOC(:,m,nbdl,ifree,kkk,kk)*drainage_per_soil(:,pref_soil_veg(m))*1E-3 &
1197                       &                                         /((z_soil(nbdl)-z_soil(nbdl-1)) * soil_mc(:,nbdl,pref_soil_veg(m))))
1198               ELSEWHERE
1199                  DOC_DRAIN(:,m,kkk,kk) = DOC_DRAIN(:,m,kkk,kk)
1200               ENDWHERE
1201               !
1202               IF (kkk .EQ. imetabo) THEN
1203                  DOC_DRAIN(:,m,iact,kk) = DOC_DRAIN(:,m,iact,kk) + DOC_DRAIN(:,m,imetabo,kk) 
1204                  DOC_DRAIN(:,m,imetabo,kk) = zero
1205               ELSEIF (kkk .EQ. istrabo) THEN                 
1206                  DOC_DRAIN(:,m,iact,kk) = DOC_DRAIN(:,m,iact,kk) + DOC_DRAIN(:,m,istrabo,kk) * (un - lignin_struc_above(:,m))       
1207                  DOC_DRAIN(:,m,islo,kk) = DOC_DRAIN(:,m,islo,kk) + DOC_DRAIN(:,m,istrabo,kk) * lignin_struc_above(:,m) 
1208                  DOC_DRAIN(:,m,istrabo,kk) = zero
1209               ELSEIF (kkk .EQ. istrbel) THEN
1210                  DOC_DRAIN(:,m,iact,kk) = DOC_DRAIN(:,m,iact,kk) + DOC_DRAIN(:,m,istrbel,kk) * (un - lignin_struc_below(:,m,nbdl))       
1211                  DOC_DRAIN(:,m,islo,kk) = DOC_DRAIN(:,m,islo,kk) + DOC_DRAIN(:,m,istrbel,kk) * lignin_struc_below(:,m,nbdl) 
1212                  DOC_DRAIN(:,m,istrbel,kk) = zero 
1213               ELSEIF (kkk .EQ. imetbel) THEN             
1214                  DOC_DRAIN(:,m,iact,kk) = DOC_DRAIN(:,m,iact,kk) + DOC_DRAIN(:,m,imetbel,kk) 
1215                  DOC_DRAIN(:,m,imetbel,kk) = zero
1216               ELSE
1217                  DOC_DRAIN(:,m,kkk,kk) = DOC_DRAIN(:,m,kkk,kk)
1218               ENDIF
1219               !
1220            ENDDO !End loop over soil and litter pool         
1221            !! TF-DOC
1222            WHERE (veget_max(:,m) .GT. zero)
1223               DOC_FLOOD(:,m,iact,kk) = DOC_FLOOD(:,m,iact,kk) + wet_dep_flood(:,m,kk)*undemi/veget_max(:,m)
1224               DOC_FLOOD(:,m,islo,kk) = DOC_FLOOD(:,m,islo,kk) + wet_dep_flood(:,m,kk)*undemi/veget_max(:,m)
1225               DOC_FLOOD(:,m,iact,kk) = DOC_FLOOD(:,m,iact,kk) + floodcarbon_input(:,m,iact,kk)*dt
1226               DOC_FLOOD(:,m,islo,kk) = DOC_FLOOD(:,m,islo,kk) + floodcarbon_input(:,m,islo,kk)*dt
1227            ELSEWHERE
1228               DOC_FLOOD(:,m,iact,kk) = DOC_FLOOD(:,m,iact,kk)
1229            ENDWHERE
1230            !
1231            DO l=1,sro_bottom
1232               WHERE (veget_max(:,m) .GT. zero)
1233                  DOC_FLOOD(:,m,islo,kk) = DOC_FLOOD(:,m,islo,kk) + CUE_coef(:) * fluxtot_flood(:,islow,kk,l)
1234                  DOC_FLOOD(:,m,ipas,kk) = DOC_FLOOD(:,m,ipas,kk) + CUE_coef(:) * fluxtot_flood(:,ipassive,kk,l)
1235                  DOC_FLOOD(:,m,iact,kk) = DOC_FLOOD(:,m,iact,kk) + CUE_coef(:) * fluxtot_flood(:,iactive,kk,l)
1236               ELSEWHERE
1237                  DOC_FLOOD(:,m,iact,kk) = DOC_FLOOD(:,m,iact,kk)
1238               ENDWHERE
1239            ENDDO !l=1,sro_bottom
1240            !
1241         ENDDO !End loop over elements
1242
1243
1244         DO kk = 1,nelements !Loop over elements
1245            DOC_EXP(:,m,irunoff,:,kk) = DOC_RUN(:,m,:,kk)/dt
1246            DOC_EXP(:,m,iflooded,:,kk) = DOC_FLOOD(:,m,:,kk)/dt
1247            DOC_EXP(:,m,idrainage,:,kk) = DOC_DRAIN(:,m,:,kk)/dt
1248         ENDDO !End loop over elements
1249         ! We substract runoff, floodout and drainage from DOC pools
1250         DO kkk = 1,npool !Loop over soil and litter pool
1251            DO kk = 1, nelements
1252               
1253               !For the first layer
1254               WHERE (soilwater_31mm(:,pref_soil_veg(m)) .GT. zero)
1255                  DOC(:,m,1,ifree,kkk,kk) = DOC(:,m,1,ifree,kkk,kk) - MIN(DOC(:,m,1,ifree,kkk,kk), &
1256                       &                    lat_exp * soilDOC_corr(:) * (runoff_per_soil(:,pref_soil_veg(m))*1E-3/soilwater_31mm(:,pref_soil_veg(m))) &
1257                       &                  * DOC(:,m,1,ifree,kkk,kk))
1258               ELSEWHERE
1259                  DOC(:,m,1,ifree,kkk,kk) = DOC(:,m,1,ifree,kkk,kk)
1260               ENDWHERE
1261               !And for the layers 2-5
1262               IF (sro_bottom .GT. 1) THEN
1263                  DO l=2,sro_bottom
1264                     WHERE (soilwater_31mm(:,pref_soil_veg(m)) .GT. zero)
1265                        DOC(:,m,l,ifree,kkk,kk) = DOC(:,m,l,ifree,kkk,kk) - MIN(DOC(:,m,l,ifree,kkk,kk), &
1266                             &                    lat_exp * soilDOC_corr(:) * (runoff_per_soil(:,pref_soil_veg(m))*1E-3/soilwater_31mm(:,pref_soil_veg(m))) &
1267                             &                  * DOC(:,m,l,ifree,kkk,kk))
1268                     ELSEWHERE
1269                        DOC(:,m,l,ifree,kkk,kk) = DOC(:,m,l,ifree,kkk,kk)
1270                     ENDWHERE
1271                  ENDDO
1272               ENDIF !sro_bottom .GT. 1
1273               !For the last layer
1274               WHERE (soil_mc(:,nbdl,pref_soil_veg(m)) .GT. zero)
1275                  DOC(:,m,nbdl,ifree,kkk,kk) = DOC(:,m,nbdl,ifree,kkk,kk) - &
1276                       &                       MIN(DOC(:,m,nbdl,ifree,kkk,kk),&
1277                       &                       lat_exp * DOC(:,m,nbdl,ifree,kkk,kk)*drainage_per_soil(:,pref_soil_veg(m))*1E-3 &
1278                       &                       /((z_soil(nbdl)-z_soil(nbdl-1)) * soil_mc(:,nbdl,pref_soil_veg(m))))
1279               ELSEWHERE
1280                  DOC(:,m,nbdl,ifree,kkk,kk) = DOC(:,m,nbdl,ifree,kkk,kk)
1281               ENDWHERE
1282            ENDDO ! End loop over elements
1283         ENDDO !End loop over soil and litter pool
1284
1285      !! 2.3.10 Transport of particulate organic C (i.e. active, slow and passive
1286      !! pools) following the second Fick's Law (O¿Brien and Stout, 1978; Wynn
1287      !! et al., 2005). Represent the bioturbation and is generally associated
1288      !! with transport  plant debris (Elzein and Balesdent, 1995; Bruun et al.,
1289      !! 2007; Braakhekke et al., 2011, Guenet et al., 2013)
1290
1291      ! Definition of the diffusion rate from Guenet et al., (2013) BG for the
1292      ! moment but must depends on clay, ph organic C, or whatever...
1293     
1294
1295carbon_old(:,:,m,:)=zero
1296carbon_old_buffer(:,:,m,:)=zero
1297carbon_old(:,:,m,:)=carbon(:,:,m,:)
1298
1299     DO kk = 1, ncarb ! Loop over soil and litter pools
1300
1301        carbon_flux(:,kk,m,1) = Dif_coef(:)*ABS(carbon_old(:,kk,m,1)/(z_soil(1))-carbon_old(:,kk,m,2)/(z_soil(2)-z_soil(1)))/(z_soil(2)-z_soil(1))
1302
1303      DO l= 2, nbdl-1 ! Loop over soil layers
1304        carbon_flux(:,kk,m,l) = Dif_coef(:)*ABS(carbon_old(:,kk,m,l)/(z_soil(l)-z_soil(l-1))-carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l)))/(z_soil(l+1)-z_soil(l))
1305      ENDDO
1306        carbon_flux_old(:,kk,m,:) = carbon_flux(:,kk,m,:)
1307
1308       !Below we checked if in case that, in a given layer, you have diffusion
1309       !in the above and below litters, both fluxes are not higher than the
1310       !stocks of the given layer.
1311
1312      DO l =1, nbdl-2
1313           WHERE (carbon_old(:,kk,m,l)/(z_soil(l)-z_soil(l-1)) .LT. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l)) .AND. &
1314                  carbon_old(:,kk,m,l+2)/(z_soil(l+2)-z_soil(l+1)) .LT. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l)) .AND. &
1315                  carbon_flux(:,kk,m,l)+carbon_flux(:,kk,m,l+1) .GT. carbon_old(:,kk,m,l+1))
1316            carbon_flux(:,kk,m,l) = carbon_old(:,kk,m,l+1)*(carbon_flux(:,kk,m,l)/(carbon_flux(:,kk,m,l)+carbon_flux(:,kk,m,l+1)))
1317            carbon_flux(:,kk,m,l+1) = carbon_old(:,kk,m,l+1)*(carbon_flux(:,kk,m,l+1)/(carbon_flux(:,kk,m,l)+carbon_flux(:,kk,m,l+1)))
1318           ELSEWHERE (carbon_old(:,kk,m,l)/(z_soil(l)-z_soil(l-1)) .GE. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l)) .AND. &
1319                      carbon_old(:,kk,m,l+2)/(z_soil(l+2)-z_soil(l+1)) .LT. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l)) .AND. &
1320                      carbon_old(:,kk,m,l+1) + carbon_flux(:,kk,m,l) - carbon_flux(:,kk,m,l+1) .LE. min_stomate)
1321            carbon_flux(:,kk,m,l) = carbon_flux(:,kk,m,l)
1322            carbon_flux(:,kk,m,l+1) = carbon_old(:,kk,m,l+1) + carbon_flux(:,kk,m,l)
1323           ELSEWHERE (carbon_old(:,kk,m,l)/(z_soil(l)-z_soil(l-1)) .LT. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l)) .AND. &
1324                      carbon_old(:,kk,m,l+2)/(z_soil(l+2)-z_soil(l+1)) .GE. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l)) .AND. &
1325                      carbon_old(:,kk,m,l+1) - carbon_flux(:,kk,m,l) + carbon_flux(:,kk,m,l+1) .LE. min_stomate)
1326            carbon_flux(:,kk,m,l) = carbon_old(:,kk,m,l+1) + carbon_flux(:,kk,m,l+1)
1327            carbon_flux(:,kk,m,l+1) = carbon_flux(:,kk,m,l+1)
1328           ELSEWHERE
1329            carbon_flux(:,kk,m,l) = carbon_flux(:,kk,m,l)
1330            carbon_flux(:,kk,m,l+1) = carbon_flux(:,kk,m,l+1)
1331           ENDWHERE
1332      ENDDO
1333
1334      DO l= 1, nbdl-1 ! Loop over soil layers
1335           WHERE ((carbon_old(:,kk,m,l)/(z_soil(l)-z_soil(l-1)) .LT. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l))) .AND. &
1336                 ((carbon_flux(:,kk,m,l) - carbon_old(:,kk,m,l+1)) .GE. min_stomate))
1337            carbon(:,kk,m,l) = carbon_old(:,kk,m,l) + carbon_old(:,kk,m,l+1) + carbon_old_buffer(:,kk,m,l)
1338            carbon_old_buffer(:,kk,m,l+1) = -carbon_old(:,kk,m,l+1) 
1339            carbon(:,kk,m,l+1) = zero
1340           ELSEWHERE ((carbon_old(:,kk,m,l)/(z_soil(l)-z_soil(l-1)) .LT. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l))) .AND. &
1341                  ((carbon_old(:,kk,m,l+1) - carbon_flux(:,kk,m,l)) .GT. min_stomate))
1342            carbon(:,kk,m,l) = carbon_old(:,kk,m,l) + carbon_flux(:,kk,m,l) + carbon_old_buffer(:,kk,m,l)
1343            carbon_old_buffer(:,kk,m,l+1) =  - carbon_flux(:,kk,m,l)
1344            carbon(:,kk,m,l+1) = carbon_old(:,kk,m,l+1) - carbon_flux(:,kk,m,l)
1345           ELSEWHERE ((carbon_old(:,kk,m,l)/(z_soil(l)-z_soil(l-1)) .GT. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l))) .AND. &
1346                  ((carbon_flux(:,kk,m,l) - carbon_old(:,kk,m,l)) .GE. min_stomate))
1347            carbon(:,kk,m,l) =zero + carbon_old_buffer(:,kk,m,l)
1348            carbon_old_buffer(:,kk,m,l+1) = carbon_old(:,kk,m,l)
1349            carbon(:,kk,m,l+1) = carbon_old(:,kk,m,l+1) + carbon_old(:,kk,m,l)
1350           ELSEWHERE ((carbon_old(:,kk,m,l)/(z_soil(l)-z_soil(l-1)) .GT. carbon_old(:,kk,m,l+1)/(z_soil(l+1)-z_soil(l))) .AND. &
1351                  ((carbon_old(:,kk,m,l) - carbon_flux(:,kk,m,l)) .GT. min_stomate))
1352            carbon(:,kk,m,l) = carbon_old(:,kk,m,l) - carbon_flux(:,kk,m,l) + carbon_old_buffer(:,kk,m,l)
1353            carbon_old_buffer(:,kk,m,l+1) = carbon_flux(:,kk,m,l)
1354            carbon(:,kk,m,l+1) = carbon_old(:,kk,m,l+1) + carbon_flux(:,kk,m,l)
1355           ELSEWHERE
1356            carbon(:,kk,m,l) = carbon_old(:,kk,m,l)
1357            carbon(:,kk,m,l+1) = carbon_old(:,kk,m,l+1)
1358           ENDWHERE
1359        ENDDO ! End loop over soil layers
1360
1361
1362
1363      ENDDO ! End loop over carbon pools
1364  ENDDO ! End loop over PFTs
1365
1366!! 3. Check mass balance closure
1367IF (ld_doc) THEN   
1368    !! 3.1 Calculate components of the mass balance
1369    pool_end(:,:,:) = zero 
1370
1371!Ending carbon pool
1372    DO m = 1, nvm
1373       DO i = 1,ncarb
1374          DO l = 1, nbdl
1375                 pool_end(:,m,icarbon) = pool_end(:,m,icarbon) + &
1376                           ( carbon(:,i,m,l)* veget_max(:,m) ) 
1377          ENDDO
1378       ENDDO
1379    ENDDO
1380
1381!Ending DOC pool
1382    DO m = 1, nvm
1383       DO i = 1,npool
1384          DO l = 1, nbdl
1385             DO j = 1, ndoc
1386                    pool_end(:,m,icarbon) = pool_end(:,m,icarbon) + &
1387                              ( DOC(:,m,l,j,i,icarbon) * veget_max(:,m)) 
1388             ENDDO
1389          ENDDO
1390       ENDDO
1391       pool_end(:,m,icarbon) = pool_end(:,m,icarbon) +  interception_storage(:,m,icarbon)
1392    ENDDO
1393
1394!WRITE(numout,*) " pool_start(:,m,icarbon) soil=",  pool_start(:,4,icarbon)
1395!WRITE(numout,*) " pool_end(:,m,icarbon) soil=",  pool_end(:,4,icarbon)
1396
1397 !! 3.2 Calculate mass balance
1398    !  Note that soilcarbon is transfered to other pools but that the pool
1399    !  was not updated. We should not account for it in ::pool_end
1400    check_intern(:,:,:,:) = zero
1401    check_intern(:,:,iatm2land,icarbon) = zero
1402
1403   
1404    check_intern(:,:,iland2atm,icarbon) = -un * (resp_hetero_soil(:,:)+resp_flood_soil(:,:)) * &
1405         veget_max(:,:) * dt 
1406    DO i=1,npool
1407       check_intern(:,:,ilat2out,icarbon) = check_intern(:,:,ilat2out,icarbon) - un * (DOC_EXP(:,:,irunoff,i,icarbon)+ &
1408            &         DOC_EXP(:,:,iflooded,i,icarbon)+ DOC_EXP(:,:,idrainage,i,icarbon)) * veget_max(:,:) * dt
1409    ENDDO !i=1,npool
1410    check_intern(:,:,ilat2in,icarbon) = un * zero
1411    check_intern(:,:,ipoolchange,icarbon) = -un * (pool_end(:,:,icarbon) - &
1412         pool_start(:,:,icarbon))
1413   
1414
1415    closure_intern = zero
1416
1417    DO i = 1,nmbcomp
1418       closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + &
1419           check_intern(:,:,i,icarbon)
1420    ENDDO
1421
1422!! 3.3 Write verdict
1423    DO i = 1,npts
1424       IF (SUM(closure_intern(i,:,icarbon)) .GT. min_stomate .OR. &
1425            SUM(closure_intern(i,:,icarbon)) .LT. -min_stomate) THEN
1426          WRITE(numout,*) 'Error: mass balance is not closed in stomate_soilcarbon.f90 at pixel ', i
1427          WRITE(numout,*) '   Difference is, ', SUM(closure_intern(i,:,icarbon))
1428          WRITE(numout,*) '   Vegetation is, ', veget_max(i,:)
1429
1430          WRITE(numout,*) '   deposition, ', SUM(wet_dep_ground(i,:,icarbon))
1431          WRITE(numout,*) '   deposition flood, ', SUM(wet_dep_flood(i,:,icarbon))
1432
1433          WRITE(numout,*) '   soilcarbon_input, ', SUM(SUM(soilcarbon_input(i,:,:,:,icarbon),DIM=2),DIM=2)
1434
1435          WRITE(numout,*) '   DOC_EXP, ',  SUM(SUM(DOC_EXP(i,:,:,:,icarbon),DIM=2),DIM=2)
1436          WRITE(numout,*) '   FLOOD_carbon_input, ', SUM(floodcarbon_input(i,:,:,icarbon),DIM=2)
1437
1438          WRITE(numout,*) '   resp_hetero_soil(:,:), ', SUM(resp_hetero_soil(i,:))
1439          WRITE(numout,*) '   resp_flood_soil(:,:), ', SUM(resp_flood_soil(i,:))
1440       ENDIF
1441    ENDDO !i = 1,npts
1442ENDIF
1443
1444!! 4 Writing history files
1445
1446    CALL xios_orchidee_send_field("dry_dep_canopy",dry_dep_canopy(:,:,icarbon)/dt)
1447    CALL xios_orchidee_send_field("DOC_precip2canopy",DOC_precip2canopy(:,:,icarbon)/dt)
1448    CALL xios_orchidee_send_field("DOC_precip2ground",DOC_precip2ground(:,:,icarbon)/dt)
1449    CALL xios_orchidee_send_field("DOC_canopy2ground",DOC_canopy2ground(:,:,icarbon)/dt)
1450!    CALL xios_orchidee_send_field("EXP_DOC_RUNOFF_ref",DOC_EXP(:,:,irunoff,islo,kk)+DOC_EXP(:,:,irunoff,ipas,kk))
1451!    CALL xios_orchidee_send_field("EXP_DOC_RUNOFF_lab",DOC_EXP(:,:,irunoff,iact,kk))
1452!    CALL xios_orchidee_send_field("EXP_DOC_DRAIN_ref",DOC_EXP(:,:,idrainage,islo,kk)+DOC_EXP(:,:,idrainage,ipas,kk))
1453!    CALL xios_orchidee_send_field("EXP_DOC_DRAIN_lab",DOC_EXP(:,:,idrainage,iact,kk))
1454!    CALL xios_orchidee_send_field("EXP_DOC_FLOOD_ref",DOC_EXP(:,:,iflooded,islo,kk)+DOC_EXP(:,:,iflooded,ipas,kk))
1455!    CALL xios_orchidee_send_field("EXP_DOC_FLOOD_lab",DOC_EXP(:,:,iflooded,iact,kk))
1456    CALL xios_orchidee_send_field("interception_storage",interception_storage(:,:,icarbon))
1457    CALL xios_orchidee_send_field("precip2canopy",precip2canopy(:,:)/dt)
1458    CALL xios_orchidee_send_field("precip2ground",precip2ground(:,:)/dt)
1459    CALL xios_orchidee_send_field("canopy2ground",canopy2ground(:,:)/dt)
1460    CALL xios_orchidee_send_field("runoff_forest",runoff_per_soil(:,2)/dt)
1461    CALL xios_orchidee_send_field("CLAY",clay(:))
1462    CALL xios_orchidee_send_field("SOIL_pH",soil_pH(:))
1463    CALL xios_orchidee_send_field("KD",Kd(:,:))
1464    CALL xios_orchidee_send_field("POOR_SOILS",poor_soils(:))
1465    IF (printlev>=4) WRITE(numout,*) 'Leaving soilcarbon'
1466   
1467  END SUBROUTINE soilcarbon
1468
1469END MODULE stomate_soilcarbon
Note: See TracBrowser for help on using the repository browser.