Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_stomate/stomate_litter.f90 @ 7852

Last change on this file since 7852 was 7615, checked in by sebastiaan.luyssaert, 2 years ago
Enhanced consistency of variable names: input has been changed in n_input throughout the code and the variable name vegstress introduced in sechiba is now also used in stomate. Enhnaced computational consistency: Pgap_cumul is used in stomate rather than recalculating it before calculating light_tran_to_floor_season. Edited getin_p while checking the code (but no real changes were made) and added several missing stomate and sechiba variables to age_class_distr to improve the 1+1=2 issue when comparing a model run with against a run without age classes. Finally: Write warning 10b in allocation to the history file rather than the out_orchidee file
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 137.9 KB

Line
1	! ================================================================================================================================
2	! MODULE : stomate_litter
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Update litter and lignine content after litter fall and
10	!! calculating litter decomposition.
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL$
20	!! $Date$
21	!! $Revision$
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE stomate_litter
26
27	! modules used:
28
29	USE ioipsl_para
30	USE stomate_data
31	USE constantes
32	USE constantes_soil
33	USE pft_parameters
34	USE function_library, ONLY : check_mass_balance, biomass_to_lai, cc_to_biomass
35	USE grid, ONLY : area
36
37	IMPLICIT NONE
38
39	! private & public routines
40
41	PRIVATE
42	PUBLIC littercalc,littercalc_clear, deadleaf
43
44	LOGICAL, SAVE :: firstcall_litter = .TRUE. !! first call
45	!$OMP THREADPRIVATE(firstcall_litter)
46	INTEGER(i_std), SAVE :: printlev_loc !! Local level of text output for current module
47	!$OMP THREADPRIVATE(printlev_loc)
48
49	CONTAINS
50
51	!! ================================================================================================================================
52	!! SUBROUTINE : littercalc_clear
53	!!
54	!!\BRIEF Set the flag ::firstcall_litter to .TRUE. and as such activate section
55	!! 1.1 of the subroutine littercalc (see below).
56	!!
57	!! DESCRIPTION : None
58	!!
59	!! RECENT CHANGE(S) : None
60	!!
61	!! MAIN OUTPUT VARIABLE(S) : None
62	!!
63	!! REFERENCE(S) : None
64	!!
65	!! FLOWCHART : None
66	!! \n
67	!_ ================================================================================================================================
68
69	SUBROUTINE littercalc_clear
70	firstcall_litter =.TRUE.
71	END SUBROUTINE littercalc_clear
72
73
74	!! ================================================================================================================================
75	!! SUBROUTINE : littercalc
76	!!
77	!!\BRIEF Calculation of the litter decomposition and therefore of the
78	!! heterotrophic respiration from litter following Parton et al. (1987).
79	!!
80	!! DESCRIPTION : The littercal routine splits the litter in 4 pools:
81	!! aboveground metaboblic, aboveground structural, belowground metabolic and
82	!! belowground structural. the fraction (F) of plant material going to metabolic
83	!! and structural is defined following Parton et al. (1987)
84	!! \latexonly
85	!! \input{littercalc1.tex}
86	!! \endlatexonly
87	!! \n
88	!! where L is the lignin content of the plant carbon pools considered and CN
89	!! its CN ratio. L and CN are fixed parameters for each plant carbon pools,
90	!! therefore it is the ratio between each plant carbon pool within a PFT, which
91	!! controlled the part of the total litter, that will be considered as
92	!! recalcitrant (i.e. structural litter) or labile (i.e. metabolic litter).\n
93	!!
94	!! The routine calculates the fraction of aboveground litter which is metabolic
95	!! or structural (the litterpart variable) which is then used in lpj_fire.f90.\n
96	!!
97	!! In the section 2, the routine calculate the new plant material entering the
98	!! litter pools by phenological death of plants organs (corresponding to the
99	!! variable turnover) and by fire, herbivory and others non phenological causes
100	!! (variable bm_to_litter). This calculation is first done for each PFT and then
101	!! the values calculated for each PFT are added up. Following the same approach
102	!! the lignin content of the total structural litter is calculated and will be
103	!! then used as a factor control of the decomposition of the structural litter
104	!! (lignin_struc) in the section 5.1.2. A test is performed to avoid that we add
105	!! more lignin than structural litter. Finally, the variable litterpart is
106	!! updated.\n
107	!!
108	!! In the section 3 and 4 the temperature and the moisture controlling the
109	!! decomposition are calculated for above and belowground. For aboveground
110	!! litter, air temperature and litter moisture are calculated in sechiba and used
111	!! directly. For belowground, soil temperature and moisture are also calculated
112	!! in sechiba but are modulated as a function of the soil depth. The modulation
113	!! is a multiplying factor exponentially distributed between 0 (in depth) and 1
114	!! in surface.\n
115	!!
116	!! Then, in the section 5, the routine calculates the structural litter decomposition
117	!! (C) following first order kinetics following Parton et al. (1987).
118	!! \latexonly
119	!! \input{littercalc2.tex}
120	!! \endlatexonly
121	!! \n
122	!! with k the decomposition rate of the structural litter.
123	!! k corresponds to
124	!! \latexonly
125	!! \input{littercalc3.tex}
126	!! \endlatexonly
127	!! \n
128	!! with littertau the turnover rate, T a function of the temperature and M a function of
129	!! the moisture described below.\n
130	!!
131	!! Then, the fraction of dead leaves (DL) composed by aboveground structural litter is
132	!! calculated as following
133	!! \latexonly
134	!! \input{littercalc4.tex}
135	!! \endlatexonly
136	!! \n
137	!! with k the decomposition rate of the structural litter previously
138	!! described.\n
139	!!
140	!! In the section 5.1, the fraction of decomposed structural litter
141	!! incorporated to the soil (Input) and its associated heterotrophic respiration are
142	!! calculated. For structural litter, the C decomposed could go in the active
143	!! soil carbon pool or in the slow carbon, as described in
144	!! stomate_soilcarbon.f90.\n
145	!! \latexonly
146	!! \input{littercalc5.tex}
147	!! \endlatexonly
148	!! \n
149	!! with f a parameter describing the fraction of structural litter incorporated
150	!! into the considered soil carbon pool, C the amount of litter decomposed and L
151	!! the amount of lignin in the litter. The litter decomposed which is not
152	!! incorporated into the soil is respired.\n
153	!!
154	!! In the section 5.2, the fraction of decomposed metabolic litter
155	!! incorporated to the soil and its associated heterotrophic respiration are
156	!! calculated with the same approaches presented for 5.1 but no control factor
157	!! depending on the lignin content are used.\n
158	!!
159	!! In the section 6 the dead leaf cover is calculated through a call to the
160	!! deadleaf subroutine presented below.\n
161	!!
162	!! In the section 7, if the flag SPINUP_ANALYTIC is set to true, we fill MatrixA
163	!! and VectorB following Lardy(2011).
164	!!
165	!! MAIN OUTPUT VARIABLES: ::deadleaf_cover, ::resp_hetero_litter
166	!! ::control_temp, ::control_moist
167	!!
168	!! REFERENCES:
169	!! - Parton, WJ, Schimel, DS, Cole, CV, and Ojima, DS. 1987. Analysis
170	!! of factors controlling soil organic matter levels in Great Plains
171	!! grasslands. Soil Science Society of America journal (USA)
172	!! (51):1173-1179.
173	!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
174	!! Environmental Modelling & Software (2011), doi:10.1016\|j.envsoft.2011.05.016
175	!!
176	!! FLOWCHART :
177	!! \latexonly
178	!! \includegraphics(scale=0.5){littercalcflow.jpg}
179	!! \endlatexonly
180	!! \n
181	!_ ================================================================================================================================
182
183	SUBROUTINE littercalc (npts, turnover, bm_to_litter, tree_bm_to_litter, &
184	veget_max, tsurf, stempdiag, shumdiag, litterhum, som, clay, silt, &
185	soil_n_min, n_input, month, harvest_pool_acc, litter, dead_leaves, lignin_struc, &
186	lignin_wood, n_mineralisation, deadleaf_cover, resp_hetero_litter, &
187	som_input, control_temp, control_moist, n_fungivores, &
188	MatrixA, VectorB, CN_target, CN_som_litter_longterm, tau_CN_longterm, &
189	ld_redistribute, circ_class_biomass, circ_class_n, tsoil_decomp)
190
191	!! 0. Variable and parameter declaration
192
193	!! 0.1 Input variables
194
195	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels
196	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: turnover !! Turnover rates of plant biomass
197	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
198	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: bm_to_litter !! Conversion of biomass to litter
199	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
200	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: tree_bm_to_litter !! Conversion of biomass to litter
201	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
202	REAL(r_std),DIMENSION(:,:),INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
203	!! defined in the input vegetation map
204	!! @tex $(m^2 m^{-2})$ @endtex
205	REAL(r_std), DIMENSION(:), INTENT(in) :: tsurf !! Temperature (K) at the surface
206	REAL(r_std), DIMENSION(:,:), INTENT(in) :: stempdiag !! Soil temperature (K)
207	REAL(r_std), DIMENSION(:,:), INTENT(in) :: shumdiag !! Daily soil humidity of each soil layer
208	!! (unitless)
209	REAL(r_std), DIMENSION(:), INTENT(in) :: litterhum !! Daily litter humidity (unitless)
210	REAL(r_std), DIMENSION(:), INTENT(in) :: clay !! Clay fraction (unitless, 0-1)
211	REAL(r_std), DIMENSION(:), INTENT(in) :: silt !! Silt fraction (unitless, 0-1)
212	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass components of the model tree
213	!! within a circumference class
214	!! @tex $(g C ind^{=1})$ @endtex
215	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of individuals in each circ class
216	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: n_input !! nitrogen inputs into the soil (gN/m**2/day)
217	INTEGER(i_std), INTENT(in) :: month !! month number required for n_input (1-12)
218
219
220	!! 0.2 Output variables
221
222	REAL(r_std), DIMENSION(:), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves
223	!! over all PFTs (0-1, unitless)
224	REAL(r_std), DIMENSION(:,:), INTENT(out) :: resp_hetero_litter !! Litter heterotrophic respiration. The unit
225	!! is given by m^2 of ground.
226	!! @tex $(gC dt_sechiba one\_day^{-1}) m^{-2})$ @endtex
227	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: som_input !! Quantity of Carbon (or Nitrogen) going into SOM pools
228	!! from litter decomposition. The unit is
229	!! given by m^2 of ground
230	!! @tex $(gC(orN) m^{-2} dt\_slow^{-1})$ @endtex
231	REAL(r_std), DIMENSION(:,:), INTENT(out) :: control_temp !! Temperature control of heterotrophic
232	!! respiration, above and below (0-1,
233	!! unitless)
234	REAL(r_std), DIMENSION(:,:), INTENT(out) :: control_moist !! Moisture control of heterotrophic
235	!! respiration (0.25-1, unitless)
236	REAL(r_std), DIMENSION(:,:), INTENT(out) :: n_fungivores !! Fraction of N released for plant uptake
237	!! due to fungivore consumption.
238	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: MatrixA !! Matrix containing the fluxes between the
239	!! carbon pools per sechiba time step
240	!! @tex $(gC.m^2.day^{-1})$ @endtex
241	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: VectorB !! Vector containing the litter increase per
242	!! sechiba time step
243	!! @tex $(gC m^{-2})$ @endtex
244	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: CN_target !! C to N ratio of SOM flux from one pool to another (gN m-2 dt-1)
245	REAL(r_std), DIMENSION(:), INTENT(out) :: ld_redistribute !! logical set to redistribute som and litter
246	REAL(r_std), DIMENSION(:), INTENT(out) :: tsoil_decomp !! Temperature used for decompostition in soil (K)
247
248
249	!! 0.3 Modified variables
250
251	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: som !! SOM pools: active, slow, or passive, \f$(gC m^{2})$\f
252	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! Metabolic and structural litter,above and
253	!! below ground. The unit is given by m^2 of
254	!! ground @tex $(gC m^{-2})$ @endtex
255	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
256	!! metabolic and structural in
257	!! @tex $(gC m^{-2})$ @endtex
258	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! Ratio Lignin content in structural litter,
259	!! above and below ground, (0-1, unitless)
260	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_wood !! Ratio Lignin/Carbon in woody litter,
261	!! above and below ground, (0-1, unitless)
262	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: n_mineralisation !! Mineral N pool (gN m-2)
263	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: CN_som_litter_longterm !! Longterm CN ratio of litter and som pools (gC/gN)
264	REAL(r_std), INTENT(inout) :: tau_CN_longterm !! Counter used for calculating the longterm CN_ratio of som and litter pools [seconds]
265	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: harvest_pool_acc !! The wood and biomass havested by humans
266	REAL(r_std), DIMENSION(:,:,:),INTENT(inout) :: soil_n_min !! mineral nitrogen in the soil (gN/m**2)
267
268	!! 0.4 Local variables
269
270	REAL(r_std) :: dt !! Number of sechiba(fast processes) time-step per day
271	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: litterfrac !! The fraction of leaves, wood, etc. that
272	!! goes into metabolic and structural
273	!! litterpools (0-1, unitless)
274	!$OMP THREADPRIVATE(litterfrac)
275	REAL(r_std) :: rpc !! Integration constant for vertical root
276	!! profiles (unitless)
277	REAL(r_std), SAVE, DIMENSION(nlitt) :: litter_dec_fac !! Turnover time in litter pools (days)
278	!$OMP THREADPRIVATE(litter_dec_fac)
279	REAL(r_std), SAVE, DIMENSION(nlitt,ncarb,nlevs) :: frac_soil !! Fraction of litter that goes into soil
280	!! (litter -> carbon, above and below). The
281	!! remaining part goes to the atmosphere
282	!$OMP THREADPRIVATE(frac_soil)
283	REAL(r_std), DIMENSION(npts) :: soilhum_decomp !! Humidity used for decompostition in soil
284	!! (unitless)
285	REAL(r_std), DIMENSION(npts) :: fd !! Fraction of structural or metabolic litter
286	!! decomposed (unitless)
287	REAL(r_std), DIMENSION(npts,nelements) :: qd !! Quantity of structural or metabolic litter
288	!! decomposed @tex $(gC m^{-2})$ @endtex
289	!! @tex $(gC m^{-2})$ @endtex
290	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: old_struc !! Old structural litter, above and below
291	!! @tex $(gC m^{-2})$ @endtex
292	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: old_woody !! Old woody litter, above and below
293	!! @tex $(gC m^{-2})$ @endtex
294	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements) :: litter_inc !! Increase of metabolic and structural
295	!! litter, above and below ground. The unit
296	!! is given by m^2 of ground.
297	!! @tex $(gC m^{-2})$ @endtex
298	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: lignin_struc_inc !! Lignin increase in structural litter,
299	!! above and below ground. The unit is given
300	!! by m^2 of ground.
301	!! @tex $(gC m^{-2})$ @endtex
302	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: lignin_wood_inc !! Lignin increase in woody litter,
303	!! above and below ground. The unit is given
304	!! by m^2 of ground.
305	!! @tex $(gC m^{-2})$ @endtex
306	REAL(r_std), DIMENSION(npts) :: zdiff_min !! Intermediate field for looking for minimum
307	!! of what? this is not used in the code.
308	!! [??CHECK] could we delete it?
309	CHARACTER(LEN=10), DIMENSION(nlitt) :: litter_str !! Messages to write output information about
310	!! the litter
311	CHARACTER(LEN=22), DIMENSION(nparts) :: part_str !! Messages to write output information about
312	!! the plant
313	CHARACTER(LEN=7), DIMENSION(ncarb) :: carbon_str !! Messages to write output information about
314	!! the soil carbon
315	CHARACTER(LEN=5), DIMENSION(nlevs) :: level_str !! Messages to write output information about
316	!! the level (aboveground or belowground litter)
317	INTEGER(i_std) :: ilitt,ilev,ibdl !! Indices (unitless)
318	INTEGER(i_std) :: ipts,ivm, j !! Indices (unitless)
319	INTEGER(i_std) :: icarb,inspec,iele !! Indices (unitless)
320	INTEGER(i_std) :: ipar,imbc,iicarb !! Indices (unitless)
321	REAL(r_std), DIMENSION(npts,nvm) :: f_soil_n_min !! soil_n_min function response used for defing C to N target ratios (-)
322	REAL(r_std), DIMENSION(npts,nvm) :: f_pnc !! plant nitrogen concentration function response
323	!! used for defining C to N target ratios (-)
324	INTEGER(i_std) :: itarget !! target som pool
325	REAL(r_std), DIMENSION(npts,nvm,nparts) :: CN !! CN ratio of the litter pools
326	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
327	!! mass balance check for this routine
328	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
329	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
330	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
331	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start and end pool of this routine
332	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
333	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! Start and end pool of this routine
334	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
335	REAL(r_std), DIMENSION(npts,nvm) :: temp
336	REAL(r_std), DIMENSION(npts) :: err !! Array for error checking
337	REAL(r_std), DIMENSION(ncarb,nelements) :: som_input_old !! Used to store som_input values
338	!! in case of inadequate soil N
339	REAL(r_std), DIMENSION(ncarb,nelements) :: delta_som_input !! Used to store som_input values
340	!! in case of inadequate soil N
341	REAL(r_std), DIMENSION(nvm,nelements) :: som_input_new !! Used to recalculate som_input in
342	!! case of inadequate soil N
343	REAL(r_std) :: CN_input_total !! overall C/N ratio of som_input_total
344	REAL(r_std), DIMENSION(nvm,nelements) :: som_input_total !! Used to recalculate som_input in
345	!! case of inadequate soil N
346	REAL(r_std), DIMENSION(nelements) :: litter_total_new !! Updated litter pool in case of nitrogen
347	!! deficit
348	REAL(r_std), DIMENSION(npts,nvm) :: total_init_nitrogen !! Litter pool at the very start of the
349	!! calculations. Required to correct for the
350	!! nitrogen deficit
351	REAL(r_std) :: n_mineralisation_init !! temp variable in case of overspending
352	REAL(r_std) :: delta_hetero_litter !! reduction factor of resp_hetero
353	!! in case of nitrogen shortage
354	REAL(r_std), DIMENSION(npts,nvm,nelements) :: share_litter_struc_a !! Fractions of litter in various pools
355	REAL(r_std), DIMENSION(npts,nvm,nelements) :: share_litter_struc_b !! Fractions of litter in various pools
356	REAL(r_std), DIMENSION(npts,nvm,nelements) :: share_litter_met_a !! Fractions of litter in various pools
357	REAL(r_std), DIMENSION(npts,nvm,nelements) :: share_litter_met_b !! Fractions of litter in various pools
358	REAL(r_std), DIMENSION(npts,nvm,nelements) :: share_litter_woody_a !! Fractions of litter in various pools
359	REAL(r_std), DIMENSION(npts,nvm,nelements) :: share_litter_woody_b !! Fractions of litter in various pools
360	REAL(r_std) :: share_som_input_active !! Active pool fraction of
361	!! som_input(:,:,ivm,initrogen)
362	REAL(r_std) :: share_som_input_slow !! Slow pool fraction of
363	!! som_input(:,:,ivm,initrogen)
364	REAL(r_std) :: share_som_input_surface !! Surface pool fraction of
365	!! som_input(:,:,ivm,initrogen)
366	REAL(r_std), DIMENSION(npts,ncarb,nvm,nlevs) :: N_deficit !! Missing litter N (per pool) after pulling
367	!! additional litter N to satisfy
368	!! N_mineralisation demand (5.5)
369	REAL(r_std) :: N_deficit_total !! Total litter N deficit (5.5)
370	REAL(r_std), DIMENSION(npts,nvm) :: litter_max !! Finds largest current litter pool
371	REAL(r_std), DIMENSION(npts,nvm,ncarb,ncarb,nelements) :: flux !! Fluxes between carbon pools (gC/m^2)
372	REAL(r_std), DIMENSION(npts,nvm,ncarb,nelements) :: fluxtot !! Total flux out of carbon pools (gC/m^2)
373	REAL(r_std), DIMENSION(npts,ncarb,nvm,nelements) :: som_temp !! Temp variable
374	REAL(r_std), DIMENSION(npts,ncarb,ncarb) :: frac_carb !! Flux fractions between carbon pools
375	REAL(r_std), DIMENSION(npts,ncarb) :: frac_resp !! Flux fractions from carbon pools to the atmosphere
376	!! atmosphere (respiration) (unitless, 0-1)
377	REAL(r_std), DIMENSION(ncarb) :: som_turn !! Residence time in SOM pools (days)
378	REAL(r_std), DIMENSION(npts,nvm) :: n_mineralisation_som !! n_mineralisation that we will calculate in the
379	!! next routine
380	REAL(r_std), DIMENSION(npts,nvm) :: n_min_old !! Temp variable
381	REAL(r_std), DIMENSION(nvm,nelements) :: som_input_total_new !! Truncated value of N moved from
382	!! n_mineralisation to som_input
383	REAL(r_std), DIMENSION(npts,nvm) :: share_som_input_act_n !! Active pool fraction of
384	!! som_input(:,:,ivm,initrogen)
385	REAL(r_std), DIMENSION(npts,nvm) :: share_som_input_slow_n !! Slow pool fraction of
386	!! som_input(:,:,ivm,initrogen)
387	REAL(r_std), DIMENSION(npts,nvm) :: share_som_input_surf_n !! Surface pool fraction of
388	!! som_input(:,:,ivm,initrogen)
389	REAL(r_std), DIMENSION(npts,nvm) :: resp_hetero_soil !! Soil heterotrophic respiration (gC/m^2/day)
390	REAL(r_std), DIMENSION(npts,nvm,nelements) :: total_som_n_flux !! Total flux between all soil pools. Used to adjust
391	!! som if overspending @tex $(N gN m^{-2})$ @endtex
392	REAL(r_std), DIMENSION(npts,nvm) :: f_fungivores !! Fraction of decomposed litter consumed by
393	!! fungivores
394	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: biomass_temp !! Temporal value to calculate fungivores
395	INTEGER(i_std), DIMENSION(2) :: ix !! Row/column number of the maximal value
396	INTEGER(i_std) :: ier !! Error handling
397	REAL(r_std), DIMENSION(npts) :: error_count !! Count the number of errors
398	REAL(r_std), DIMENSION(nlitt) :: tree_litterfrac
399	REAL(r_std), DIMENSION(npts) :: residual !! temporary variable to ensure mass balance closure for som_input and hetero_resp
400	REAL(r_std) :: total !! Total of share_litter used as temporary variable (unitless,0-1)
401
402
403	!_ ================================================================================================================================
404
405	IF ( firstcall_litter ) THEN
406	! 1.1 Initialize local printlev
407	printlev_loc=get_printlev('litter')
408
409	END IF
410
411	IF (printlev_loc>=3) WRITE(numout,*) 'Entering littercalc'
412
413	!! 1. Initialisations of the different fields during the first call of the routine
414	dt = dt_sechiba/one_day
415
416	IF ( firstcall_litter ) THEN
417
418	IF ( .NOT. ALLOCATED(litterfrac) ) THEN
419	ALLOCATE(litterfrac(npts,nvm,nparts,nlitt))
420	ENDIF
421
422	!! 1.1.2 residence times in litter pools (days)
423	! .5 years, 3 years and ??30 years??(2.45)
424	litter_dec_fac(imetabolic) = turn_metabolic / one_year
425	litter_dec_fac(istructural) = turn_struct / one_year
426	litter_dec_fac(iwoody) = turn_woody / one_year
427
428	!! 1.1.5 decomposition flux fraction that goes into soil
429	! (litter -> carbon, above and below)
430	! 1-frac_soil goes into atmosphere
431	frac_soil(:,:,:) = zero
432
433	! structural litter: lignin fraction goes into slow pool + respiration,
434	! rest into active pool + respiration
435	frac_soil(istructural,isurface,iabove) = frac_soil_struct_sua
436	frac_soil(istructural,iactive,ibelow) = frac_soil_struct_ab
437	frac_soil(istructural,islow,iabove) = frac_soil_struct_sa
438	frac_soil(istructural,islow,ibelow) = frac_soil_struct_sb
439
440	! metabolic litter: all goes into active pool + respiration.
441	! Nothing into slow or passive pool.
442	frac_soil(imetabolic,isurface,iabove) = frac_soil_metab_sua
443	frac_soil(imetabolic,iactive,ibelow) = frac_soil_metab_ab
444
445	! woody litter lignin fraction should be lower than that of structural
446	! litter, else woody litter will not decompose quickly enough.
447	! In effect, this decreases the carbon use efficiency of the system to
448	! allow for faster litter decomposition. The parameter frac_woody is not
449	! well-constrained and was tested for a boreal PFT in northern Sweden.
450	frac_soil(iwoody,isurface,iabove) = frac_woody * frac_soil_struct_sua
451	frac_soil(iwoody,iactive,ibelow) = frac_woody * frac_soil_struct_ab
452	frac_soil(iwoody,islow,iabove) = frac_woody * frac_soil_struct_sa
453	frac_soil(iwoody,islow,ibelow) = frac_woody * frac_soil_struct_sb
454
455	!! 1.3 messages
456	litter_str(imetabolic) = 'metabolic'
457	litter_str(istructural) = 'structural'
458	litter_str(iwoody) = 'woody'
459
460	carbon_str(iactive) = 'active'
461	carbon_str(isurface) = 'surface'
462	carbon_str(islow) = 'slow'
463	carbon_str(ipassive) = 'passive'
464
465	level_str(iabove) = 'above'
466	level_str(ibelow) = 'below'
467
468	part_str(ileaf) = 'leaves'
469	part_str(isapabove) = 'sap above ground'
470	part_str(isapbelow) = 'sap below ground'
471	part_str(iheartabove) = 'heartwood above ground'
472	part_str(iheartbelow) = 'heartwood below ground'
473	part_str(iroot) = 'roots'
474	part_str(ifruit) = 'fruits'
475	part_str(icarbres) = 'carbohydrate reserve'
476	part_str(ilabile) = 'labile reserve'
477
478	! Debug
479	IF (printlev_loc >= 4) THEN
480	WRITE(numout,*) 'litter:'
481
482	WRITE(numout,*) ' > C/N ratios: '
483	DO ipar = 1, nparts
484	WRITE(numout,*) ' ', part_str(ipar), ': ',CN_fix(ipar)
485	ENDDO
486
487	WRITE(numout,*) ' > scaling depth for decomposition (m): ',z_decomp
488
489	WRITE(numout,*) ' > minimal carbon residence time in litter pools (d):'
490	DO ilitt = 1, nlitt
491	WRITE(numout,*) ' ',litter_str(ilitt),':',litter_dec_fac(ilitt)
492	ENDDO
493
494	WRITE(numout,*) ' > litter decomposition flux fraction that really goes '
495	WRITE(numout,*) ' into carbon pools (rest into the atmosphere):'
496	DO ilitt = 1, nlitt
497	DO ilev = 1, nlevs
498	DO icarb = 1, ncarb
499	WRITE(numout,*) ' ',litter_str(ilitt),' ',level_str(ilev),' -> ',&
500	carbon_str(icarb),':', frac_soil(ilitt,icarb,ilev)
501	ENDDO
502	ENDDO
503	ENDDO
504
505	ENDIF
506	!-
507
508	firstcall_litter = .FALSE.
509
510	ENDIF ! firstcall_litter
511
512	!! 1.1.3 litter fractions:
513	! What fraction of leaves, wood, etc. goes into metabolic and
514	! structural litterpools
515	! Initialised for when everything is printed below
516	litterfrac(:,:,:,:) = zero
517	CN(:,:,:) = zero
518	CN_target(:,:,:)= zero
519
520	!+++CHECK+++
521	! Add documentation
522	! Calculate the CN ratio. If the CN ratio cannot be calculated
523	! then use the prescribed value
524	DO ipar = 1, nparts
525
526	! Debug
527	IF(printlev_loc >=4) THEN
528	WRITE(numout,*) 'bm_to_litter carbon for ',&
529	part_str(ipar),':',bm_to_litter(test_grid,test_pft,ipar,icarbon)
530	WRITE(numout,*) 'bm_to_litter nitrogen for ',&
531	part_str(ipar),':',bm_to_litter(test_grid,test_pft,ipar,initrogen)
532	WRITE(numout,*) 'turnover carbon for ',&
533	part_str(ipar),':',turnover(test_grid,test_pft,ipar,icarbon)
534	WRITE(numout,*) 'turnover nitrogen for ',&
535	part_str(ipar),':',turnover(test_grid,test_pft,ipar,initrogen)
536	ENDIF
537	!-
538
539	! Note that all the carbon and nitrogen contained in tree_bm_to_litter
540	! is also contained in bm_to_litter. tree_bm_to_litter is only needed
541	! to ensure that some of the woody litter that moved to bare soil, crop/grass
542	! PFTs during LCC is correctly dealt with. tree_bm_to_litter should
543	! NOT be accounted for in most of the calculations.
544	WHERE((bm_to_litter(:,:,ipar,initrogen)+turnover(:,:,ipar,initrogen)).GT.min_stomate)
545	CN(:,:,ipar) = &
546	(bm_to_litter(:,:,ipar,icarbon)+turnover(:,:,ipar,icarbon)) / &
547	(bm_to_litter(:,:,ipar,initrogen)+turnover(:,:,ipar,initrogen))
548
549	ELSEWHERE
550
551	CN(:,:,ipar) = CN_fix(ipar)
552
553	ENDWHERE
554
555	DO ivm = 1,nvm
556
557	IF ( ((ipar == isapabove) .OR. (ipar == isapbelow) .OR. &
558	(ipar == iheartabove) .OR. (ipar == iheartbelow)) .AND. &
559	ivm == ibare_sechiba) THEN
560
561	WHERE (bm_to_litter(:,ivm,ipar,icarbon).GT.10*EPSILON(un))
562
563	! PFT 1 with woody litter will be treated as
564	! a forest. The threshold needs to be much lower than min_stomate.
565	! If not, mass is being lost. Because this calculation is repeated
566	! 48 times, small errors (1e-9) soon exceed 1e-8 and make the model
567	! crash in the consistency cross-check for nbp.
568	litterfrac(:,ivm,ipar,iwoody) = 1.0
569	litterfrac(:,ivm,ipar,imetabolic) = zero
570	litterfrac(:,ivm,ipar,istructural) = zero
571
572	ENDWHERE
573
574	ELSEIF ( ((ipar == isapabove) .OR. (ipar == isapbelow) .OR. &
575	(ipar == iheartabove) .OR. (ipar == iheartbelow)) &
576	.AND. is_tree(ivm) ) THEN
577
578	! Woody components of forest PFTs
579	litterfrac(:,ivm,ipar,iwoody) = 1.0
580	litterfrac(:,ivm,ipar,imetabolic) = zero
581	litterfrac(:,ivm,ipar,istructural) = zero
582
583	ELSE
584
585	! PFT 1 without woody litter and all non-forest PFTs
586	IF( (ipar == icarbres) .OR. (ipar == ilabile)) THEN
587
588	litterfrac(:,ivm,ipar,imetabolic) = 1.0
589
590	ELSE
591
592	litterfrac(:,ivm,ipar,imetabolic) = &
593	MAX(metabolic_ref_frac - metabolic_LN_ratio * &
594	LC(ivm,ipar) * CN(:,ivm,ipar),zero)
595
596	ENDIF
597
598	litterfrac(:,ivm,ipar,istructural) = &
599	1. - litterfrac(:,ivm,ipar,imetabolic)
600	litterfrac(:,ivm,ipar,iwoody) = zero
601
602	ENDIF
603
604	END DO
605
606	END DO
607
608	! Error checking
609	IF(err_act.GT.1)THEN
610	DO ipts = 1,npts
611	DO ivm = 1,nvm
612	DO ipar = 1, nparts
613	IF((bm_to_litter(ipts,ivm,ipar,initrogen) + &
614	turnover(ipts,ivm,ipar,initrogen)).GT.min_stomate)THEN
615	IF (ABS(SUM(litterfrac(ipts,ivm,ipar,:))-1.0).GT.10*EPSILON(un)) THEN
616	WRITE(numout,*) 'litterfrac - should equal to 1, ', ivm, ipar, &
617	SUM(litterfrac(ipts,ivm,ipar,:))
618	WRITE(numout,*) 'CN, ', CN(ipts,ivm,ipar)
619	WRITE(numout,*) 'CN_target (ncarb), ',CN_target(ipts,ivm,:)
620	CALL ipslerr(3,'stomate_litter',&
621	'litterfrac do not add up to 1','','')
622	ENDIF
623	ENDIF
624	ENDDO
625	ENDDO
626	ENDDO
627	ENDIF
628
629	! Error checking
630	! soil_n_min should never be zero. It is a nitrogen pool
631	! expressed in gN m-2
632	temp(:,:)=zero
633	WHERE(soil_n_min(:,:,iammonium).LT.zero)
634	temp(:,:) = un
635	ENDWHERE
636
637	IF (SUM(SUM(temp(:,:),2)).GT.zero) THEN
638	WRITE(numout,*)'ERROR: soil_n_min iammonium is negative',&
639	soil_n_min(:,:,iammonium)
640	CALL ipslerr(3,'stomate_litter',&
641	'At least one soil_n_min value',&
642	'is negative in stomate_litter','')
643	ENDIF
644
645	temp(:,:)=zero
646	WHERE(soil_n_min(:,:,initrate).LT.zero)
647	temp(:,:) = un
648	ENDWHERE
649
650	IF (SUM(SUM(temp(:,:),2)).GT.zero) THEN
651	WRITE(numout,*)'ERROR: soil_n_min initrate is negative', &
652	soil_n_min(:,:,initrate)
653	CALL ipslerr(3,'stomate_litter',&
654	'At least one soil_n_min value',&
655	'is negative in stomate_litter','')
656	ENDIF
657
658	! Error checking
659	! This code is OK only when the land cover change and the age
660	! class distribution code are working properly. Better to test
661	! it in those subroutines than to enforce it here. If there are
662	! problems with land cover change or age class distribution the
663	! mass balance for soil_n_min may be violated here.
664	DO ipts = 1, npts
665	DO ivm = 1, nvm
666	DO inspec = 1, nnspec
667	IF((veget_max(ipts,ivm).LT.min_stomate) .AND. &
668	(soil_n_min(ipts,ivm,inspec).GT. min_stomate)) THEN
669	WRITE(numout,*) 'PFT,',ivm,', pixel:', ipts,', n-species:',inspec
670	WRITE(numout,*) 'veget_max,',veget_max(ipts,ivm)
671	WRITE(numout,*) 'soil_n_min, ',soil_n_min(ipts,ivm,inspec)
672	CALL ipslerr_p(3,'stomate_litter',&
673	'soil_n_min present in pixels without veget_max',&
674	'','')
675	ENDIF
676	ENDDO
677	ENDDO
678	ENDDO
679
680	!+++CHECK+++
681	! ADD documentation. What justifies the use of 2.? Should
682	! that parameter be externalized? What does that parameter
683	! represent? The MIN does NOT protect against negative
684	! soil_n_values. Negative values should no longer occur!
685	f_soil_n_min(:,:) = &
686	MIN((soil_n_min(:,:,iammonium) + soil_n_min(:,:,initrate)), 2.)
687	!+++++++++++
688
689	! C to N target ratios of different pools
690	CN_target(:,:,:)=zero
691
692	! CN ratio ranges between 3 and 15 for active pool
693	! Figure 4 of Parton et al. (1993) show lower CN values for active
694	! pool of 2 rather than 3
695	CN_target(:,:,iactive)= CN_target_iactive_ref + &
696	CN_target_iactive_Nmin * f_soil_n_min(:,:)
697
698	! CN ratio ranges between 12 and 20 for slow pool
699	CN_target(:,:,islow)= CN_target_islow_ref + &
700	CN_target_islow_Nmin * f_soil_n_min(:,:)
701
702	! CN ratio ranges between 7 and 10 for passive pool
703	! OCN uses a fixed value of 9. (don't know why)
704	! Figure 4 of Parton et al. (1993) show lower CN values
705	! for passive pool of 3 rather than 7
706	CN_target(:,:,ipassive)= CN_target_ipassive_ref + &
707	CN_target_ipassive_Nmin * f_soil_n_min(:,:)
708
709	! Litter nitrogen content (%)
710	WHERE(litter(:,imetabolic,:,iabove,icarbon)+&
711	litter(:,istructural,:,iabove,icarbon).GT.zero)
712
713	!+++CHECK+++
714	! Why do we multiply with 2.?, why do we take the
715	! minimun between the calculation and 2. Do those
716	! two 2s represent the same parameter?
717	f_pnc(:,:)=MIN((litter(:,imetabolic,:,iabove,initrogen)+&
718	litter(:,istructural,:,iabove,initrogen)) / &
719	(2.*(litter(:,imetabolic,:,iabove,icarbon)+&
720	litter(:,istructural,:,iabove,icarbon))) * 100., 2.)
721	!+++++++++++
722
723	ELSEWHERE
724
725	f_pnc(:,:)= zero
726
727	ENDWHERE
728
729	!+++CHECK+++
730	! Add documentation
731	CN_target(:,:,isurface) = CN_target_isurface_ref + &
732	CN_target_isurface_pnc * f_pnc(:,:)
733	!+++++++++++
734
735	!! 1.4 set output to zero
736	deadleaf_cover(:) = zero
737	resp_hetero_litter(:,:) = zero
738	som_input(:,:,:,:) = zero
739	n_fungivores(:,:) = zero
740	f_fungivores(:,:) = zero
741
742	!! 1.5 Initialize check for mass balance closure
743	! The mass balance is calculated at the end of this routine
744	! in section 6. Initial biomass and harvest pool all other
745	! relevant pools were just set to zero.
746	IF (err_act.GT.1) THEN
747
748	pool_start(:,:,:) = zero
749	pool_start(:,:,initrogen) = pool_start(:,:,initrogen) + &
750	n_mineralisation(:,:) * veget_max(:,:)
751	pool_start(:,:,initrogen) = pool_start(:,:,initrogen) + &
752	n_fungivores(:,:) * veget_max(:,:)
753	pool_start(:,:,initrogen) = pool_start(:,:,initrogen) + &
754	(soil_n_min(:,:,iammonium) + soil_n_min(:,:,initrate))*veget_max(:,:)
755
756	DO iele = 1,nelements
757	! The litter pool
758	DO ilitt = 1,nlitt
759	DO ilev = 1,nlevs
760	pool_start(:,:,iele) = pool_start(:,:,iele) + &
761	(litter(:,ilitt,:,ilev,iele) * veget_max(:,:))
762	ENDDO
763	ENDDO
764
765	! Pools added to the litter
766	! Note that all the carbon and nitrogen conatined in tree_bm_to_litter
767	! is also contained in bm_to_litter. tree_bm_to_litter is only needed
768	! to ensure that some of the woody litter that moved to a crop/grass
769	! PFT during LCC is correctly dealt with. tree_bm_to_litter should
770	! NOT be accounted for in most of the calculations.
771	DO ipar = 1,nparts
772	pool_start(:,:,iele) = pool_start(:,:,iele) + &
773	(turnover(:,:,ipar,iele) + &
774	bm_to_litter(:,:,ipar,iele)) * veget_max(:,:)
775	ENDDO
776
777	DO ivm = 1, nvm
778	pool_start(:,ivm,iele) = pool_start(:,ivm,iele) + &
779	SUM(harvest_pool_acc(:,ivm,:,iele),2)/area(:)
780	END DO
781
782	ENDDO
783
784	ENDIF ! err_act.GT.1
785
786	! Store the initial available N. In case there is overspending of nitrogen
787	! this value can be used to adjust the spending
788	WHERE (veget_max(:,:).GT.min_stomate)
789	total_init_nitrogen(:,:) = n_mineralisation(:,:) + &
790	SUM(SUM(litter(:,:,:,:,initrogen),4),2) + SUM(turnover(:,:,:,initrogen),3) + &
791	SUM(bm_to_litter(:,:,:,initrogen),3)
792	ELSEWHERE
793	total_init_nitrogen(:,:) = zero
794	ENDWHERE
795
796	! fungivores were introduced to avoid the problem that soil mineral nitrogen is
797	! used up by the decomposition after all plants within a PFT die. A certain portion
798	! (f_fungivores) of the nitrogen pool from the litter was fed to the vegetation
799	! using n_fungivores. But not to give extra
800	! nitrogen even in optimum condition, dynamic f_fungivores is more reasonable
801	! than fixed f_fungivores. Before calculating n_fungivores f_fungivoes is
802	! calculated using total biomass and total litter of carbon.
803	! Underlying idea is to use predator prey relationship between fungi and
804	! fungivores. When there is a lot of litter, there will be a lot of fungi and
805	! so we expect a lot of fungivores. Nitrogen excreted by the fungivores will
806	! be used by plants. This implementation would benefit from a more solid
807	! process understanding and more literature on the topic.
808	biomass_temp(:,:,:,:) = cc_to_biomass(npts,nvm,circ_class_biomass(:,:,:,:,:),&
809	circ_class_n(:,:,:))
810
811	WHERE ((SUM(biomass_temp(:,:,:,icarbon),3) + SUM(SUM(litter(:,:,:,:,icarbon),4),2) + &
812	SUM(bm_to_litter(:,:,:,icarbon),3) ) .LT. min_stomate)
813	f_fungivores(:,:) = zero
814	ELSEWHERE
815	f_fungivores(:,:) = (SUM(SUM(litter(:,:,:,:,icarbon),4),2) + &
816	SUM(bm_to_litter(:,:,:,icarbon),3))/(SUM(biomass_temp(:,:,:,icarbon),3) + &
817	SUM(SUM(litter(:,:,:,:,icarbon),4),2) + SUM(bm_to_litter(:,:,:,icarbon),3) )
818	ENDWHERE
819
820	!+++HACK+++
821	! Some initial tests show that f_fungivores result in a too high N-availability
822	! at the start of a simulation and that following unidentified model developments
823	! the model is now able to regrow vegetation right after a clear cut.
824	! f_fungivores was left in the code just in case we jumped to conclusion too quickly
825	! and we want to use it anyway but a recent proposal that got funded will add
826	! biomass pools of the different functional groups into the model. At that point
827	! the crude current implementation of f_fungivores will have to be revised anyway.
828	f_fungivores = zero
829	!++++++++++
830
831	!! 2. Add biomass to different litterpools (per m^2 of ground)
832
833	!! 2.1 first, save old structural litter (needed for lignin fractions).
834	! above/below
835	DO ilev = 1, nlevs !Loop over litter levels (above and below ground)
836
837	DO ivm = 1,nvm !Loop over PFTs
838
839	old_struc(:,ivm,ilev) = litter(:,istructural,ivm,ilev,icarbon)
840	old_woody(:,ivm,ilev) = litter(:,iwoody,ivm,ilev,icarbon)
841
842	ENDDO
843
844	ENDDO
845
846	!! 2.2 update litter, dead leaves, and lignin content in structural litter
847	litter_inc(:,:,:,:,:) = zero
848	lignin_struc_inc(:,:,:) = zero
849	lignin_wood_inc(:,:,:) = zero
850
851	tree_litterfrac(iwoody) = 1.0
852	tree_litterfrac(imetabolic) = zero
853	tree_litterfrac(istructural) = zero
854
855	DO ivm = 1,nvm !Loop over PFTs
856
857	!! 2.2.1 litter
858	DO ilitt = 1, nlitt !Loop over litter pools (metabolic and structural)
859
860	DO iele = 1, nelements ! Loop over element pools (carbon and nitrogen)
861
862	!! 2.2.2 calculate litter increase (per m^2 of ground).
863	! Only a given fracion of fruit turnover is directly
864	! coverted into litter. Litter increase for each PFT,
865	! structural and metabolic, above/below
866	litter_inc(:,ilitt,ivm,iabove,iele) = &
867	litterfrac(:,ivm,ileaf,ilitt) * bm_to_litter(:,ivm,ileaf,iele) + &
868	litterfrac(:,ivm,isapabove,ilitt) * &
869	(bm_to_litter(:,ivm,isapabove,iele)-tree_bm_to_litter(:,ivm,isapabove,iele)) + &
870	litterfrac(:,ivm,iheartabove,ilitt) * &
871	(bm_to_litter(:,ivm,iheartabove,iele)-tree_bm_to_litter(:,ivm,iheartabove,iele)) + &
872	tree_litterfrac(ilitt) * tree_bm_to_litter(:,ivm,isapabove,iele) + &
873	tree_litterfrac(ilitt) * tree_bm_to_litter(:,ivm,iheartabove,iele) + &
874	litterfrac(:,ivm,ifruit,ilitt) * bm_to_litter(:,ivm,ifruit,iele) + &
875	litterfrac(:,ivm,icarbres,ilitt) * bm_to_litter(:,ivm,icarbres,iele) + &
876	litterfrac(:,ivm,ilabile,ilitt) * bm_to_litter(:,ivm,ilabile,iele) + &
877	litterfrac(:,ivm,ileaf,ilitt) * turnover(:,ivm,ileaf,iele) + &
878	litterfrac(:,ivm,isapabove,ilitt) * turnover(:,ivm,isapabove,iele) + &
879	litterfrac(:,ivm,iheartabove,ilitt) * turnover(:,ivm,iheartabove,iele) + &
880	litterfrac(:,ivm,ifruit,ilitt) * turnover(:,ivm,ifruit,iele) + &
881	litterfrac(:,ivm,icarbres,ilitt) * turnover(:,ivm,icarbres,iele)+ &
882	litterfrac(:,ivm,ilabile,ilitt) * turnover(:,ivm,ilabile,iele)
883
884	litter_inc(:,ilitt,ivm,ibelow,iele) = &
885	litterfrac(:,ivm,isapbelow,ilitt) * &
886	(bm_to_litter(:,ivm,isapbelow,iele)-tree_bm_to_litter(:,ivm,isapbelow,iele)) + &
887	litterfrac(:,ivm,iheartbelow,ilitt) * &
888	(bm_to_litter(:,ivm,iheartbelow,iele)-tree_bm_to_litter(:,ivm,iheartbelow,iele)) + &
889	tree_litterfrac(ilitt) * tree_bm_to_litter(:,ivm,isapbelow,iele) + &
890	tree_litterfrac(ilitt) * tree_bm_to_litter(:,ivm,iheartbelow,iele) + &
891	litterfrac(:,ivm,iroot,ilitt) * bm_to_litter(:,ivm,iroot,iele) + &
892	litterfrac(:,ivm,isapbelow,ilitt) * turnover(:,ivm,isapbelow,iele) + &
893	litterfrac(:,ivm,iheartbelow,ilitt) * turnover(:,ivm,iheartbelow,iele) + &
894	litterfrac(:,ivm,iroot,ilitt) * turnover(:,ivm,iroot,iele)
895
896	ENDDO
897
898	IF ( ilitt /= iwoody ) THEN
899
900	!! 2.2.3 dead leaves, for soil cover.
901	dead_leaves(:,ivm,ilitt) = &
902	dead_leaves(:,ivm,ilitt) + litterfrac(:,ivm,ileaf,ilitt) * &
903	( bm_to_litter(:,ivm,ileaf,icarbon) + &
904	turnover(:,ivm,ileaf,icarbon) )
905
906	ENDIF
907
908	IF(ivm == ibare_sechiba.OR.is_tree(ivm))THEN
909
910	! Note that a bare soil pixel may have some left
911	! overs of the litter of a previous forest. So
912	! bare soil is treated as forest here.
913
914	IF ( ilitt == istructural ) THEN
915
916	!! 2.2.4 lignin increase in structural litter
917	lignin_struc_inc(:,ivm,iabove) = &
918	LC(ivm,ileaf) * bm_to_litter(:,ivm,ileaf,icarbon) + &
919	LC(ivm,ifruit) * bm_to_litter(:,ivm,ifruit,icarbon) + &
920	LC(ivm,icarbres) * bm_to_litter(:,ivm,icarbres,icarbon) + &
921	LC(ivm,ilabile) * bm_to_litter(:,ivm,ilabile,icarbon) + &
922	LC(ivm,ileaf) * turnover(:,ivm,ileaf,icarbon) + &
923	LC(ivm,ifruit) * turnover(:,ivm,ifruit,icarbon) + &
924	LC(ivm,icarbres) * turnover(:,ivm,icarbres,icarbon) + &
925	LC(ivm,ilabile) * turnover(:,ivm,ilabile,icarbon)
926
927	lignin_struc_inc(:,ivm,ibelow) = &
928	LC(ivm,iroot) * bm_to_litter(:,ivm,iroot,icarbon) + &
929	LC(ivm,iroot) * turnover(:,ivm,iroot,icarbon)
930
931	ELSEIF ( ilitt == iwoody) THEN
932
933	!! 2.2.4 lignin increase in woody litter
934	lignin_wood_inc(:,ivm,iabove) = &
935	LC(ivm,isapabove) * turnover(:,ivm,isapabove,icarbon) + &
936	LC(ivm,iheartabove) * turnover(:,ivm,iheartabove,icarbon)
937
938	lignin_wood_inc(:,ivm,ibelow) = &
939	LC(ivm,isapbelow) * turnover(:,ivm,isapbelow,icarbon) + &
940	LC(ivm,iheartbelow) * turnover(:,ivm,iheartbelow,icarbon)
941
942	ENDIF ! istructural or iwoody
943
944	ELSE
945
946	! None woody vegetation
947	IF (ilitt == istructural) THEN
948
949	!! 2.2.4 lignin increase in structural litter
950	lignin_struc_inc(:,ivm,iabove) = lignin_struc_inc(:,ivm,iabove) + &
951	LC(ivm,isapabove) * turnover(:,ivm,isapabove,icarbon) + &
952	LC(ivm,iheartabove) * turnover(:,ivm,iheartabove,icarbon)
953
954	lignin_struc_inc(:,ivm,ibelow) = lignin_struc_inc(:,ivm,ibelow) + &
955	LC(ivm,isapbelow)*turnover(:,ivm,isapbelow,icarbon) + &
956	LC(ivm,iheartbelow)*turnover(:,ivm,iheartbelow,icarbon)
957
958	ENDIF ! istructural
959
960	ENDIF ! is_tree
961
962	IF (ilitt == istructural) THEN
963
964	lignin_wood_inc(:,ivm,iabove) = &
965	LC(ivm,isapabove) * tree_bm_to_litter(:,ivm,isapabove,icarbon) + &
966	LC(ivm,iheartabove) * tree_bm_to_litter(:,ivm,iheartabove,icarbon)
967
968	lignin_wood_inc(:,ivm,ibelow) = &
969	LC(ivm,isapbelow) * tree_bm_to_litter(:,ivm,isapbelow,icarbon) + &
970	LC(ivm,iheartbelow) * tree_bm_to_litter(:,ivm,iheartbelow,icarbon)
971
972	lignin_struc_inc(:,ivm,iabove) = lignin_struc_inc(:,ivm,iabove) + &
973	LC(ivm,isapabove) * (bm_to_litter(:,ivm,isapabove,icarbon)-tree_bm_to_litter(:,ivm,isapabove,icarbon)) + &
974	LC(ivm,iheartabove) * (bm_to_litter(:,ivm,iheartabove,icarbon)-tree_bm_to_litter(:,ivm,iheartabove,icarbon))
975
976	lignin_struc_inc(:,ivm,ibelow) = lignin_struc_inc(:,ivm,ibelow) + &
977	LC(ivm,isapbelow) * (bm_to_litter(:,ivm,isapbelow,icarbon)-tree_bm_to_litter(:,ivm,isapbelow,icarbon)) + &
978	LC(ivm,iheartbelow) * (bm_to_litter(:,ivm,iheartbelow,icarbon)-bm_to_litter(:,ivm,iheartbelow,icarbon))
979
980	END IF
981
982	ENDDO ! #nlitt
983
984	! Error checking
985	DO iele = 1,nelements
986	! bm_to_litter and turnover_daily should be completely transfered to litter_inc
987	error_count(:) = zero
988	WHERE ( ABS(SUM(SUM(litter_inc(:,:,ivm,:,iele),2)) - &
989	SUM(bm_to_litter(:,ivm,:,iele),2) - &
990	SUM(turnover(:,ivm,:,iele),2)) .GT. 10000*EPSILON(un) )
991	error_count(:) = un
992	END WHERE
993	IF (SUM(error_count(:)).GT.zero) THEN
994	DO ipts = 1,npts
995	IF ( ABS(SUM(SUM(litter_inc(ipts,:,ivm,:,iele),2)) - &
996	SUM(bm_to_litter(ipts,ivm,:,iele)) - &
997	SUM(turnover(ipts,ivm,:,iele))) .GT. 10000*EPSILON(un) ) THEN
998	! Too much carbon/nitrogen went missing. Seems too big too ignore
999	WRITE(numout,*) 'ipts, ivm, ielement, ', ipts, ivm, iele
1000	WRITE(numout,*) 'Sum of litter_inc, ', &
1001	SUM(SUM(litter_inc(ipts,:,ivm,:,iele),2))
1002	WRITE(numout,*) 'Sum of bm_to_litter and turnover, ', &
1003	SUM(bm_to_litter(ipts,ivm,:,iele)) - SUM(turnover(ipts,ivm,:,iele))
1004	WRITE(numout,*) 'Diff, ', SUM(SUM(litter_inc(ipts,:,ivm,:,iele),2)) - &
1005	SUM(bm_to_litter(ipts,ivm,:,iele)) - SUM(turnover(ipts,ivm,:,iele))
1006	CALL ipslerr(3,'littercalc','more than 1e-10 carbon or nitrogen', &
1007	'has gone missing','this might hint at a mass balance problem')
1008	END IF
1009	END DO ! error_count .GT. zero
1010	END IF
1011
1012	! The mismatch is likely a precision error. Fix the issue rather than
1013	! stopping the model
1014	error_count(:) = zero
1015	WHERE ( ABS(SUM(SUM(litter_inc(:,:,ivm,:,iele),2)) - &
1016	SUM(bm_to_litter(:,ivm,:,iele),2) - &
1017	SUM(turnover(:,ivm,:,iele),2)) .GT. 100*EPSILON(un) )
1018	error_count(:) = un
1019	END WHERE
1020	IF (SUM(error_count(:)).GT.zero) THEN
1021	DO ipts = 1,npts
1022	IF ( ABS(SUM(SUM(litter_inc(ipts,:,ivm,:,iele),2)) - &
1023	SUM(bm_to_litter(ipts,ivm,:,iele)) - &
1024	SUM(turnover(ipts,ivm,:,iele))) .GT. 100*EPSILON(un) ) THEN
1025	! Looks like a precision error, fix it.
1026	ix = MAXLOC(litter_inc(ipts,:,ivm,:,iele))
1027	litter_inc(ipts,ix(1),ivm,ix(2),iele) = &
1028	litter_inc(ipts,ix(1),ivm,ix(2),iele) + &
1029	SUM(bm_to_litter(ipts,ivm,:,iele)) + &
1030	SUM(turnover(ipts,ivm,:,iele)) - &
1031	SUM(SUM(litter_inc(ipts,:,ivm,:,iele),2))
1032	IF (litter_inc(ipts,ix(1),ivm,ix(2),iele).LT. zero) THEN
1033	! If this really was a precision error, it should be
1034	! a small correction compared to the actual value of
1035	! litter_inc. If litter_inc becomes zero, the correction
1036	! wasn that small after all. Check.
1037	WRITE(numout,*) 'ipts, ivm, ielement, ', ipts, ivm, iele
1038	WRITE(numout,*) 'Sum of corrected litter_inc, ', &
1039	SUM(SUM(litter_inc(ipts,:,ivm,:,iele),2))
1040	WRITE(numout,*) 'Sum of bm_to_litter and turnover, ', &
1041	SUM(bm_to_litter(ipts,ivm,:,iele)) - &
1042	SUM(turnover(ipts,ivm,:,iele))
1043	CALL ipslerr(3,'littercalc','Trying to correct what was', &
1044	'thought to be a precision error',&
1045	'The negative values suggests there is a real problem')
1046	END IF
1047	END IF
1048	END DO ! npts
1049	ENDIF ! error_count .GT. zero
1050	END DO ! nelements
1051	!-
1052
1053	ENDDO !#nvm
1054
1055	!! 2.2.5 add new litter (struct/met, above/below)
1056	! litter_inc should be a positive number but due to all the multiplications
1057	! above it is possible that very small negative numbers are generated (10e-21)
1058	! When added to a very small number for litter, the total litter could become
1059	! a negative number representing zero. That could fail some of the quality
1060	! tests further down. This is a safe place to enhance numerical consistency of
1061	! the code.
1062	WHERE (litter_inc(:,:,:,:,:).LT.10*EPSILON(un))
1063	litter_inc(:,:,:,:,:) = zero
1064	ENDWHERE
1065	litter(:,:,:,:,:) = litter(:,:,:,:,:) + litter_inc(:,:,:,:,:)
1066
1067	!! 2.2.6 for security: can't add more lignin than structural
1068	!! litter (above/below)
1069	DO ilev = 1, nlevs !Loop over litter levels (above and below ground)
1070
1071	DO ivm = 1,nvm !Loop over PFTs
1072
1073	lignin_struc_inc(:,ivm,ilev) = MIN( lignin_struc_inc(:,ivm,ilev), &
1074	litter_inc(:,istructural,ivm,ilev,icarbon) )
1075	lignin_wood_inc(:,ivm,ilev) = MIN( lignin_wood_inc(:,ivm,ilev), &
1076	litter_inc(:,iwoody,ivm,ilev,icarbon) )
1077
1078	ENDDO ! nvm
1079
1080	ENDDO ! nlevs
1081
1082	!! 2.2.7 new lignin content: add old lignin and lignin increase, divide by
1083	!! total structural litter (above/below)
1084	DO ilev = 1, nlevs !Loop over litter levels (above and below ground)
1085
1086	DO ivm = 1,nvm !Loop over PFTs
1087
1088	WHERE( litter(:,istructural,ivm,ilev,icarbon) .GT. min_stomate )
1089
1090	lignin_struc(:,ivm,ilev) = ( lignin_struc(:,ivm,ilev) * &
1091	old_struc(:,ivm,ilev) + lignin_struc_inc(:,ivm,ilev) ) / &
1092	litter(:,istructural,ivm,ilev,icarbon)
1093
1094	ELSEWHERE
1095
1096	lignin_struc(:,ivm,ilev) = zero
1097
1098	ENDWHERE
1099
1100	WHERE ( litter(:,iwoody,ivm,ilev,icarbon) .GT. min_stomate )
1101
1102	lignin_wood(:,ivm,ilev) = ( lignin_wood(:,ivm,ilev) * &
1103	old_woody(:,ivm,ilev) + lignin_wood_inc(:,ivm,ilev) ) / &
1104	litter(:,iwoody,ivm,ilev,icarbon)
1105
1106	ELSEWHERE
1107
1108	lignin_wood(:,ivm,ilev) = zero
1109
1110	ENDWHERE
1111
1112	ENDDO ! ivm
1113	ENDDO ! ilev
1114
1115	!! 2.2.8 Add the manure into the metabolic above ground pool if enough C is
1116	! harvested over the pixel. If not enough C is harvested a part goes as
1117	! litter and the other part goes at mineral N. Here we assume that N is not
1118	! limiting, only the C harvested is limiting
1119	! Note that here we only treat input for imanure because it is added to the
1120	! litter pool. The other nitrogen input are added to the soil and therefore
1121	! taken care of in stomate_soilcarbon (nitrogen_dynamics).
1122	DO ivm=1,nvm
1123	WHERE((veget_max(:,ivm).GT.min_stomate) .AND. &
1124	(harvest_pool_acc(:,ivm,1,icarbon) .GE. n_input(:,ivm,month,imanure)cn_ratio_manuredtveget_max(:,ivm)area(:)))
1125
1126	litter(:,imetabolic,ivm,iabove,icarbon) = litter(:,imetabolic,ivm,iabove,icarbon) + &
1127	n_input(:,ivm,month,imanure)cn_ratio_manuredt
1128	litter(:,imetabolic,ivm,iabove,initrogen) = litter(:,imetabolic,ivm,iabove,initrogen) + &
1129	n_input(:,ivm,month,imanure)*dt
1130
1131	harvest_pool_acc(:,ivm,1,icarbon) = harvest_pool_acc(:,ivm,1,icarbon) - &
1132	n_input(:,ivm,month,imanure)cn_ratio_manuredtveget_max(:,ivm)area(:)
1133
1134	! Account for the manure when litter decomposition will have to be recalculated
1135	! further down this subroutine
1136	total_init_nitrogen(:,ivm) = total_init_nitrogen(:,ivm) + n_input(:,ivm,month,imanure)*dt
1137
1138	ELSEWHERE((veget_max(:,ivm).GT.min_stomate) .AND. &
1139	(harvest_pool_acc(:,ivm,1,icarbon) .LT. n_input(:,ivm,month,imanure)cn_ratio_manuredtveget_max(:,ivm)area(:)))
1140
1141	litter(:,imetabolic,ivm,iabove,icarbon) = litter(:,imetabolic,ivm,iabove,icarbon) + &
1142	harvest_pool_acc(:,ivm,1,icarbon)/veget_max(:,ivm)/area(:)
1143	litter(:,imetabolic,ivm,iabove,initrogen) = litter(:,imetabolic,ivm,iabove,initrogen) + &
1144	harvest_pool_acc(:,ivm,1,icarbon)/veget_max(:,ivm)/area(:)/cn_ratio_manure
1145
1146	soil_n_min(:,ivm,iammonium) = soil_n_min(:,ivm,iammonium) + &
1147	(n_input(:,ivm,month,imanure)*dt - &
1148	harvest_pool_acc(:,ivm,1,icarbon)/veget_max(:,ivm)/area(:)/cn_ratio_manure)*ratio_nh4_fert
1149	soil_n_min(:,ivm,initrate) = soil_n_min(:,ivm,initrate) + &
1150	(n_input(:,ivm,month,imanure)*dt - &
1151	harvest_pool_acc(:,ivm,1,icarbon)/veget_max(:,ivm)/area(:)/cn_ratio_manure)*(1.-ratio_nh4_fert)
1152
1153	! Account for the manure when litter decomposition will have to be recalculated
1154	! further down this subroutine
1155	total_init_nitrogen(:,ivm) = total_init_nitrogen(:,ivm) + &
1156	harvest_pool_acc(:,ivm,1,icarbon)/veget_max(:,ivm)/area(:)/cn_ratio_manure
1157	harvest_pool_acc(:,ivm,1,icarbon) = zero
1158
1159	ELSEWHERE
1160
1161	! Too small fraction, the n_input is not used
1162	n_input(:,ivm,month,imanure)=zero
1163	ENDWHERE
1164	ENDDO
1165
1166
1167	!! 3. Temperature control on decay: Factor between 0 and 1
1168	!! 3.1 above: surface temperature
1169	control_temp(:,iabove) = control_temp_func(npts, tsurf)
1170
1171	!! 3.2 below: convolution of temperature and decomposer profiles
1172	!! exponential decomposer profile supposed)
1173	!! 3.2.1 rpc is an integration constant such that the integral of
1174	! the decomposer profile becomes one. This looks like the calculation
1175	! of a root profile but it uses a different constant; z_decomp instead
1176	! of humncste for the root profile. The integral goes from the top
1177	! of the soil all the way to the bottom.
1178	! All pixels have the same zdr profile and all PFTs and pixels
1179	! are assumed to have the same z_decomp so rpc is globally fixed.
1180	rpc = un / (EXP(-zdr(0)/z_decomp) - EXP(-zdr(nslm)/z_decomp))
1181
1182	!! 3.2.2 integrate over the nslm levels
1183	tsoil_decomp(:) = zero
1184
1185	DO ibdl = 1, nslm
1186
1187	! Calculate the soil temperature weighted by the assumed presence of
1188	! the decomposers. Presence of the decomposers is prescribes by the
1189	! presence = e^(-soil_depth/z_decomp) where z_decomp described the
1190	! shape of an exponentially decaying function with soil depth. stempdiag
1191	! is descritized along the nodes of the soil layers following the scheme
1192	! of De Rosnay (PhD, figure C.2 page 156), this has been calculated in
1193	! vertical_soil.f90 as zdr.
1194	tsoil_decomp(:) = &
1195	tsoil_decomp(:) + stempdiag(:,ibdl) * rpc * &
1196	( EXP(-zdr(ibdl-1)/z_decomp) - EXP(-zdr(ibdl)/z_decomp) )
1197
1198	ENDDO
1199
1200	control_temp(:,ibelow) = control_temp_func (npts, tsoil_decomp)
1201
1202	!! 4. Moisture control. Factor between 0 and 1
1203	!! 4.1 above the ground: litter humidity
1204	control_moist(:,iabove) = control_moist_func (npts, litterhum)
1205
1206	!! 4.2 below: convolution of humidity and decomposer profiles
1207	! (exponential decomposer profile supposed)
1208
1209	!! 4.2.2 integrate over the nslm levels
1210	soilhum_decomp(:) = zero
1211
1212	DO ibdl = 1, nslm !Loop over soil levels
1213
1214	! Calculate the soil humidity weighted by the assumed presence of
1215	! the decomposers. Presence of the decomposers is prescribes by the
1216	! presence = e^(-soil_depth/z_decomp) where z_decomp described the
1217	! shape of an exponentially decaying function with soil depth. stempdiag
1218	! is descritized along the nodes of the soil layers following the scheme
1219	! of De Rosnay (PhD, figure C.2 page 156), this has been calculated in
1220	! vertical_soil.f90 as zdr
1221	soilhum_decomp(:) = &
1222	soilhum_decomp(:) + shumdiag(:,ibdl) * rpc * &
1223	( EXP(-zdr(ibdl-1)/z_decomp) - EXP(-zdr(ibdl)/z_decomp) )
1224
1225	ENDDO
1226
1227	control_moist(:,ibelow) = control_moist_func (npts, soilhum_decomp)
1228
1229	! Debug
1230	IF(printlev_loc.GE.4) THEN
1231	WRITE(numout,*) 'tsurf', tsurf(test_grid)
1232	WRITE(numout,*) 'literhum', litterhum(test_grid)
1233	WRITE(numout,*) 'soil_hum_decomp', soilhum_decomp(test_grid)
1234	ENDIF
1235
1236	!! 5. fluxes from litter to carbon pools and respiration
1237
1238	!Loop over above and below ground litter
1239	DO ilev = 1, nlevs
1240
1241	DO ivm = 1,nvm
1242	IF ( ok_soil_carbon_discretization ) THEN
1243	itarget=iactive
1244	ELSE
1245	IF (ilev.EQ.iabove)THEN
1246	itarget=isurface
1247	ELSE
1248	itarget=iactive
1249	ENDIF
1250	END IF
1251
1252	!! 5.1 structural litter: goes into active and slow carbon
1253	! pools + respiration
1254
1255	!! 5.1.1 total quantity of structural litter which is decomposed
1256	fd(:) = dtlitter_dec_fac(istructural) &
1257	control_temp(:,ilev) * control_moist(:,ilev) * &
1258	exp( -litter_struct_coef * lignin_struc(:,ivm,ilev) )
1259
1260	! Ensure that fd is between 0 and 1
1261	fd(:) = MAX(zero,MIN(fd(:),un))
1262
1263	!! decompose same fraction of structural part of dead leaves.
1264	! Not exact as lignin content is not the same as that of the total
1265	! structural litter. to avoid a multiple (for ibelow and iabove)
1266	! modification of dead_leaves, we do this test to do this calculate
1267	! only ones in 1,nlev loop
1268	IF (ilev == iabove) THEN
1269
1270	dead_leaves(:,ivm,istructural) = &
1271	dead_leaves(:,ivm,istructural) * ( un - fd(:) )
1272
1273	ENDIF
1274
1275	!! 5.1.2 Calculate the amount of structural litter
1276	! that was decomposed
1277	err(:) = zero
1278	WHERE ((litter(:,istructural,ivm,ilev,icarbon).GE.zero) .AND. &
1279	(litter(:,istructural,ivm,ilev,initrogen).GE.zero))
1280
1281	! Calculate the quantity of structural litter that should
1282	! be decomposed
1283	qd(:,icarbon) = litter(:,istructural,ivm,ilev,icarbon) * fd(:)
1284	qd(:,initrogen) = litter(:,istructural,ivm,ilev,initrogen) * fd(:)
1285
1286	ELSEWHERE
1287
1288	! Keep track of how many errors there are
1289	err(:) = 1.
1290
1291	ENDWHERE
1292
1293	! Check for errors
1294	IF (SUM(err(:)) .GT. zero) THEN
1295
1296	WRITE(numout,*) 'Number of errors found: ,', SUM(err(:))
1297	WRITE(numout,*) 'Negative structural litter pools in stomate_litter'
1298	WRITE(numout,*) 'That does not make any sense'
1299	DO iele = 1,nelements
1300	DO ipts=1,npts
1301	! Notice that checking for a value below zero here will
1302	! not catch a NaN in production mode, even though that
1303	! is registered as an error above.
1304	WRITE(numout,*) 'Structural litter pool, ', &
1305	ipts, ivm, ilev, iele, litter(ipts,istructural,ivm,ilev,iele)
1306	END DO
1307	END DO
1308	CALL ipslerr(3,'stomate_litter','Negative litter pools',&
1309	'This should not be the case','')
1310	!-
1311
1312	END IF
1313
1314
1315	! Subtract the mineralized fraction from the litter but make sure
1316	! the mineralisation cannot become negative
1317	err(:) = zero
1318	WHERE ( (litter(:,istructural,ivm,ilev,icarbon).GE.qd(:,icarbon)) .AND. &
1319	(litter(:,istructural,ivm,ilev,initrogen).GE.qd(:,initrogen)))
1320
1321	! qd is calculated as litter(:,istructural,ivm,ilev,iele)*fd(:)
1322	! as long as fd is positive (which was enforced above) the
1323	! subtraction below should never result in a number that is
1324	! less than zero (except for precision issues). fd is the same
1325	! for carbon and nitrogen and therefore the ratio is preserved
1326	! during mineralization. The test conditions should be
1327	! satisfied for C and N at the same time.
1328	litter(:,istructural,ivm,ilev,icarbon) = MAX(zero, &
1329	litter(:,istructural,ivm,ilev,icarbon) - qd(:,icarbon))
1330	litter(:,istructural,ivm,ilev,initrogen) = MAX(zero, &
1331	litter(:,istructural,ivm,ilev,initrogen) - qd(:,initrogen))
1332
1333	ELSEWHERE
1334
1335	err(:) = 1.
1336
1337	ENDWHERE
1338
1339	! Check for errors
1340	IF (SUM(err(:)) .GT. zero) THEN
1341
1342	WRITE(numout,*) 'Number of errors found: ,', SUM(err(:))
1343	WRITE(numout,*) 'Mineralization exceeds litter pools in stomate_litter'
1344	WRITE(numout,*) 'This will cause problems later on'
1345	DO iele = 1,nelements
1346	WRITE(numout,*) 'Structural litter pool, qd, ', &
1347	ivm, ilev, iele, litter(:,istructural,ivm,ilev,iele),qd(:,iele)
1348	END DO
1349	CALL ipslerr(3,'stomate_litter','Structural mineralization exceeds litter pool',&
1350	'This will cause problems later on','Fix the problem now')
1351
1352	ENDIF
1353
1354	! After large-scale dieback events, so much soil mineral N becomes immobilized
1355	! to decompose litter that too little N is left for plant regrowth. To address
1356	! this, we implicitly represent the action of fungivores, which release N for
1357	! the plants and increase N turnover rates. We set aside a fraction of qd,
1358	! n_fungivores, which becomes available for plant uptake in nitrogen_dynamics.
1359	n_mineralisation(:,ivm) = n_mineralisation(:,ivm) + (1-f_fungivores(:,ivm))*qd(:,initrogen)
1360	n_fungivores(:,ivm) = n_fungivores(:,ivm) + f_fungivores(:,ivm)*qd(:,initrogen)
1361
1362
1363	!! 5.1.3 non-lignin fraction of structural litter goes into
1364	!! active (or surface) carbon pool + respiration
1365	residual(:) = qd(:,icarbon)
1366	som_input(:,itarget,ivm,icarbon) = &
1367	som_input(:,itarget,ivm,icarbon) + &
1368	frac_soil(istructural,itarget,ilev) * qd(:,icarbon) * &
1369	( 1. - lignin_struc(:,ivm,ilev) ) / dt
1370	residual(:) = residual(:) - frac_soil(istructural,itarget,ilev) * qd(:,icarbon) * &
1371	( 1. - lignin_struc(:,ivm,ilev) )
1372
1373	! Here resp_hetero_litter is divided by dt to
1374	! have a value which corresponds to the sechiba time step but
1375	! then in stomate.f90 resp_hetero_litter is multiplied by dt.
1376	! Perhaps it could be simplified.
1377	resp_hetero_litter(:,ivm) = resp_hetero_litter(:,ivm) + &
1378	( 1. - frac_soil(istructural,itarget,ilev) ) * qd(:,icarbon) * &
1379	( 1. - lignin_struc(:,ivm,ilev) ) / dt
1380	residual(:) = residual(:) - ( 1. - frac_soil(istructural,itarget,ilev) ) * qd(:,icarbon) * &
1381	( 1. - lignin_struc(:,ivm,ilev) )
1382
1383	!! 5.1.4 lignin fraction of structural litter goes into
1384	!! slow carbon pool + respiration
1385	som_input(:,islow,ivm,icarbon) = som_input(:,islow,ivm,icarbon) + &
1386	frac_soil(istructural,islow,ilev) * qd(:,icarbon) * &
1387	lignin_struc(:,ivm,ilev) / dt
1388	residual(:) = residual(:) - frac_soil(istructural,islow,ilev) * qd(:,icarbon) * &
1389	lignin_struc(:,ivm,ilev)
1390
1391	! BE CAREFUL: Here resp_hetero_litter is divided by dt to have a
1392	! value which corresponds to the sechiba time step but then in
1393	! stomate.f90 resp_hetero_litter is multiplied by dt. Perhaps it
1394	! could be simplified. Moreover, we must totally adapt the
1395	! routines to the dt_sechiba/one_day time step and avoid some
1396	! constructions that could create bug during future developments.
1397	! The last component of the heterothropic respiration could be
1398	! calculated as below but the residual will be used to avoid
1399	! mass balance problems.
1400	!resp_hetero_litter(:,ivm) = resp_hetero_litter(:,ivm) + &
1401	! ( 1. - frac_soil(istructural,islow,ilev) ) * qd(:,icarbon) * &
1402	! lignin_struc(:,ivm,ilev) / dt
1403	resp_hetero_litter(:,ivm) = resp_hetero_litter(:,ivm) + residual(:) / dt
1404
1405	! Debug
1406	IF(printlev_loc.GE.4 .AND. ivm.EQ.test_pft)THEN
1407	WRITE(numout,*) 'check 1, structural litter ',ilev
1408	WRITE(numout,*) 'dead_leaves(istructural), ',&
1409	dead_leaves(test_grid,test_pft,istructural)
1410	WRITE(numout,*) 'qd', ilev, qd(test_grid,:)
1411	WRITE(numout,*) 'resp_hetero, ',resp_hetero_litter(test_grid,test_pft)
1412	WRITE(numout,*) 'som_input - C, ', &
1413	SUM(som_input(test_grid,:,test_pft,icarbon))
1414	WRITE(numout,*) 'n_mineralisation, ',n_mineralisation(test_grid,test_pft)
1415	WRITE(numout,*) 'som_input - N, ',&
1416	SUM(som_input(test_grid,:,test_pft,initrogen))
1417	ENDIF
1418	!-
1419
1420	!! 5.2 metabolic litter goes into active carbon pool + respiration
1421	!! 5.2.1 total quantity of metabolic litter that is decomposed
1422	fd(:) = dtlitter_dec_fac(imetabolic) control_temp(:,ilev) * &
1423	control_moist(:,ilev)
1424
1425	! Ensure that fd is between 0 and 1
1426	fd(:) = MAX(zero,MIN(fd(:),un))
1427
1428	! Calculate qd, make sure it has a positive value
1429	err(:) = zero
1430	WHERE ((litter(:,imetabolic,ivm,ilev,icarbon).GE.zero) .AND. &
1431	(litter(:,imetabolic,ivm,ilev,initrogen).GE.zero))
1432
1433	! Calculate the quantity of metabolic litter that should
1434	! be decomposed
1435	qd(:,icarbon) = litter(:,imetabolic,ivm,ilev,icarbon) * fd(:)
1436
1437	! Accumulate qd(initrogen) for later use
1438	qd(:,initrogen) = litter(:,imetabolic,ivm,ilev,initrogen) * fd(:)
1439
1440	ELSEWHERE
1441
1442	! Keep track of how many errors there are
1443	err(:) = 1.
1444
1445	ENDWHERE
1446
1447	! Check for errors
1448	IF (SUM(err(:)) .GT. zero) THEN
1449
1450	WRITE(numout,*) 'Number of errors found in icarbon: ,', SUM(err(:))
1451	WRITE(numout,*) 'Negative litter pools in stomate_litter'
1452	WRITE(numout,*) 'That does not make any sense'
1453	DO iele = 1,nelements
1454	DO ipts=1,npts
1455	! Notice that checking for a value below zero here will
1456	! not catch a NaN in production mode, even though that
1457	! is registered as an error above.
1458	WRITE(numout,*) 'Metabolic litter pool, ', &
1459	ipts,ivm, ilev, iele, litter(ipts,imetabolic,ivm,ilev,iele)
1460	ENDDO
1461	END DO
1462	CALL ipslerr(3,'stomate_litter','Negative metabolic litter pools',&
1463	'This should not be the case','')
1464	!-
1465
1466	END IF
1467
1468	! Subtract the mineralization from the litter pool but make
1469	! sure that no negative values for the litter pool are created
1470	err(:) = zero
1471	WHERE ( (litter(:,imetabolic,ivm,ilev,icarbon).GE.qd(:,icarbon)) .AND. &
1472	(litter(:,imetabolic,ivm,ilev,initrogen).GE.qd(:,initrogen)))
1473
1474	! qd is calculated as litter(:,imetabolic,ivm,ilev,iele)*fd(:)
1475	! as long as fd is positive (which was enforced above) the
1476	! subtraction below should never result in a number that is
1477	! less than zero (except for precision issues). fd is the same
1478	! for carbon and nitrogen and therefore the ratio is preserved
1479	! during mineralization. The test conditions should be
1480	! satisfied for C and N at the same time.
1481	litter(:,imetabolic,ivm,ilev,icarbon) = MAX(zero, &
1482	litter(:,imetabolic,ivm,ilev,icarbon) - qd(:,icarbon))
1483	litter(:,imetabolic,ivm,ilev,initrogen) = MAX(zero, &
1484	litter(:,imetabolic,ivm,ilev,initrogen) - qd(:,initrogen))
1485
1486	ELSEWHERE
1487
1488	err(:) = 1.
1489
1490	ENDWHERE
1491
1492	! Check for errors
1493	IF (SUM(err(:)) .GT. zero) THEN
1494
1495	WRITE(numout,*) 'Number of errors found: ,', SUM(err(:))
1496	WRITE(numout,*) 'Mineralization exceeds litter pools in stomate_litter'
1497	WRITE(numout,*) 'This will cause problems latter on'
1498	DO iele = 1,nelements
1499	WRITE(numout,*) 'Metabolic litter pool, qd, ', &
1500	ivm, ilev, iele, litter(:,imetabolic,ivm,ilev,iele),qd(:,initrogen)
1501	END DO
1502	CALL ipslerr(3,'stomate_litter','Metabolic mineralization exceeds litter pool',&
1503	'This will cause problems later on','Fix the problem now')
1504	!-
1505
1506	ENDIF
1507
1508	!! 5.2.2 decompose same fraction of metabolic part of dead leaves.
1509	! to avoid a multiple (for ibelow and iabove) modification of
1510	! dead_leaves, we do this test to do this calcul only ones in
1511	! 1,nlev loop
1512	IF (ilev == iabove) THEN
1513
1514	dead_leaves(:,ivm,imetabolic) = &
1515	dead_leaves(:,ivm,imetabolic) * ( 1. - fd(:) )
1516
1517	ENDIF
1518
1519	! After large-scale dieback events, so much soil mineral N becomes immobilized
1520	! to decompose litter that too little N is left for plant regrowth. To address
1521	! this, we implicitly represent the action of fungivores, which release N for
1522	! the plants and increase N turnover rates. We set aside a fraction of qd,
1523	! n_fungivores, which becomes available for plant uptake in nitrogen_dynamics.
1524	n_mineralisation(:,ivm) = n_mineralisation(:,ivm) + (1-f_fungivores(:,ivm))*qd(:,initrogen)
1525	n_fungivores(:,ivm) = n_fungivores(:,ivm) + f_fungivores(:,ivm)*qd(:,initrogen)
1526
1527	!! 5.2.3 put decomposed litter into active
1528	! (or surface) pool + respiration
1529	residual(:) = qd(:,icarbon)
1530	som_input(:,itarget,ivm,icarbon) = &
1531	som_input(:,itarget,ivm,icarbon) + &
1532	frac_soil(imetabolic,itarget,ilev) * qd(:,icarbon) / dt
1533	residual(:) = residual(:) - frac_soil(imetabolic,itarget,ilev) * qd(:,icarbon)
1534
1535	! BE CAREFUL: Here resp_hetero_litter is divided by dt to have a
1536	! value which corresponds to the sechiba time step but then in
1537	! stomate.f90 resp_hetero_litter is multiplied by dt.
1538	! Perhaps it could be simplified. Moreover, we must totally adapt
1539	! the routines to the dtradia/one_day time step and avoid some
1540	! constructions that could create bugs during future developments.
1541	! The last component of the heterothropic respiration could be
1542	! calculated as below but the residual will be used to avoid
1543	! mass balance problems.
1544	!resp_hetero_litter(:,ivm) = resp_hetero_litter(:,ivm) + &
1545	! ( 1. - frac_soil(imetabolic,itarget,ilev) ) * qd(:,icarbon) / dt
1546	resp_hetero_litter(:,ivm) = resp_hetero_litter(:,ivm) + residual(:) / dt
1547
1548	! Debug
1549	IF(printlev_loc.GE.4 .AND. ivm.EQ.test_pft)THEN
1550	WRITE(numout,*) 'check 2, metabolic litter',ilev
1551	WRITE(numout,*) 'qd, ilev', ilev, qd(test_grid,:)
1552	WRITE(numout,*) 'litter', litter(test_grid,imetabolic,ivm,ilev,:)
1553	WRITE(numout,*) 'resp_hetero, ',resp_hetero_litter(test_grid,test_pft)
1554	WRITE(numout,*) 'som_input - C, ', &
1555	SUM(som_input(test_grid,:,test_pft,icarbon))
1556	WRITE(numout,*) 'n_mineralisation, ',n_mineralisation(test_grid,test_pft)
1557	WRITE(numout,*) 'som_input - N, ',&
1558	SUM(som_input(test_grid,:,test_pft,initrogen))
1559	WRITE(numout,*) ''
1560	ENDIF
1561	!-
1562
1563	!! 5.3 woody litter: goes into active and slow carbon
1564	! pools + respiration
1565
1566	!! 5.3.1 total quantity of woody litter which is decomposed
1567	fd(:) = dtlitter_dec_fac(iwoody) &
1568	control_temp(:,ilev) * control_moist(:,ilev) * &
1569	EXP( -3. * lignin_wood(:,ivm,ilev) )
1570
1571	! Ensure that fd is between 0 and 1
1572	fd(:) = MAX(zero,MIN(fd(:),un))
1573
1574	! Calculate qd for woody litter
1575	err(:) = zero
1576	WHERE ((litter(:,iwoody,ivm,ilev,icarbon).GE.zero) .AND. &
1577	(litter(:,iwoody,ivm,ilev,initrogen).GE.zero))
1578
1579	! Calculate the quantity of woody litter that should
1580	! be decomposed
1581	qd(:,icarbon) = litter(:,iwoody,ivm,ilev,icarbon) * fd(:)
1582	qd(:,initrogen) = litter(:,iwoody,ivm,ilev,initrogen) * fd(:)
1583
1584	ELSEWHERE
1585
1586	! Keep track of how many errors there are
1587	err(:) = 1.
1588
1589	ENDWHERE
1590
1591	! Check for errors
1592	IF (SUM(err(:)) .GT. zero) THEN
1593
1594	WRITE(numout,*) 'Number of errors found: ,', SUM(err(:))
1595	WRITE(numout,*) 'Negative litter pools in stomate_litter'
1596	WRITE(numout,*) 'That does not make any sense'
1597	DO iele = 1,nelements
1598	WRITE(numout,*) 'Woody litter pool, ', &
1599	ivm, ilev, iele, litter(:,iwoody,ivm,ilev,iele)
1600	END DO
1601	CALL ipslerr(3,'stomate_litter','Negative woody litter pools',&
1602	'This should not be the case','')
1603	!-
1604
1605	END IF
1606
1607	! Subtract the mineralization from the litter pool but make
1608	! sure that no negative values for the litter pool are created
1609	err(:) = zero
1610	WHERE ( (litter(:,iwoody,ivm,ilev,icarbon).GE.qd(:,icarbon)) .AND. &
1611	(litter(:,iwoody,ivm,ilev,initrogen).GE.qd(:,initrogen)))
1612
1613	! qd is calculated as litter(:,iwoody,ivm,ilev,iele)*fd(:)
1614	! as long as fd is positive (which was enforced above) the
1615	! subtraction below should never result in a number that is
1616	! less than zero (except for precision issues). fd is the same
1617	! for carbon and nitrogen and therefore the ratio is preserved
1618	! during mineralization. The test conditions should be
1619	! satisfied for C and N at the same time.
1620	litter(:,iwoody,ivm,ilev,icarbon) = MAX(zero, &
1621	litter(:,iwoody,ivm,ilev,icarbon) - qd(:,icarbon))
1622	litter(:,iwoody,ivm,ilev,initrogen) = MAX(zero, &
1623	litter(:,iwoody,ivm,ilev,initrogen) - qd(:,initrogen))
1624
1625	ELSEWHERE
1626
1627	err(:) = 1.
1628
1629	ENDWHERE
1630
1631	! Check for errors
1632	IF (SUM(err(:)) .GT. zero) THEN
1633
1634	WRITE(numout,*) 'Number of errors found: ,', SUM(err(:))
1635	WRITE(numout,*) 'Mineralization exceeds litter pools in stomate_litter'
1636	WRITE(numout,*) 'This will cause problems latter on'
1637	DO iele = 1,nelements
1638	WRITE(numout,*) 'Woody litter pool, qd, ', &
1639	ivm, ilev, iele, litter(:,iwoody,ivm,ilev,iele),qd(:,initrogen)
1640	END DO
1641	CALL ipslerr(3,'stomate_litter','Woody mineralization exceeds litter pool',&
1642	'This will cause problems later on','Fix the problem now')
1643	!-
1644
1645	ENDIF
1646
1647	! After large-scale dieback events, so much soil mineral N becomes immobilized
1648	! to decompose litter that too little N is left for plant regrowth. To address
1649	! this, we implicitly represent the action of fungivores, which release N for
1650	! the plants and increase N turnover rates. We set aside a fraction of qd,
1651	! n_fungivores, which becomes available for plant uptake in nitrogen_dynamics.
1652	n_mineralisation(:,ivm) = n_mineralisation(:,ivm) + (1-f_fungivores(:,ivm))*qd(:,initrogen)
1653	n_fungivores(:,ivm) = n_fungivores(:,ivm) + f_fungivores(:,ivm)*qd(:,initrogen)
1654
1655	!! 5.3.2 non-lignin fraction of woody litter goes into
1656	!! active/structural carbon pool + respiration (per time unit)
1657	residual(:) = qd(:,icarbon)
1658	som_input(:,itarget,ivm,icarbon) = som_input(:,itarget,ivm,icarbon) + &
1659	frac_soil(iwoody,itarget,ilev) * qd(:,icarbon) * &
1660	( 1. - lignin_wood(:,ivm,ilev) ) / dt
1661	residual(:) = residual(:) - frac_soil(iwoody,itarget,ilev) * qd(:,icarbon) * &
1662	( 1. - lignin_wood(:,ivm,ilev) )
1663
1664	resp_hetero_litter(:,ivm) = resp_hetero_litter(:,ivm) + &
1665	( 1. - frac_soil(iwoody,itarget,ilev) ) * qd(:,icarbon) * &
1666	( 1. - lignin_wood(:,ivm,ilev) ) / dt
1667	residual(:) = residual(:) - ( 1. - frac_soil(iwoody,itarget,ilev) ) * qd(:,icarbon) * &
1668	( 1. - lignin_wood(:,ivm,ilev) )
1669
1670	!! 5.3.3 lignin fraction of woody litter goes into
1671	!! slow carbon pool + respiration (per time unit)
1672	som_input(:,islow,ivm,icarbon) = som_input(:,islow,ivm,icarbon) + &
1673	frac_soil(iwoody,islow,ilev) * qd(:,icarbon) * &
1674	lignin_wood(:,ivm,ilev) / dt
1675	residual(:) = residual(:) - frac_soil(iwoody,islow,ilev) * qd(:,icarbon) * &
1676	lignin_wood(:,ivm,ilev)
1677
1678	! The last component of the heterothropic respiration could be
1679	! calculated as below but the residual will be used to avoid
1680	! mass balance problems.
1681	resp_hetero_litter(:,ivm) = resp_hetero_litter(:,ivm) + &
1682	( 1. - frac_soil(iwoody,islow,ilev) ) * qd(:,icarbon) * &
1683	lignin_wood(:,ivm,ilev) / dt
1684	resp_hetero_litter(:,ivm) = resp_hetero_litter(:,ivm) + residual(:) / dt
1685
1686	! Debug
1687	IF(printlev_loc.GE.4 .AND. ivm.EQ.test_pft)THEN
1688	WRITE(numout,*) 'check 3, ',ilev
1689	WRITE(numout,*) 'qd, icarbon, ilev', ilev, qd(test_grid,icarbon)
1690	WRITE(numout,*) 'resp_hetero, ',resp_hetero_litter(test_grid,test_pft)
1691	WRITE(numout,*) 'som_input - C, ', &
1692	SUM(som_input(test_grid,:,test_pft,icarbon))
1693	WRITE(numout,*) 'n_mineralisation, ',n_mineralisation(test_grid,test_pft)
1694	WRITE(numout,*) 'som_input - N, ',&
1695	SUM(som_input(test_grid,:,test_pft,initrogen))
1696	WRITE(numout,*) ''
1697	ENDIF
1698	!-
1699
1700	ENDDO ! nvm
1701
1702	ENDDO ! nlevs
1703
1704	! 5.4 Calculate som_input for nitrogen
1705	! This far the som input for icarbon has been recorded. Now
1706	! use the CN ratios to calculate the som input for nitrogen.
1707	! Note the difference in units between ::som_input (gN m-2 s-1)
1708	! and ::n_mineralisation, (gN m-2).
1709	som_input(:,iactive,:,initrogen) = &
1710	som_input(:,iactive,:,icarbon)/CN_target(:,:,iactive)
1711	som_input(:,islow,:,initrogen) = &
1712	som_input(:,islow,:,icarbon)/CN_target(:,:,islow)
1713	som_input(:,isurface,:,initrogen) = &
1714	som_input(:,isurface,:,icarbon)/CN_target(:,:,isurface)
1715
1716	! Ideally we take the nitrogen contained in som_input from the
1717	! n_mineralisation pool. Before we can do so we have to check whether
1718	! there is enough nitrogen to do so. If not we will have to decompose
1719	! less litter than we would like to. That is what is being checked
1720	! in the rest of this subroutine. If we have to adjust som_input, we
1721	! may need the following factors later:
1722	DO ipts = 1, npts
1723
1724	DO ivm = 1, nvm
1725
1726	n_min_old(ipts,ivm) = n_mineralisation(ipts,ivm)
1727	som_input_total(ivm,initrogen) = &
1728	som_input(ipts,iactive,ivm,initrogen) + &
1729	som_input(ipts,islow,ivm,initrogen) + &
1730	som_input(ipts,isurface,ivm,initrogen)
1731
1732	IF (som_input_total(ivm,initrogen).GT.zero) THEN
1733
1734	share_som_input_act_n(ipts,ivm) = &
1735	som_input(ipts,iactive,ivm,initrogen) &
1736	/ som_input_total(ivm,initrogen)
1737	share_som_input_slow_n(ipts,ivm) = &
1738	som_input(ipts,islow,ivm,initrogen) &
1739	/ som_input_total(ivm,initrogen)
1740	share_som_input_surf_n(ipts,ivm) = &
1741	som_input(ipts,isurface,ivm,initrogen) &
1742	/ som_input_total(ivm,initrogen)
1743	ENDIF
1744
1745	ENDDO
1746
1747	ENDDO
1748
1749	!! 5.5 Forward checking of n_mineralization
1750	! Up to here we calculated the som_input from the litter point of view but
1751	! decomposing so much litter may require more nitrogen than currently
1752	! available. If that is the case we will truncate som_input(:,:,:,initrogen)
1753	! to a value that will not result in conflicts in stomate_soilcarbon.
1754
1755	! 5.5.1 We first need to pre-determine the value of n_mineralisation(:,:)
1756	! we will have in the next routine (som_dynamics) after the fluxes are
1757	! subtracted. This means we have to port some of the code from som_dynamics
1758	! to here. The following block of code is copied from som_dynamics:
1759
1760	!! Initialise
1761	resp_hetero_soil(:,:) = zero
1762	total_som_n_flux(:,:,:) = zero
1763	n_mineralisation_som(:,:) = zero
1764	frac_carb(:,:,:) = zero
1765	som_temp(:,:,:,:) = zero
1766
1767	! Get soil "constants"
1768	! Flux fractions between carbon pools: depend on clay content,
1769	! recalculated each time
1770	! From active pool: depends on clay content
1771	frac_carb(:,iactive,ipassive) = active_to_pass_ref_frac + &
1772	active_to_pass_clay_frac*clay(:)
1773	frac_carb(:,iactive,islow) = un - frac_carb(:,iactive,ipassive) - &
1774	(active_to_co2_ref_frac - active_to_co2_clay_silt_frac*(clay(:)+silt(:)))
1775
1776	! From slow pool
1777	frac_carb(:,islow,ipassive) = slow_to_pass_ref_frac + &
1778	slow_to_pass_clay_frac*clay(:)
1779
1780	! OCN doesn't use Parton 1993 formulation for
1781	! frac_carb(:,islow,ipassive)
1782	! but the one of 1987 : ie = 0.03 ....
1783	frac_carb(:,islow,iactive) = un - frac_carb(:,islow,ipassive) - &
1784	slow_to_co2_ref_frac
1785
1786	! From passive pool
1787	frac_carb(:,ipassive,iactive) = pass_to_active_ref_frac
1788	frac_carb(:,ipassive,islow) = pass_to_slow_ref_frac
1789
1790	! From surface pool
1791	frac_carb(:,isurface,islow) = surf_to_slow_ref_frac
1792
1793	!! Determine the respiration fraction : what's left after
1794	! subtracting all the 'pool-to-pool' flux fractions
1795	! Diagonal elements of frac_carb are zero
1796	frac_resp(:,:) = un - frac_carb(:,:,isurface) - frac_carb(:,:,iactive) - &
1797	frac_carb(:,:,islow) - frac_carb(:,:,ipassive)
1798
1799	!! Turnover in SOM pools (in days)
1800	! som_turn_ipool is the turnover (in years)
1801	! It is weighted by Temp and Humidity function later
1802	som_turn(iactive) = som_turn_iactive / one_year
1803	som_turn(islow) = som_turn_islow / one_year
1804	som_turn(ipassive) = som_turn_ipassive / one_year
1805	som_turn(isurface) = som_turn_isurface / one_year
1806
1807	! Update the SOM stocks with the different soil carbon input that
1808	! we would like to add from litter decomposition.
1809	som_temp(:,:,:,:) = som(:,:,:,:) + som_input(:,:,:,:) * dt
1810
1811	! Fluxes between carbon reservoirs, and to the atmosphere (respiration)
1812	! Calculate fluxes out of pools
1813	! Loop over carbon pools from which the flux comes
1814	DO ivm = 1,nvm
1815
1816	DO icarb = 1, ncarb
1817
1818	DO iele = 1, nelements
1819
1820	! Determine total flux out of pool [gN m-2]
1821	fluxtot(:,ivm,icarb,iele) = dtsom_turn(icarb) &
1822	som_temp(:,icarb,ivm,iele) * control_moist(:,ibelow) * &
1823	control_temp(:,ibelow) * decomp_factor(ivm)
1824
1825	! Flux from active pools depends on clay content [gN m-2]
1826	IF ( icarb .EQ. iactive ) THEN
1827	fluxtot(:,ivm,icarb,iele) = fluxtot(:,ivm,icarb,iele) * &
1828	( un - som_turn_iactive_clay_frac * clay(:) )
1829	ENDIF
1830
1831	! Update the loss in each carbon pool [gN m-2]
1832	som_temp(:,icarb,ivm,iele) = som_temp(:,icarb,ivm,iele) - &
1833	fluxtot(:,ivm,icarb,iele)
1834
1835	ENDDO ! nelements
1836
1837	! Fluxes towards the other pools (icarb -> iicarb)
1838	! Loop over the SOM pools where the flux goes
1839	DO iicarb = 1, ncarb
1840
1841	! Carbon flux [gN m-2]
1842	flux(:,ivm,icarb,iicarb,icarbon) = frac_carb(:,icarb,iicarb) * &
1843	fluxtot(:,ivm,icarb,icarbon)
1844
1845	! Nitrogen flux - Function of the C stock of the 'departure'
1846	! pool and of the C to N target ratio of the 'arrival' pool [gN
1847	! m-2]
1848	flux(:,ivm,icarb,iicarb,initrogen) = frac_carb(:,icarb,iicarb) * &
1849	fluxtot(:,ivm,icarb,icarbon) / &
1850	CN_target(:,ivm,iicarb)
1851
1852	ENDDO ! receiving SOM pools
1853
1854	ENDDO ! donating SOM pools
1855
1856	!! Respiration
1857	! Fluxtot is in [gN m-2], resp_hetero_soil is in [gN m-2 dt-1]
1858	! divide by dt.
1859	resp_hetero_soil(:,ivm) = &
1860	( frac_resp(:,iactive) * fluxtot(:,ivm,iactive,icarbon) + &
1861	frac_resp(:,islow) * fluxtot(:,ivm,islow,icarbon) + &
1862	frac_resp(:,ipassive) * fluxtot(:,ivm,ipassive,icarbon) + &
1863	frac_resp(:,isurface) * fluxtot(:,ivm,isurface,icarbon) ) / dt
1864
1865	!! Add fluxes to active, slow, and passive pools
1866	! Loop over SOM pools
1867	DO icarb = 1, ncarb
1868
1869	! Calculate the components of the som fluxes [gN m-2]
1870	! or [gC m-2]
1871	som_temp(:,icarb,ivm,initrogen) = som_temp(:,icarb,ivm,initrogen) + &
1872	flux(:,ivm,iactive,icarb,initrogen) + &
1873	flux(:,ivm,ipassive,icarb,initrogen) + &
1874	flux(:,ivm,islow,icarb,initrogen) + &
1875	flux(:,ivm,isurface,icarb,initrogen)
1876
1877	! Calculate the total flux [gN m-2]
1878	total_som_n_flux(:,ivm,initrogen) = total_som_n_flux(:,ivm,initrogen) + &
1879	flux(:,ivm,iactive,icarb,initrogen) + &
1880	flux(:,ivm,ipassive,icarb,initrogen) + &
1881	flux(:,ivm,islow,icarb,initrogen) + &
1882	flux(:,ivm,isurface,icarb,initrogen)
1883
1884	ENDDO ! Loop over SOM pools
1885
1886	ENDDO ! End loop over PFTs
1887
1888
1889	! At this point we calculated the fluxes that will occur if we use
1890	! the som_input calculated above in stomate_litter. We will now check
1891	! whether this will result in conflicts. If it results in conflicts
1892	! som_input will be adjusted. if it does not result in conflicts we
1893	! keep som_input as calculated above.
1894
1895	! Calculate n_mineralisation for som_input
1896	DO ivm = 1, nvm
1897
1898	DO ipts = 1, npts
1899
1900	! No redistrubtion needed unless specified in the code
1901	ld_redistribute(ivm) = zero
1902
1903	! Only check for existing PFTs
1904	IF (veget_max(ipts,ivm).LE.zero) CYCLE
1905
1906	! To decompose the som + som_input from the litter we
1907	! need to take nitrogen from the mineralised pool. Note that
1908	! litter decomposition does not contribute directly to the
1909	! passive component of the soil organic matter [gN m-2 dt-1]
1910	! multiply with dt to obtain [gN m-2].
1911	som_input_total(ivm,:) = (som_input(ipts,iactive,ivm,:) + &
1912	som_input(ipts,islow,ivm,:) + &
1913	som_input(ipts,isurface,ivm,:))*dt
1914
1915	! Store the positive component of n_mineralisation calculated
1916	! above. This n_mineralisation is consistent with som_input.
1917	! This is the mineralisation we need to decompose the amount
1918	! of litter initialily estimated in stomate_litter
1919	n_mineralisation_init = n_mineralisation(ipts,ivm)
1920
1921	! Calculate the n_mineralisation pool after all internal
1922	! n fluxes from som_input(initrogen) are accounted for. This
1923	! looks like a worst case scenario.
1924	! n_mineralisation can be both positive or negative
1925	n_mineralisation(ipts,ivm) = n_mineralisation_init - &
1926	som_input_total(ivm,initrogen)
1927
1928	! Check whether the decomposition estimates above are realistic in
1929	! terms of the nitrogen available for immobilisation if
1930	! n_mineralisation is negative.
1931	IF ((soil_n_min(ipts,ivm,initrate)+soil_n_min(ipts,ivm,iammonium)) + &
1932	n_mineralisation_init.LT.zero) THEN
1933
1934	CALL ipslerr(3,'stomate_litter',&
1935	'First time that this problem is encountered',&
1936	'Debug the code first before continuing','Really!')
1937
1938	!+++CHECK+++
1939	! Possible solution but it has never been checked
1940	! All the C and N contained in som_input will need to be moved
1941	! back into the litter pool. Furthermore the n_mineralisation
1942	! needs to be decreased as is the heterotrophic respiration
1943	DO iele = 1,nelements
1944	litter_total_new(iele) = &
1945	SUM(SUM(litter(ipts,:,ivm,:,iele),1)) + som_input_total(ivm,iele)
1946	som_input_total(ivm,iele) = zero
1947	ENDDO
1948
1949	! Adjust resp_hetero_litter proportionaly.
1950	delta_hetero_litter = (resp_hetero_litter(ipts,ivm)- &
1951	(resp_hetero_litter(ipts,ivm)/n_mineralisation_init) * &
1952	(soil_n_min(ipts,ivm,initrate)+soil_n_min(ipts,ivm,iammonium)))*dt
1953
1954	! Update the litter pools
1955	litter_total_new(initrogen) = litter_total_new(initrogen) + &
1956	n_mineralisation_init - (soil_n_min(ipts,ivm,initrate)+soil_n_min(ipts,ivm,iammonium))
1957	litter_total_new(icarbon) = litter_total_new(icarbon) + delta_hetero_litter
1958
1959	! Set flag to redistribute the new pools
1960	ld_redistribute(ivm) = 1.
1961	!+++++++++++
1962
1963	ELSEIF ((((soil_n_min(ipts,ivm,initrate)+soil_n_min(ipts,ivm,iammonium)) + &
1964	n_mineralisation(ipts,ivm).LT.zero) .AND. &
1965	((soil_n_min(ipts,ivm,initrate)+soil_n_min(ipts,ivm,iammonium)) + &
1966	n_mineralisation_init.GT.zero)))THEN
1967
1968	! If the above conditions are satified there exists a solution.
1969	! The new n_mineralisation should get a value between its initial
1970	! and its current value.
1971
1972	! Initialize
1973	litter_total_new(:) = zero
1974
1975	! If we have persistent negative values, however, we
1976	! can reach a situation where we completely deplete our soil_n_min
1977	! pools to satisfy the demand for immobilization. To remedy this,
1978	! we'll truncate som_input based on the current
1979	! n_mineralisation pool, but we'll limit the occurrence of negative
1980	! values. Instead we'll set a minimum (min_n) for n_mineralisation
1981	! (e.g. .00001 g/m^2) so we can keep some N in the soil.
1982	IF(n_min_old(ipts,ivm).GT.min_n) THEN
1983
1984	! Later in this code n_mineralisation will be calculated as
1985	! n_mineralisation(ipts,ivm) = n_min_old(ipts,ivm) - som_input_new(ivm,initrogen)
1986	! So, if the n_min_old is GT then min_n we will truncate the
1987	! n_mineralisation at min_n. This requires that som_input_new is
1988	! calculated as below.
1989	som_input_new(ivm,initrogen) = n_min_old(ipts,ivm) - &
1990	min_n
1991
1992	ELSE
1993
1994	! If n_mineralisation is less than min_n, take a fraction of the
1995	! remaining n_mineralisation for som_input. This is arbitrarily
1996	! set to 50% for now. We are just trying to avoid getting values
1997	! that equal to zero.
1998	som_input_new(ivm,initrogen) = 0.5*n_min_old(ipts,ivm)
1999
2000	ENDIF
2001
2002	! Error checking
2003	IF (som_input_new(ivm,initrogen).LT.zero) THEN
2004
2005	WRITE(numout,*) 'ERROR: should not be here ',&
2006	(soil_n_min(ipts,ivm,initrate)+soil_n_min(ipts,ivm,iammonium)),&
2007	ivm, n_mineralisation(ipts,ivm),&
2008	som_input_total(ivm,initrogen)
2009	CALL ipslerr(3,'stomate_litter',&
2010	'Seems that the nitrogen can be satisfied',&
2011	'Should not have ended up in this part of the code','')
2012
2013	ENDIF
2014
2015	! The adjusted som_input carbon [gC m-2] should be
2016	! som_input_new(ivm,icarbon) = som_input_new(ivm,initrogen) / CN_input_total
2017	! However, When too little N is available to reach the target C:N ratio of
2018	! the receiving pool, the equation above results in values that are too
2019	! low for som_input_total(:,icarbon), meaning carbon litter can not
2020	! decompose. som_input_total(:,icarbon) still needs to be constrainted,
2021	! however, or else the C:N ratio in the surface soil pool will become
2022	! unrealistic.
2023	som_input_new(ivm,icarbon) = som_input_total(ivm,icarbon)
2024
2025	! Protect against dividing by zero. som_input_new(ivm,initrogen) may be zero here in
2026	! cases where n_mineralisation is zero at the beginning of the run, and then there
2027	! is some new litter input which forces us to this point in the code. This can
2028	! happen if there is no litter for a long time, but to date we have only seen
2029	! this in cases where the GPP is zero for a long time. So I will flag it as a
2030	! warning.
2031	IF(som_input_new(ivm,initrogen) .NE. zero)THEN
2032	IF((som_input_new(ivm,icarbon)/som_input_new(ivm,initrogen)) .GT. max_cn) THEN
2033	som_input_new(ivm,icarbon) = som_input_new(ivm,initrogen) * max_cn
2034	ENDIF
2035	ELSE
2036	! leave the new soil organic matter carbon input as it was, but
2037	! write a warning. This warning is not very useful for crops
2038	! so it is only written for trees.
2039	IF (natural(ivm)) THEN
2040	WRITE(numout,*) 'WARNING: ivm,,ipts ',ivm,ipts
2041	WRITE(numout,*) 'WARNING: som_input_new(ivm,initrogen), &
2042	& n_min_old(ipts,ivm)',som_input_new(ivm,initrogen),&
2043	n_min_old(ipts,ivm)
2044	WRITE(numout,*) 'WARNING: min_n,&
2045	n_mineralisation(ipts,ivm),som_input_new(ivm,icarbon)',&
2046	min_n, n_mineralisation(ipts,ivm),&
2047	som_input_new(ivm,icarbon)
2048	CALL ipslerr(2,'stomate_litter','Perhaps no n_mineralization &
2049	& pools, but new litter','Previously only seen in cases &
2050	& where GPP is zero for long term','')
2051	ENDIF !natural
2052	ENDIF
2053
2054	! Calculate the n_mineralisation pool after all internal
2055	! n fluxes from som_input(initrogen) are accounted for.
2056	! n_mineralisation can be both positive or negative
2057	n_mineralisation(ipts,ivm) = n_min_old(ipts,ivm) - &
2058	som_input_new(ivm,initrogen)
2059
2060	! The nitrogen that was not mineralized should be added back
2061	! into the litter pool. The carbon no longer included in
2062	! som_input should also be added back into the litter pool.
2063	litter_total_new(initrogen) = total_init_nitrogen(ipts,ivm) - &
2064	n_mineralisation(ipts,ivm) - som_input_new(ivm,initrogen) - &
2065	n_fungivores(ipts,ivm)
2066
2067	litter_total_new(icarbon) = SUM(SUM(litter(ipts,:,ivm,:,icarbon),1)) + &
2068	som_input_total(ivm,icarbon) - som_input_new(ivm,icarbon)
2069
2070	! Som_input was adjusted to ensure that after accounting for
2071	! n_mineralisation there is enough soil nitrogen left to
2072	! account for nitrogen immobilisation. When litter is
2073	! transformed into som_input, part of the carbon is lost as
2074	! heterotrophic respiration. For each unit of som_input a
2075	! certain number given by ::resp_ratio is respired
2076	IF ((som_input_total(ivm,icarbon) .GT. zero) .AND. &
2077	(resp_hetero_litter(ipts,ivm).GT.zero)) THEN
2078
2079	! The change will be used to adjust the litter pool. The units should
2080	! be [gC m-2]
2081	delta_hetero_litter = &
2082	(som_input_total(ivm,icarbon) - som_input_new(ivm,icarbon)) * &
2083	(resp_hetero_litter(ipts,ivm)/som_input_total(ivm,icarbon))*dt
2084
2085	! Respiration is a flux, the units should be [gC m-2 dt-1]
2086	resp_hetero_litter(ipts,ivm) = som_input_new(ivm,icarbon) * &
2087	resp_hetero_litter(ipts,ivm) / som_input_total(ivm,icarbon)
2088
2089	ELSE
2090
2091	WRITE(numout,*) 'ERROR: surprise when recalculating &
2092	& resp_hetero_litter, ', som_input_new(ivm,icarbon), &
2093	resp_hetero_litter(ipts,ivm), som_input_total(ivm,icarbon)
2094	CALL ipslerr(3,'stomate_litter','problems with the sign',&
2095	'when calculating delta_hetero_litter','')
2096
2097	ENDIF
2098
2099	! The carbon that is not respired stays in the litter pool.
2100	litter_total_new(icarbon) = litter_total_new(icarbon) + &
2101	delta_hetero_litter
2102
2103	! Set flag to redistribute the new pools
2104	ld_redistribute(ivm) = 1.
2105
2106
2107	ELSEIF ((soil_n_min(ipts,ivm,initrate)+soil_n_min(ipts,ivm,iammonium)) + &
2108	n_mineralisation(ipts,ivm).GE.zero) THEN
2109
2110	! Set flag to redistribute the new pools
2111	ld_redistribute(ivm) = 0
2112
2113	ELSE
2114
2115	WRITE(numout,*) 'ERROR: unexpected case'
2116	WRITE(numout,*) 'mineralisation wanted, ',n_mineralisation(ipts,ivm)
2117	WRITE(numout,*) 'mineralisation_old, ',n_mineralisation_init
2118	WRITE(numout,*) 'soil_n_min, ',(soil_n_min(ipts,ivm,initrate)+soil_n_min(ipts,ivm,iammonium))
2119	CALL ipslerr(3,'stomate_litter','Unexpected case',&
2120	'Check the logic of IF-statement','')
2121
2122	ENDIF
2123
2124	IF (ld_redistribute(ivm) == 1.) THEN
2125
2126	! We are in, so reset the flag
2127	ld_redistribute(ivm) = 0
2128
2129	! Share of each component to the total flux
2130	DO iele = 1,nelements
2131
2132	IF (som_input_total(ivm,iele).GT.zero) THEN
2133
2134	share_som_input_active = som_input(ipts,iactive,ivm,iele) / &
2135	(som_input_total(ivm,iele) / dt)
2136	share_som_input_slow = som_input(ipts,islow,ivm,iele) / &
2137	(som_input_total(ivm,iele) / dt)
2138	share_som_input_surface = MAX(1 - share_som_input_active - &
2139	share_som_input_slow, zero)
2140
2141	ELSE
2142
2143	WRITE(numout,*) 'ERROR: new problem. som_input = zero &
2144	& cannot redistribute the overspended nitrogen'
2145	CALL ipslerr(3,'stomate_litter','som_input = zero',&
2146	'cannot redistribute the overspended N or C','')
2147
2148	ENDIF
2149
2150	! Error checking
2151	IF(err_act.GT.1)THEN
2152
2153	IF ( abs(share_som_input_surface + share_som_input_active + &
2154	share_som_input_slow - 1) .GT. 10*EPSILON(un)) THEN
2155	WRITE(numout,*) 'ERROR: sum of shares should be 0 ', &
2156	ipts, ivm, share_som_input_surface + &
2157	share_som_input_active + share_som_input_slow
2158	CALL ipslerr (3,'stomate_litter','sum of shares',&
2159	'share_som_input_xxx','')
2160
2161	ENDIF
2162
2163	ENDIF ! err_act.GT.1
2164
2165	! Store old som_inputs in a temporary variable
2166	DO icarb = 1, ncarb
2167	som_input_old(icarb,iele) = som_input(ipts,icarb,ivm,iele)*dt
2168	ENDDO
2169
2170	! Truncate som_input values
2171	som_input(ipts,iactive,ivm,iele) = (share_som_input_active * &
2172	som_input_new(ivm,iele))
2173	som_input(ipts,islow,ivm,iele) = (share_som_input_slow * &
2174	som_input_new(ivm,iele))
2175	som_input(ipts,isurface,ivm,iele) = (share_som_input_surface * &
2176	som_input_new(ivm,iele))
2177
2178	ENDDO
2179
2180	! Convert to the initial units to avoid problems when closing the mass balance
2181	som_input(ipts,:,ivm,:) = som_input(ipts,:,ivm,:) / dt
2182
2183	! Distribute delta_som_input back over the litter pools
2184	! Calculate the share of each litter pool.
2185	DO iele = 1,nelements
2186
2187	! Calculate the share of each litter pool. In principle all values should be
2188	! larger than zero so the max is to prevent against zero stored as very small
2189	! negative numbers, e.g. -1e-17
2190	share_litter_struc_a(ipts,ivm,iele) = MAX(litter(ipts,istructural,ivm,iabove,iele) / &
2191	SUM(SUM(litter(ipts,:,ivm,:,iele),1)),zero)
2192	share_litter_struc_b(ipts,ivm,iele) = MAX(zero,litter(ipts,istructural,ivm,ibelow,iele) / &
2193	SUM(SUM(litter(ipts,:,ivm,:,iele),1)),zero)
2194	share_litter_woody_a(ipts,ivm,iele) = MAX(litter(ipts,iwoody,ivm,iabove,iele) / &
2195	SUM(SUM(litter(ipts,:,ivm,:,iele),1)),zero)
2196	share_litter_woody_b(ipts,ivm,iele) = MAX(litter(ipts,iwoody,ivm,ibelow,iele) / &
2197	SUM(SUM(litter(ipts,:,ivm,:,iele),1)),zero)
2198
2199	! All PFTs should have metabolic litter. Use the metabolic litter to ensure
2200	! mass balance closure. Not all PFTs have woody litter so the woody litter
2201	! cannot be used as the residual term. Note that the initial calculation of
2202	! share_litter_met_b as the residual of the other components could result
2203	! in a small negative number.
2204	share_litter_met_a(ipts,ivm,iele) = MAX(litter(ipts,imetabolic,ivm,iabove,iele) / &
2205	SUM(SUM(litter(ipts,:,ivm,:,iele),1)),zero)
2206	share_litter_met_b(ipts,ivm,iele) = un - (share_litter_struc_a(ipts,ivm,iele) + &
2207	share_litter_struc_b(ipts,ivm,iele) + share_litter_met_a(ipts,ivm,iele) + &
2208	share_litter_woody_a(ipts,ivm,iele) + share_litter_woody_b(ipts,ivm,iele))
2209
2210	! Correct precision errors caused by very small negative values for share_litter_met_b
2211	IF (share_litter_met_b(ipts,ivm,iele).LT.zero) THEN
2212
2213	IF (share_litter_met_b(ipts,ivm,iele).LT.-min_stomate) THEN
2214	CALL ipslerr(3,'stomate_litter','share_litter_met_b is negative',&
2215	'the value is too negative to be a precision error','')
2216	ELSE
2217
2218	! Value between 0 and -min_stomate. Correct it. The aboveground wood
2219	! pool is expected to be larger than the belowground pool so transfer
2220	! the precision issue to the aboverground pool.
2221	CALL ipslerr(2,'stomate_litter','share_litter_met_b is slightly negative',&
2222	'if this happens rarely it is OK','If it happens a lot debug!')
2223
2224	! Debug
2225	IF(printlev_loc.GT.4) WRITE(numout,*) 'initial share_litter_xxx, ', &
2226	share_litter_struc_a(ipts,ivm,iele), share_litter_struc_b(ipts,ivm,iele), &
2227	share_litter_woody_a(ipts,ivm,iele), share_litter_woody_b(ipts,ivm,iele), &
2228	share_litter_met_a(ipts,ivm,iele), share_litter_met_b(ipts,ivm,iele)
2229	!-
2230
2231	! If share_litter_met_a _ share_litter_met_b is positive it is easy to solve the
2232	! the precision problem
2233	IF(share_litter_met_a(ipts,ivm,iele)+share_litter_met_b(ipts,ivm,iele).GT.zero)THEN
2234
2235	! Move the precision error into share_litter_met_a
2236	share_litter_met_a(ipts,ivm,iele) = share_litter_met_a(ipts,ivm,iele) + &
2237	share_litter_met_b(ipts,ivm,iele)
2238	share_litter_met_b(ipts,ivm,iele) = zero
2239
2240	ELSE
2241
2242	! Difficult case. Nevertheless mass needs to be conserved. Set the negative value
2243	! to zero and then try to rescale all other share_litter_xxxs
2244	share_litter_met_b(ipts,ivm,iele) = zero
2245	total = MAX(zero, share_litter_struc_a(ipts,ivm,iele) + &
2246	share_litter_struc_b(ipts,ivm,iele) + share_litter_woody_a(ipts,ivm,iele) + &
2247	share_litter_woody_b(ipts,ivm,iele) + share_litter_met_a(ipts,ivm,iele))
2248
2249	IF(total.GT.zero)THEN
2250
2251	! Rescaling is possible
2252	share_litter_struc_a(ipts,ivm,iele)=share_litter_struc_a(ipts,ivm,iele)/total
2253	share_litter_struc_b(ipts,ivm,iele)=share_litter_struc_b(ipts,ivm,iele)/total
2254	share_litter_woody_a(ipts,ivm,iele)=share_litter_woody_a(ipts,ivm,iele)/total
2255	share_litter_woody_b(ipts,ivm,iele)=share_litter_woody_b(ipts,ivm,iele)/total
2256	share_litter_met_a(ipts,ivm,iele)=share_litter_met_a(ipts,ivm,iele)/total
2257
2258	! Debug
2259	IF(printlev_loc.GT.4) WRITE(numout,*) 'final share_litter_xxx, ', &
2260	share_litter_struc_a(ipts,ivm,iele), share_litter_struc_b(ipts,ivm,iele), &
2261	share_litter_woody_a(ipts,ivm,iele), share_litter_woody_b(ipts,ivm,iele), &
2262	share_litter_met_a(ipts,ivm,iele), share_litter_met_b(ipts,ivm,iele)
2263	!-
2264
2265	ELSE
2266
2267	WRITE(numout,*) 'ipts, ivm, iele, ', ipts,ivm, iele
2268	WRITE(numout,*) 'share_litter_xxx_a/b, ', share_litter_met_a(ipts,ivm,iele), &
2269	share_litter_met_b(ipts,ivm,iele), share_litter_struc_a(ipts,ivm,iele), &
2270	share_litter_struc_b(ipts,ivm,iele), share_litter_woody_a(ipts,ivm,iele), &
2271	share_litter_woody_b(ipts,ivm,iele)
2272	WRITE(numout,*) 'SUM(litter), ',SUM(SUM(litter(ipts,:,ivm,:,iele),1))
2273	CALL ipslerr(3,'stomate_litter','all share_litter_xxx_a/b turn out to be zero',&
2274	'This is very strange and should have resulted in a divide by zero earlier',&
2275	'')
2276
2277	END IF ! total.GT.zero
2278
2279	END IF ! share_litter_met_a+share_litter_met_b.EQ.zero
2280
2281	END IF ! share_litter_met_b.LT.-min_stomate
2282
2283	END IF ! share_litter_met_b.LT.zero
2284
2285	! Calculate the new litter value. When being decomposed, these values will not result
2286	! in negative soil_n_min values that cannot be accounted for through immobilisation.
2287	litter(ipts,istructural,ivm,iabove,iele) = share_litter_struc_a(ipts,ivm,iele) * &
2288	litter_total_new(iele)
2289	litter(ipts,istructural,ivm,ibelow,iele) = share_litter_struc_b(ipts,ivm,iele) * &
2290	litter_total_new(iele)
2291	litter(ipts,imetabolic,ivm,iabove,iele) = share_litter_met_a(ipts,ivm,iele) * &
2292	litter_total_new(iele)
2293	litter(ipts,imetabolic,ivm,ibelow,iele) = share_litter_met_b(ipts,ivm,iele) * &
2294	litter_total_new(iele)
2295	litter(ipts,iwoody,ivm,iabove,iele) = share_litter_woody_a(ipts,ivm,iele) * &
2296	litter_total_new(iele)
2297	litter(ipts,iwoody,ivm,ibelow,iele) = share_litter_woody_b(ipts,ivm,iele) * &
2298	litter_total_new(iele)
2299
2300	ENDDO ! nelements
2301
2302	ENDIF !ld_redistribute
2303
2304	ENDDO ! nvm
2305
2306	ENDDO ! npts
2307
2308	!! Reset variables
2309	resp_hetero_soil(:,:) = zero
2310	total_som_n_flux(:,:,:) = zero
2311
2312	!! 6. Check mass balance closure
2313
2314	!! 6.1 Calculate components of the mass balance
2315	! The carbon in turnover and bm_to_litter is moved to different pools but
2316	! ::turnover and ::bm_to_litter are never updated. They should either be
2317	! set to zero or not be included in the calculation of ::pool_end
2318	IF (err_act.GT.1) THEN
2319
2320	pool_end(:,:,:) = zero
2321	pool_end(:,:,initrogen) = pool_end(:,:,initrogen) + &
2322	n_mineralisation(:,:) * veget_max(:,:)
2323
2324	pool_end(:,:,initrogen) = pool_end(:,:,initrogen) + &
2325	(soil_n_min(:,:,iammonium) + soil_n_min(:,:,initrate))*veget_max(:,:)
2326
2327	DO ivm = 1,nvm
2328	pool_end(:,ivm,initrogen) = pool_end(:,ivm,initrogen) + &
2329	n_fungivores(:,ivm) * veget_max(:,ivm)
2330	ENDDO
2331
2332	DO iele = 1,nelements
2333
2334	! The litter pool
2335	DO ilitt = 1,nlitt
2336	DO ilev = 1,nlevs
2337	pool_end(:,:,iele) = pool_end(:,:,iele) + &
2338	(litter(:,ilitt,:,ilev,iele) * veget_max(:,:))
2339	ENDDO
2340	ENDDO
2341
2342	DO ivm = 1, nvm
2343	pool_end(:,ivm,iele) = pool_end(:,ivm,iele) + &
2344	SUM(harvest_pool_acc(:,ivm,:,iele),2)/area(:)
2345	END DO
2346	ENDDO
2347
2348	!! 6.2 Calculate mass balance
2349	check_intern(:,:,:,:) = zero
2350	check_intern(:,:,iland2atm,icarbon) = -un * resp_hetero_litter(:,:) * &
2351	veget_max(:,:) * dt
2352
2353	DO iele = 1,nelements
2354	check_intern(:,:,ilat2out,iele) = -un * &
2355	SUM(som_input(:,:,:,iele),2) * veget_max(:,:) * dt
2356	check_intern(:,:,ipoolchange,iele) = -un * (pool_end(:,:,iele) - &
2357	pool_start(:,:,iele))
2358	ENDDO
2359
2360	check_intern(:,:,iatm2land,initrogen) = &
2361	check_intern(:,:,iatm2land,initrogen) + &
2362	n_input(:,:,month,imanure) * dt * veget_max(:,:)
2363
2364	closure_intern = zero
2365	DO imbc = 1,nmbcomp
2366	DO iele=1,nelements
2367
2368	closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
2369	check_intern(:,:,imbc,iele)
2370	ENDDO
2371
2372	ENDDO
2373
2374	!! 6.3 Write the verdict
2375
2376	CALL check_mass_balance("littercalc", closure_intern, npts, pool_end, &
2377	pool_start, veget_max, 'pft')
2378
2379	ENDIF ! err_act.GT.1
2380
2381
2382	!! 7. calculate fraction of total soil covered by dead leaves
2383
2384	CALL deadleaf (npts, veget_max, dead_leaves, deadleaf_cover)
2385
2386
2387	!! 8. (Quasi-)Analytical Spin-up : Start filling MatrixA
2388
2389	IF (spinup_analytic) THEN
2390
2391	MatrixA(:,:,:,:) = zero
2392	VectorB(:,:,:) = zero
2393
2394
2395	DO ivm = 1,nvm
2396
2397	!- MatrixA : carbon fluxes leaving the litter
2398
2399	MatrixA(:,ivm ,istructural_above,istructural_above)= - dtlitter_dec_fac(istructural) &
2400	control_temp(:,iabove) * control_moist(:,iabove) * exp( -litter_struct_coef * lignin_struc(:,ivm ,iabove) )
2401
2402	MatrixA(:,ivm ,istructural_below,istructural_below) = - dtlitter_dec_fac(istructural) &
2403	control_temp(:,ibelow) * control_moist(:,ibelow) * exp( -litter_struct_coef * lignin_struc(:,ivm ,ibelow) )
2404
2405	MatrixA(:,ivm ,imetabolic_above,imetabolic_above) = - dtlitter_dec_fac(imetabolic) &
2406	control_temp(:,iabove) * control_moist(:,iabove)
2407
2408	MatrixA(:,ivm ,imetabolic_below,imetabolic_below) = - dtlitter_dec_fac(imetabolic) &
2409	control_temp(:,ibelow) * control_moist(:,ibelow)
2410
2411	! Flux leaving the woody above litter pool :
2412	MatrixA(:, ivm , iwoody_above, iwoody_above) = - dt * litter_dec_fac(iwoody) * control_temp(:,iabove) * &
2413	control_moist(:,iabove) * exp( -3. * lignin_wood(:,ivm ,iabove) )
2414
2415	! Flux leaving the woody below litter pool :
2416	MatrixA(:, ivm , iwoody_below, iwoody_below) = - dt * litter_dec_fac(iwoody) * control_temp(:,ibelow) * &
2417	control_moist(:,ibelow) * exp( -3. * lignin_wood(:,ivm ,ibelow))
2418
2419	! Flux received by the carbon surface from the woody above litter pool :
2420	MatrixA(:, ivm , isurface_pool, iwoody_above) = frac_soil(iwoody, isurface, iabove) * &
2421	dt litter_dec_fac(iwoody) &
2422	control_temp(:,iabove) * &
2423	control_moist(:,iabove) * &
2424	exp( -3. * lignin_wood(:,ivm ,iabove) ) * ( 1. - lignin_wood(:,ivm ,iabove) )
2425
2426	! Flux received by the carbon active from the woody below litter pool :
2427	MatrixA(:, ivm , iactive_pool, iwoody_below) = frac_soil(iwoody, iactive, ibelow) * &
2428	dt litter_dec_fac(iwoody) &
2429	control_temp(:,ibelow) * &
2430	control_moist(:,ibelow) * &
2431	exp( -3. * lignin_wood(:,ivm ,ibelow) ) * ( 1. - lignin_wood(:,ivm ,ibelow) )
2432
2433	! Flux received by the carbon slow from the woody above litter pool :
2434	MatrixA(:, ivm , islow_pool, iwoody_above) = frac_soil(iwoody, islow, iabove) * &
2435	dt litter_dec_fac(iwoody) &
2436	control_temp(:,iabove) * &
2437	control_moist(:,iabove) * &
2438	exp( -3. * lignin_wood(:,ivm ,iabove) ) * lignin_wood(:,ivm ,iabove)
2439
2440	! Flux received by the carbon slow from the woody below litter pool :
2441	MatrixA(:, ivm , islow_pool, iwoody_below) = frac_soil(iwoody, islow, ibelow) * &
2442	dt litter_dec_fac(iwoody) &
2443	control_temp(:,ibelow) * &
2444	control_moist(:,ibelow) * &
2445	exp( -3. * lignin_wood(:,ivm ,ibelow) ) * lignin_wood(:,ivm ,ibelow)
2446
2447
2448	!- MatrixA : carbon fluxes between the litter and the pools (the rest of the matrix is filled in stomate_soilcarbon.f90)
2449	MatrixA(:,ivm ,isurface_pool,istructural_above) = frac_soil(istructural,isurface,iabove) * &
2450	dtlitter_dec_fac(istructural) &
2451	control_temp(:,iabove) * control_moist(:,iabove) * &
2452	exp( -litter_struct_coef * lignin_struc(:,ivm ,iabove) ) * &
2453	( 1. - lignin_struc(:,ivm ,iabove) )
2454
2455
2456	MatrixA(:,ivm ,iactive_pool,istructural_below) = frac_soil(istructural,iactive,ibelow) * &
2457	dtlitter_dec_fac(istructural) &
2458	control_temp(:,ibelow) * control_moist(:,ibelow) * &
2459	exp( -litter_struct_coef * lignin_struc(:,ivm ,ibelow) ) * &
2460	( 1. - lignin_struc(:,ivm ,ibelow) )
2461
2462	MatrixA(:,ivm ,isurface_pool,imetabolic_above) = frac_soil(imetabolic,isurface,iabove) * &
2463	dtlitter_dec_fac(imetabolic) control_temp(:,iabove) * control_moist(:,iabove)
2464
2465	MatrixA(:,ivm ,iactive_pool,imetabolic_below) = frac_soil(imetabolic,iactive,ibelow) * &
2466	dtlitter_dec_fac(imetabolic) control_temp(:,ibelow) * control_moist(:,ibelow)
2467
2468	MatrixA(:,ivm ,islow_pool,istructural_above) = frac_soil(istructural,islow,iabove) * &
2469	dtlitter_dec_fac(istructural) &
2470	control_temp(:,iabove) * control_moist(:,iabove) * &
2471	exp( -litter_struct_coef * lignin_struc(:,ivm ,iabove) )* &
2472	lignin_struc(:,ivm ,iabove)
2473
2474
2475	MatrixA(:,ivm ,islow_pool,istructural_below) = frac_soil(istructural,islow,ibelow) * &
2476	dtlitter_dec_fac(istructural) &
2477	control_temp(:,ibelow) * control_moist(:,ibelow) * &
2478	exp( -litter_struct_coef * lignin_struc(:,ivm ,ibelow) )* &
2479	lignin_struc(:,ivm ,ibelow)
2480
2481
2482	!- VectorB : carbon input -
2483
2484	VectorB(:,ivm ,istructural_above) = litter_inc(:,istructural,ivm ,iabove,icarbon)
2485	VectorB(:,ivm ,istructural_below) = litter_inc(:,istructural,ivm ,ibelow,icarbon)
2486	VectorB(:,ivm ,imetabolic_above) = litter_inc(:,imetabolic,ivm ,iabove,icarbon)
2487	VectorB(:,ivm ,imetabolic_below) = litter_inc(:,imetabolic,ivm ,ibelow,icarbon)
2488	VectorB(:,ivm ,iwoody_above) = litter_inc(:,iwoody,ivm ,iabove,icarbon)
2489	VectorB(:,ivm ,iwoody_below) = litter_inc(:,iwoody,ivm ,ibelow,icarbon)
2490
2491	IF (printlev>=4) WRITE(numout,*) 'We filled MatrixA and VectorB'
2492
2493	WHERE(litter(:,istructural,ivm ,iabove,initrogen) .GT. min_stomate)
2494	CN_som_litter_longterm(:,ivm ,istructural_above) = &
2495	( CN_som_litter_longterm(:,ivm ,istructural_above) * (tau_CN_longterm-dt) &
2496	+ litter(:,istructural,ivm ,iabove,icarbon)/litter(:,istructural,ivm ,iabove,initrogen) * dt)/ (tau_CN_longterm)
2497	ENDWHERE
2498
2499	WHERE(litter(:,istructural,ivm ,ibelow,initrogen) .GT. min_stomate)
2500	CN_som_litter_longterm(:,ivm ,istructural_below) = &
2501	( CN_som_litter_longterm(:,ivm ,istructural_below) * (tau_CN_longterm-dt) &
2502	+ litter(:,istructural,ivm ,ibelow,icarbon)/litter(:,istructural,ivm ,ibelow,initrogen) * dt)/ (tau_CN_longterm)
2503	ENDWHERE
2504
2505	WHERE(litter(:,imetabolic,ivm ,iabove,initrogen) .GT. min_stomate)
2506	CN_som_litter_longterm(:,ivm ,imetabolic_above) = &
2507	( CN_som_litter_longterm(:,ivm ,imetabolic_above) * (tau_CN_longterm-dt) &
2508	+ litter(:,imetabolic,ivm ,iabove,icarbon)/litter(:,imetabolic,ivm ,iabove,initrogen) * dt)/ (tau_CN_longterm)
2509	ENDWHERE
2510
2511	WHERE(litter(:,imetabolic,ivm ,ibelow,initrogen) .GT. min_stomate)
2512	CN_som_litter_longterm(:,ivm ,imetabolic_below) = &
2513	( CN_som_litter_longterm(:,ivm ,imetabolic_below) * (tau_CN_longterm-dt) &
2514	+ litter(:,imetabolic,ivm ,ibelow,icarbon)/litter(:,imetabolic,ivm ,ibelow,initrogen) * dt)/ (tau_CN_longterm)
2515	ENDWHERE
2516
2517	WHERE(litter(:,iwoody,ivm ,iabove,initrogen) .GT. min_stomate)
2518	CN_som_litter_longterm(:,ivm ,iwoody_above) = &
2519	( CN_som_litter_longterm(:,ivm ,iwoody_above) * (tau_CN_longterm-dt) &
2520	+ litter(:,iwoody,ivm ,iabove,icarbon)/litter(:,iwoody,ivm ,iabove,initrogen) * dt)/ (tau_CN_longterm)
2521	ENDWHERE
2522
2523	WHERE(litter(:,iwoody,ivm ,ibelow,initrogen) .GT. min_stomate)
2524	CN_som_litter_longterm(:,ivm ,iwoody_below) = &
2525	( CN_som_litter_longterm(:,ivm ,iwoody_below) * (tau_CN_longterm-dt) &
2526	+ litter(:,iwoody,ivm ,ibelow,icarbon)/litter(:,iwoody,ivm ,ibelow,initrogen) * dt)/ (tau_CN_longterm)
2527	ENDWHERE
2528
2529	ENDDO ! Loop over # PFTs
2530
2531
2532	ENDIF ! spinup analytic
2533
2534	IF (printlev_loc>=4) WRITE(numout,*) 'Leaving littercalc'
2535
2536	END SUBROUTINE littercalc
2537
2538
2539	!! ==============================================================================================================================\n
2540	!! SUBROUTINE : deadleaf
2541	!!
2542	!>\BRIEF This routine calculates the deadleafcover.
2543	!!
2544	!! DESCRIPTION : It first calculates the lai corresponding to the dead leaves (LAI) using
2545	!! the dead leaves carbon content (DL) the specific leaf area (sla) and the
2546	!! maximal coverage fraction of a PFT (veget_max) using the following equations:
2547	!! \latexonly
2548	!! \input{deadleaf1.tex}
2549	!! \endlatexonly
2550	!! \n
2551	!! Then, the dead leaf cover (DLC) is calculated as following:\n
2552	!! \latexonly
2553	!! \input{deadleaf2.tex}
2554	!! \endlatexonly
2555	!! \n
2556	!!
2557	!! RECENT CHANGE(S) : None
2558	!!
2559	!! MAIN OUTPUT VARIABLE: ::deadleaf_cover
2560	!!
2561	!! REFERENCE(S) : None
2562	!!
2563	!! FLOWCHART : None
2564	!! \n
2565	!_ ================================================================================================================================
2566
2567	SUBROUTINE deadleaf (npts, veget_max, dead_leaves, deadleaf_cover)
2568
2569	!! 0. Variable and parameter declaration
2570
2571	!! 0.1 Input variables
2572
2573	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels (unitless)
2574	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
2575	!! metabolic and structural
2576	!! @tex $(gC m^{-2})$ @endtex
2577	REAL(r_std),DIMENSION(:,:),INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT defined in
2578	!! the input vegetation map
2579	!! @tex $(m^2 m^{-2})$ @endtex
2580
2581	!! 0.2 Output variables
2582
2583	REAL(r_std), DIMENSION(:), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves over all PFTs
2584	!! (0-1, unitless)
2585
2586	!! 0.3 Modified variables
2587
2588	!! 0.4 Local variables
2589
2590	REAL(r_std), DIMENSION(npts) :: dead_lai !! LAI of dead leaves @tex $(m^2 m^{-2})$ @endtex
2591	INTEGER(i_std) :: ivm !! Index (unitless)
2592	REAL(r_std), DIMENSION(npts) :: total_dead_leaves !! imetabolic + istructural
2593	!_ ================================================================================================================================
2594
2595	!! 1. LAI of dead leaves
2596
2597	IF (printlev>=3) WRITE(numout,*) 'Entering deadleaf'
2598
2599	! Debug
2600	IF (printlev_loc>=4) THEN
2601	WRITE(numout,*) 'dead_leaves, imetabolic ',dead_leaves(test_grid,test_pft,imetabolic)
2602	WRITE(numout,*) 'dead_leaves, istructural ',dead_leaves(test_grid,test_pft,istructural)
2603	END IF
2604	!-
2605
2606	dead_lai(:) = zero
2607
2608	DO ivm = 1,nvm !Loop over PFTs
2609	total_dead_leaves(:) = dead_leaves(:,ivm,imetabolic) + dead_leaves(:,ivm,istructural)
2610	dead_lai(:) = dead_lai(:) + biomass_to_lai( total_dead_leaves, npts, ivm) &
2611	* veget_max(:,ivm)
2612	ENDDO
2613
2614	!! 2. fraction of soil covered by dead leaves
2615
2616	IF (printlev_loc>=4) WRITE(numout,*) 'dead_lai, ', dead_lai(test_grid)
2617
2618
2619	deadleaf_cover(:) = un - exp( - 0.5 * dead_lai(:) )
2620
2621	IF (printlev>=4) WRITE(numout,*) 'Leaving deadleaf'
2622
2623	END SUBROUTINE deadleaf
2624
2625
2626	!! ================================================================================================================================
2627	!! FUNCTION : control_moist_func
2628	!!
2629	!>\BRIEF Calculate moisture control for litter and soil C decomposition
2630	!!
2631	!! DESCRIPTION : Calculate moisture control factor applied
2632	!! to litter decomposition and to soil carbon decomposition in
2633	!! stomate_soilcarbon.f90 using the following equation: \n
2634	!! \latexonly
2635	!! \input{control_moist_func1.tex}
2636	!! \endlatexonly
2637	!! \n
2638	!! with M the moisture control factor and soilmoisutre, the soil moisture
2639	!! calculated in sechiba.
2640	!! Then, the function is ranged between Moistcont_min and 1:\n
2641	!! \latexonly
2642	!! \input{control_moist_func2.tex}
2643	!! \endlatexonly
2644	!! \n
2645	!! RECENT CHANGE(S) : None
2646	!!
2647	!! RETURN VALUE : ::moistfunc_result
2648	!!
2649	!! REFERENCE(S) : None
2650	!!
2651	!! FLOWCHART : None
2652	!! \n
2653	!_ ================================================================================================================================
2654
2655	FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)
2656
2657	!! 0. Variable and parameter declaration
2658
2659	!! 0.1 Input variables
2660
2661	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixel (unitless)
2662	REAL(r_std), DIMENSION(:), INTENT(in) :: moist_in !! relative humidity (unitless)
2663
2664	!! 0.2 Output variables
2665
2666	REAL(r_std), DIMENSION(npts) :: moistfunc_result !! Moisture control factor (0.25-1, unitless)
2667
2668	!! 0.3 Modified variables
2669
2670	!! 0.4 Local variables
2671
2672	!_ ================================================================================================================================
2673
2674	moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
2675	moistfunc_result(:) = MAX( moistcont_min, MIN( un, moistfunc_result(:) ) )
2676
2677	END FUNCTION control_moist_func
2678
2679
2680	!! ================================================================================================================================
2681	!! FUNCTION : control_temp_func
2682	!!
2683	!>\BRIEF Calculate temperature control for litter and soild C decomposition
2684	!!
2685	!! DESCRIPTION : Calculate temperature control factor applied
2686	!! to litter decomposition and to soil carbon decomposition in
2687	!! stomate_soilcarbon.f90 using the following equation: \n
2688	!! \latexonly
2689	!! \input{control_temp_func1.tex}
2690	!! \endlatexonly
2691	!! \n
2692	!! with T the temperature control factor, temp the temperature in Kelvin of
2693	!! the air (for aboveground litter) or of the soil (for belowground litter
2694	!! and soil)
2695	!! Then, the function is limited in its maximal range to 1:\n
2696	!! \latexonly
2697	!! \input{control_temp_func2.tex}
2698	!! \endlatexonly
2699	!! \n
2700	!! RECENT CHANGE(S) : None
2701	!!
2702	!! RETURN VALUE: ::tempfunc_result
2703	!!
2704	!! REFERENCE(S) : None
2705	!!
2706	!! FLOWCHART : None
2707	!! \n
2708	!_ ================================================================================================================================
2709
2710	FUNCTION control_temp_func (npts, temp_in) RESULT (tempfunc_result)
2711
2712	!! 0. Variable and parameter declaration
2713
2714	!! 0.1 Input variables
2715	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of land pixels (unitless)
2716	REAL(r_std), DIMENSION(:), INTENT(in) :: temp_in !! Temperature (K)
2717
2718	!! 0.2 Output variables
2719	REAL(r_std), DIMENSION(npts) :: tempfunc_result !! Temperature control factor (0-1, unitless)
2720
2721	!! 0.3 Modified variables
2722
2723	!! 0.4 Local variables
2724
2725	!_ ================================================================================================================================
2726
2727	tempfunc_result(:) = exp( soil_Q10 * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
2728	tempfunc_result(:) = MIN( un, tempfunc_result(:) )
2729
2730	END FUNCTION control_temp_func
2731
2732	END MODULE stomate_litter

Note: See TracBrowser for help on using the repository browser.

Download in other formats: