Context Navigation

source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_stomate/stomate_litter.f90 @ 8367

Last change on this file since 8367 was 7245, checked in by nicolas.vuichard, 3 years ago
improve Carbon mass balance closure. See ticket #785
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 60.3 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: biomass_to_lai
35
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
47	CONTAINS
48
49	!! ================================================================================================================================
50	!! SUBROUTINE : littercalc_clear
51	!!
52	!!\BRIEF Set the flag ::firstcall_litter to .TRUE. and as such activate section
53	!! 1.1 of the subroutine littercalc (see below).
54	!!
55	!! DESCRIPTION : None
56	!!
57	!! RECENT CHANGE(S) : None
58	!!
59	!! MAIN OUTPUT VARIABLE(S) : None
60	!!
61	!! REFERENCE(S) : None
62	!!
63	!! FLOWCHART : None
64	!! \n
65	!_ ================================================================================================================================
66
67	SUBROUTINE littercalc_clear
68	firstcall_litter =.TRUE.
69	END SUBROUTINE littercalc_clear
70
71
72	!! ================================================================================================================================
73	!! SUBROUTINE : littercalc
74	!!
75	!!\BRIEF Calculation of the litter decomposition and therefore of the
76	!! heterotrophic respiration from litter following Parton et al. (1987).
77	!!
78	!! DESCRIPTION : The littercal routine splits the litter in 4 pools:
79	!! aboveground metaboblic, aboveground structural, belowground metabolic and
80	!! belowground structural. the fraction (F) of plant material going to metabolic
81	!! and structural is defined following Parton et al. (1987)
82	!! \latexonly
83	!! \input{littercalc1.tex}
84	!! \endlatexonly
85	!! \n
86	!! where L is the lignin content of the plant carbon pools considered and CN
87	!! its CN ratio. L and CN are fixed parameters for each plant carbon pools,
88	!! therefore it is the ratio between each plant carbon pool within a PFT, which
89	!! controlled the part of the total litter, that will be considered as
90	!! recalcitrant (i.e. structural litter) or labile (i.e. metabolic litter).\n
91	!!
92	!! The routine calculates the fraction of aboveground litter which is metabolic
93	!! or structural (the litterpart variable) which is then used in lpj_fire.f90.\n
94	!!
95	!! In the section 2, the routine calculate the new plant material entering the
96	!! litter pools by phenological death of plants organs (corresponding to the
97	!! variable turnover) and by fire, herbivory and others non phenological causes
98	!! (variable bm_to_litter). This calculation is first done for each PFT and then
99	!! the values calculated for each PFT are added up. Following the same approach
100	!! the lignin content of the total structural litter is calculated and will be
101	!! then used as a factor control of the decomposition of the structural litter
102	!! (lignin_struc) in the section 5.1.2. A test is performed to avoid that we add
103	!! more lignin than structural litter. Finally, the variable litterpart is
104	!! updated.\n
105	!!
106	!! In the section 3 and 4 the temperature and the moisture controlling the
107	!! decomposition are calculated for above and belowground. For aboveground
108	!! litter, air temperature and litter moisture are calculated in sechiba and used
109	!! directly. For belowground, soil temperature and moisture are also calculated
110	!! in sechiba but are modulated as a function of the soil depth. The modulation
111	!! is a multiplying factor exponentially distributed between 0 (in depth) and 1
112	!! in surface.\n
113	!!
114	!! Then, in the section 5, the routine calculates the structural litter decomposition
115	!! (C) following first order kinetics following Parton et al. (1987).
116	!! \latexonly
117	!! \input{littercalc2.tex}
118	!! \endlatexonly
119	!! \n
120	!! with k the decomposition rate of the structural litter.
121	!! k corresponds to
122	!! \latexonly
123	!! \input{littercalc3.tex}
124	!! \endlatexonly
125	!! \n
126	!! with littertau the turnover rate, T a function of the temperature and M a function of
127	!! the moisture described below.\n
128	!!
129	!! Then, the fraction of dead leaves (DL) composed by aboveground structural litter is
130	!! calculated as following
131	!! \latexonly
132	!! \input{littercalc4.tex}
133	!! \endlatexonly
134	!! \n
135	!! with k the decomposition rate of the structural litter previously
136	!! described.\n
137	!!
138	!! In the section 5.1, the fraction of decomposed structural litter
139	!! incorporated to the soil (Input) and its associated heterotrophic respiration are
140	!! calculated. For structural litter, the C decomposed could go in the active
141	!! soil carbon pool or in the slow carbon, as described in
142	!! stomate_soilcarbon.f90.\n
143	!! \latexonly
144	!! \input{littercalc5.tex}
145	!! \endlatexonly
146	!! \n
147	!! with f a parameter describing the fraction of structural litter incorporated
148	!! into the considered soil carbon pool, C the amount of litter decomposed and L
149	!! the amount of lignin in the litter. The litter decomposed which is not
150	!! incorporated into the soil is respired.\n
151	!!
152	!! In the section 5.2, the fraction of decomposed metabolic litter
153	!! incorporated to the soil and its associated heterotrophic respiration are
154	!! calculated with the same approaches presented for 5.1 but no control factor
155	!! depending on the lignin content are used.\n
156	!!
157	!! In the section 6 the dead leaf cover is calculated through a call to the
158	!! deadleaf subroutine presented below.\n
159	!!
160	!! In the section 7, if the flag SPINUP_ANALYTIC is set to true, we fill MatrixA
161	!! and VectorB following Lardy(2011).
162	!!
163	!! MAIN OUTPUT VARIABLES: ::deadleaf_cover, ::resp_hetero_litter
164	!! ::control_temp, ::control_moist
165	!!
166	!! REFERENCES:
167	!! - Parton, WJ, Schimel, DS, Cole, CV, and Ojima, DS. 1987. Analysis
168	!! of factors controlling soil organic matter levels in Great Plains
169	!! grasslands. Soil Science Society of America journal (USA)
170	!! (51):1173-1179.
171	!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
172	!! Environmental Modelling & Software (2011), doi:10.1016\|j.envsoft.2011.05.016
173	!!
174	!! FLOWCHART :
175	!! \latexonly
176	!! \includegraphics(scale=0.5){littercalcflow.jpg}
177	!! \endlatexonly
178	!! \n
179	!_ ================================================================================================================================
180
181	SUBROUTINE littercalc (npts, &
182	turnover, bm_to_litter, tree_bm_to_litter,&
183	veget_cov_max, tsurf, tsoil, soilhum, litterhum, soil_n_min, &
184	input, harvest_above, litter, dead_leaves, lignin_struc, &
185	lignin_wood, n_mineralisation, deadleaf_cover, resp_hetero_litter, &
186	som_input, control_temp, control_moist, &
187	MatrixA, VectorB, CN_target, CN_som_litter_longterm, tau_CN_longterm, &
188	tsoil_decomp)
189
190	!! 0. Variable and parameter declaration
191
192	!! 0.1 Input variables
193
194	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels
195	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: turnover !! Turnover rates of plant biomass
196	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
197	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: bm_to_litter !! Conversion of biomass to litter
198	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
199	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: tree_bm_to_litter!! Conversion of biomass to litter
200	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
201	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_cov_max !! PFT "Maximal" coverage fraction of a PFT
202	!! defined in the input vegetation map
203	!! @tex $(m^2 m^{-2})$ @endtex
204	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf !! Temperature (K) at the surface
205	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: tsoil !! Soil temperature (K)
206	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: soilhum !! Daily soil humidity of each soil layer
207	!! (unitless)
208	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhum !! Daily litter humidity (unitless)
209	REAL(r_std), DIMENSION(npts,nvm,ninput), INTENT(in) :: input !! nitrogen inputs into the soil (gN/m**2/day)
210
211	!! 0.2 Output variables
212
213	REAL(r_std), DIMENSION(npts), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves
214	!! over all PFTs (0-1, unitless)
215	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_hetero_litter !! Litter heterotrophic respiration. The unit
216	!! is given by m^2 of ground.
217	!! @tex $(gC dt_sechiba one\_day^{-1}) m^{-2})$ @endtex
218	REAL(r_std), DIMENSION(npts,ncarb,nvm,nelements), INTENT(out) :: som_input !! Quantity of Carbon (or Nitrogen) going into SOM pools
219	!! from litter decomposition. The unit is
220	!! given by m^2 of ground
221	!! @tex $(gC(orN) m^{-2} dt\_slow^{-1})$ @endtex
222	REAL(r_std), DIMENSION(npts,nlevs), INTENT(out) :: control_temp !! Temperature control of heterotrophic
223	!! respiration, above and below (0-1,
224	!! unitless)
225	REAL(r_std), DIMENSION(npts,nlevs), INTENT(out) :: control_moist !! Moisture control of heterotrophic
226	!! respiration (0.25-1, unitless)
227	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(out) :: MatrixA !! Matrix containing the fluxes between the
228	!! carbon pools per sechiba time step
229	!! @tex $(gC.m^2.day^{-1})$ @endtex
230	REAL(r_std), DIMENSION(npts,nvm,nbpools), INTENT(out) :: VectorB !! Vector containing the litter increase per
231	!! sechiba time step
232	!! @tex $(gC m^{-2})$ @endtex
233	REAL(r_std), DIMENSION(npts,nvm,ncarb), INTENT(out) :: CN_target !! C to N ratio of SOM flux from one pool to another (gN m-2 dt-1)
234	REAL(r_std), DIMENSION(npts), INTENT(out) :: tsoil_decomp !! Temperature used for decompostition in soil (K)
235
236
237	!! 0.3 Modified variables
238
239	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter !! Metabolic and structural litter,above and
240	!! below ground. The unit is given by m^2 of
241	!! ground @tex $(gC m^{-2})$ @endtex
242	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
243	!! metabolic and structural in
244	!! @tex $(gC m^{-2})$ @endtex
245	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_struc !! Ratio Lignin content in structural litter,
246	!! above and below ground, (0-1, unitless)
247	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_wood !! Ratio Lignin/Carbon in woody litter,
248	!! above and below ground, (0-1, unitless)
249	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: n_mineralisation !! Mineral N pool (gN m-2)
250
251	REAL(r_std), DIMENSION(npts,nvm,nbpools), INTENT(inout) :: CN_som_litter_longterm !! Longterm CN ratio of litter and som pools (gC/gN)
252	REAL(r_std), INTENT(inout) :: tau_CN_longterm !! Counter used for calculating the longterm CN_ratio of som and litter pools [seconds]
253	REAL(r_std), DIMENSION(npts,nelements), INTENT(inout) :: harvest_above !! Harvest above ground biomass for
254	!! agriculture @tex $(gC m^{-2})$ @endtex
255	REAL(r_std), DIMENSION(npts,nvm,nnspec),INTENT(inout) :: soil_n_min !! mineral nitrogen in the soil (gN/m**2)
256
257	!! 0.4 Local variables
258
259	REAL(r_std) :: dt !! Number of sechiba(fast processes) time-step per day
260	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: litterfrac !! The fraction of leaves, wood, etc. that
261	!! goes into metabolic and structural
262	!! litterpools (0-1, unitless)
263	!$OMP THREADPRIVATE(litterfrac)
264	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: z_soil !! Soil levels (m)
265	!$OMP THREADPRIVATE(z_soil)
266	REAL(r_std), DIMENSION(npts) :: rpc !! Integration constant for vertical root
267	!! profiles (unitless)
268	REAL(r_std), SAVE, DIMENSION(nlitt) :: litter_turn !! Turnover time in litter pools (days)
269	!$OMP THREADPRIVATE(litter_turn)
270	REAL(r_std), SAVE, DIMENSION(nlitt,ncarb,nlevs) :: frac_soil !! Fraction of litter that goes into soil
271	!! (litter -> carbon, above and below). The
272	!! remaining part goes to the atmosphere
273	!$OMP THREADPRIVATE(frac_soil)
274	REAL(r_std), DIMENSION(npts) :: soilhum_decomp !! Humidity used for decompostition in soil
275	!! (unitless)
276	REAL(r_std), DIMENSION(npts) :: fd !! Fraction of structural or metabolic litter
277	!! decomposed (unitless)
278	REAL(r_std), DIMENSION(npts,nelements) :: qd !! Quantity of structural or metabolic litter
279	!! decomposed @tex $(gC m^{-2})$ @endtex
280	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: old_struc !! Old structural litter, above and below
281	!! @tex $(gC m^{-2})$ @endtex
282	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: old_woody !! Old woody litter, above and below
283	!! @tex $(gC m^{-2})$ @endtex
284	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements) :: litter_inc !! Increase of metabolic and structural
285	!! litter, above and below ground. The unit
286	!! is given by m^2 of ground.
287	!! @tex $(gC m^{-2})$ @endtex
288	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: lignin_struc_inc !! Lignin increase in structural litter,
289	!! above and below ground. The unit is given
290	!! by m^2 of ground.
291	!! @tex $(gC m^{-2})$ @endtex
292	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: lignin_wood_inc !! Lignin increase in woody litter,
293	!! above and below ground. The unit is given
294	!! by m^2 of ground.
295	!! @tex $(gC m^{-2})$ @endtex
296	REAL(r_std), DIMENSION(npts) :: zdiff_min !! Intermediate field for looking for minimum
297	!! of what? this is not used in the code.
298	!! [??CHECK] could we delete it?
299	CHARACTER(LEN=10), DIMENSION(nlitt) :: litter_str !! Messages to write output information about
300	!! the litter
301	CHARACTER(LEN=22), DIMENSION(nparts) :: part_str !! Messages to write output information about
302	!! the plant
303	CHARACTER(LEN=7), DIMENSION(ncarb) :: carbon_str !! Messages to write output information about
304	!! the soil carbon
305	CHARACTER(LEN=5), DIMENSION(nlevs) :: level_str !! Messages to write output information about
306	!! the level (aboveground or belowground litter)
307	INTEGER(i_std) :: i,j,k,l,m !! Indices (unitless)
308	REAL(r_std), DIMENSION(npts,nvm) :: f_soil_n_min !! soil_n_min function response used for defing C to N target ratios (-)
309	REAL(r_std), DIMENSION(npts,nvm) :: f_pnc !! plant nitrogen concentration function response used for defining C to N target ratios (-)
310	INTEGER(i_std) :: itarget !! target som pool
311	REAL(r_std), DIMENSION(npts,nvm,nparts) :: CN !! CN ratio of the litter pools
312	REAL(r_std), DIMENSION(npts,nelements) :: summ
313	INTEGER(i_std) :: ier !! Error handling
314	REAL(r_std), DIMENSION(nlitt) :: tree_litterfrac
315	!_ ================================================================================================================================
316
317	IF (printlev>=3) WRITE(numout,*) 'Entering littercalc'
318
319	!! 1. Initialisations of the different fields during the first call of the routine
320	dt = dt_sechiba/one_day
321	IF ( firstcall_litter ) THEN
322
323
324	IF ( .NOT. ALLOCATED(litterfrac) ) THEN
325	ALLOCATE(litterfrac(npts,nvm,nparts,nlitt))
326	ENDIF
327
328
329	!! 1.1.4 residence times in litter pools (days)
330	litter_turn(imetabolic) = turn_metabolic / one_year ! .5 years
331	litter_turn(istructural) = turn_struct / one_year ! 3 years
332	litter_turn(iwoody) = turn_woody / one_year !!!!???? 30 years (2.45)
333
334	!! 1.1.5 decomposition flux fraction that goes into soil
335	! (litter -> carbon, above and below)
336	! 1-frac_soil goes into atmosphere
337	frac_soil(:,:,:) = zero
338
339	! structural litter: lignin fraction goes into slow pool + respiration,
340	! rest into active pool + respiration
341	frac_soil(istructural,isurface,iabove) = frac_soil_struct_sua
342	frac_soil(istructural,iactive,ibelow) = frac_soil_struct_ab
343	frac_soil(istructural,islow,iabove) = frac_soil_struct_sa
344	frac_soil(istructural,islow,ibelow) = frac_soil_struct_sb
345
346	! metabolic litter: all goes into active pool + respiration.
347	! Nothing into slow or passive pool.
348	frac_soil(imetabolic,isurface,iabove) = frac_soil_metab_sua
349	frac_soil(imetabolic,iactive,ibelow) = frac_soil_metab_ab
350
351
352	!! 1.2 soil levels
353	ALLOCATE(z_soil(0:nslm), stat=ier)
354	IF ( ier /= 0 ) CALL ipslerr_p(3,'littercalc','Pb in allocate of z_soil','','')
355
356	z_soil(0) = zero
357	z_soil(1:nslm) = zlt(1:nslm)
358
359
360	!! 1.3 messages
361	litter_str(imetabolic) = 'metabolic'
362	litter_str(istructural) = 'structural'
363	litter_str(iwoody) = 'woody'
364
365	carbon_str(iactive) = 'active'
366	carbon_str(isurface) = 'surface'
367	carbon_str(islow) = 'slow'
368	carbon_str(ipassive) = 'passive'
369
370	level_str(iabove) = 'above'
371	level_str(ibelow) = 'below'
372
373	part_str(ileaf) = 'leaves'
374	part_str(isapabove) = 'sap above ground'
375	part_str(isapbelow) = 'sap below ground'
376	part_str(iheartabove) = 'heartwood above ground'
377	part_str(iheartbelow) = 'heartwood below ground'
378	part_str(iroot) = 'roots'
379	part_str(ifruit) = 'fruits'
380	part_str(icarbres) = 'carbohydrate reserve'
381	part_str(ilabile) = 'labile reserve'
382
383
384	IF (printlev >=2) THEN
385	WRITE(numout,*) 'litter:'
386	WRITE(numout,*) ' > C/N ratios: '
387	DO k = 1, nparts
388	WRITE(numout,*) ' ', part_str(k), ': ',CN_fix(k)
389	ENDDO
390
391	WRITE(numout,*) ' > Lignine/C ratios: '
392	DO k = 1, nparts
393	WRITE(numout,*) ' ', part_str(k), ': ',LC(:,k)
394	ENDDO
395
396	WRITE(numout,*) ' > scaling depth for decomposition (m): ',z_decomp
397
398	WRITE(numout,*) ' > minimal carbon residence time in litter pools (d):'
399	DO m = 1, nlitt
400	WRITE(numout,*) ' ',litter_str(m),':',litter_turn(m)
401	ENDDO
402
403	WRITE(numout,*) ' > litter decomposition flux fraction that really goes '
404	WRITE(numout,*) ' into carbon pools (rest into the atmosphere):'
405	DO m = 1, nlitt
406	DO l = 1, nlevs
407	DO k = 1, ncarb
408	WRITE(numout,*) ' ',litter_str(m),' ',level_str(l),' -> ',&
409	carbon_str(k),':', frac_soil(m,k,l)
410	ENDDO
411	ENDDO
412	ENDDO
413	END IF ! printlev >=2
414	firstcall_litter = .FALSE.
415
416	ENDIF
417
418
419	!! 1.1.3 litter fractions:
420	!! what fraction of leaves, wood, etc. goes into metabolic and structural litterpools
421
422	! PFT 1 needs to be initialised for when everything is printed below
423	CN(:,:,:) = zero
424
425	DO k = 1, nparts
426	WHERE((bm_to_litter(:,:,k,initrogen)+turnover(:,:,k,initrogen)).GT.min_stomate)
427	CN(:,:,k) = (bm_to_litter(:,:,k,icarbon)+turnover(:,:,k,icarbon))/ &
428	(bm_to_litter(:,:,k,initrogen)+turnover(:,:,k,initrogen))
429	ELSEWHERE
430	CN(:,:,k) = CN_fix(k)
431	ENDWHERE
432	DO j = 1,nvm
433
434	IF ( ((k == isapabove) .OR. (k == isapbelow) .OR. (k == iheartabove) .OR. (k == iheartbelow)) &
435	.AND. is_tree(j) ) THEN
436	litterfrac(:,j,k,iwoody) = 1.
437	litterfrac(:,j,k,imetabolic) = zero
438	litterfrac(:,j,k,istructural) = zero
439	ELSE
440	IF( (k == icarbres) .OR. (k == ilabile)) THEN
441	litterfrac(:,j,k,imetabolic) = 1.0
442	ELSE
443	litterfrac(:,j,k,imetabolic) = MAX(metabolic_ref_frac - metabolic_LN_ratio * LC(j,k) * CN(:,j,k),zero)
444	ENDIF
445	litterfrac(:,j,k,istructural) = 1. - litterfrac(:,j,k,imetabolic)
446	litterfrac(:,j,k,iwoody) = 0.0
447
448	ENDIF
449
450	END DO
451
452	END DO
453
454	tree_litterfrac(iwoody)=1.
455	tree_litterfrac(imetabolic) = zero
456	tree_litterfrac(istructural) = zero
457
458
459	CN_target(:,:,:)=0.0
460
461	! C to N target ratios of differnt pools
462	f_soil_n_min = MIN((soil_n_min(:,:,iammonium) + soil_n_min(:,:,initrate)), 2.)
463	! CN ratio ranges between 3 and 15 for active pool
464	! Figure 4 of Parton et al. (1993) show lower CN values for active pool of 2 rather than 3
465	CN_target(:,:,iactive)= CN_target_iactive_ref + CN_target_iactive_Nmin * f_soil_n_min(:,:)
466	! CN ratio ranges between 12 and 20 for slow pool
467	CN_target(:,:,islow)= CN_target_islow_ref + CN_target_islow_Nmin * f_soil_n_min(:,:)
468	! CN ratio ranges between 7 and 10 for passive pool
469	! OCN uses a fixed value of 9. (don't know why)
470	! Figure 4 of Parton et al. (1993) show lower CN values for passive pool of 3 rather than 7
471	CN_target(:,:,ipassive)= CN_target_ipassive_ref + CN_target_ipassive_Nmin * f_soil_n_min(:,:)
472
473
474	! plant nitrogen content (%)
475	WHERE(litter(:,imetabolic,:,iabove,icarbon)+litter(:,istructural,:,iabove,icarbon).GT.0.)
476	f_pnc(:,:)=MIN((litter(:,imetabolic,:,iabove,initrogen)+litter(:,istructural,:,iabove,initrogen)) / &
477	(2.(litter(:,imetabolic,:,iabove,icarbon)+litter(:,istructural,:,iabove,icarbon))) 100., 2.)
478	ELSEWHERE
479	f_pnc(:,:)=0.
480	ENDWHERE
481
482	CN_target(:,:,isurface) = CN_target_isurface_ref + CN_target_isurface_pnc * f_pnc(:,:)
483
484
485
486	!! 1.4 set output to zero
487	deadleaf_cover(:) = zero
488	resp_hetero_litter(:,:) = zero
489	som_input(:,:,:,:) = zero
490
491
492	!! 2. Add biomass to different litterpools (per m^2 of ground)
493
494	!! 2.1 first, save old structural litter (needed for lignin fractions).
495	! above/below
496	DO l = 1, nlevs !Loop over litter levels (above and below ground)
497	DO m = 1,nvm !Loop over PFTs
498
499	old_struc(:,m,l) = litter(:,istructural,m,l,icarbon)
500	old_woody(:,m,l) = litter(:,iwoody,m,l,icarbon)
501	ENDDO
502	ENDDO
503
504
505	!! 2.2 update litter, dead leaves, and lignin content in structural litter
506	litter_inc(:,:,:,:,:) = zero
507	lignin_struc_inc(:,:,:) = zero
508	lignin_wood_inc(:,:,:) = zero
509
510	DO j = 1,nvm !Loop over PFTs
511
512	!! 2.2.1 litter
513	DO k = 1, nlitt !Loop over litter pools (metabolic and structural)
514
515	DO l = 1, nelements ! Loop over element pools (carbon and nitrogen)
516	!! 2.2.2 calculate litter increase (per m^2 of ground).
517	! Only a given fracion of fruit turnover is directly coverted into litter.
518	! Litter increase for each PFT, structural and metabolic, above/below
519	litter_inc(:,k,j,iabove,l) = &
520	litterfrac(:,j,ileaf,k) * bm_to_litter(:,j,ileaf,l) + &
521	litterfrac(:,j,isapabove,k) * (bm_to_litter(:,j,isapabove,l)-tree_bm_to_litter(:,j,isapabove,l)) + &
522	litterfrac(:,j,iheartabove,k) * (bm_to_litter(:,j,iheartabove,l)-tree_bm_to_litter(:,j,iheartabove,l)) + &
523	tree_litterfrac(k) * tree_bm_to_litter(:,j,isapabove,l) + &
524	tree_litterfrac(k) * tree_bm_to_litter(:,j,iheartabove,l) + &
525	litterfrac(:,j,ifruit,k) * bm_to_litter(:,j,ifruit,l) + &
526	litterfrac(:,j,icarbres,k) * bm_to_litter(:,j,icarbres,l) + &
527	litterfrac(:,j,ilabile,k) * bm_to_litter(:,j,ilabile,l) + &
528	litterfrac(:,j,ileaf,k) * turnover(:,j,ileaf,l) + &
529	litterfrac(:,j,isapabove,k) * turnover(:,j,isapabove,l) + &
530	litterfrac(:,j,iheartabove,k) * turnover(:,j,iheartabove,l) + &
531	litterfrac(:,j,ifruit,k) * turnover(:,j,ifruit,l) + &
532	litterfrac(:,j,icarbres,k) * turnover(:,j,icarbres,l)+ &
533	litterfrac(:,j,ilabile,k) * turnover(:,j,ilabile,l)
534
535	litter_inc(:,k,j,ibelow,l) = &
536	litterfrac(:,j,isapbelow,k) * (bm_to_litter(:,j,isapbelow,l)-tree_bm_to_litter(:,j,isapbelow,l)) + &
537	litterfrac(:,j,iheartbelow,k) * (bm_to_litter(:,j,iheartbelow,l)-tree_bm_to_litter(:,j,iheartbelow,l)) + &
538	tree_litterfrac(k) * tree_bm_to_litter(:,j,isapbelow,l) + &
539	tree_litterfrac(k) * tree_bm_to_litter(:,j,iheartbelow,l) + &
540	litterfrac(:,j,iroot,k) * bm_to_litter(:,j,iroot,l) + &
541	litterfrac(:,j,isapbelow,k) * turnover(:,j,isapbelow,l) + &
542	litterfrac(:,j,iheartbelow,k) * turnover(:,j,iheartbelow,l) + &
543	litterfrac(:,j,iroot,k) * turnover(:,j,iroot,l)
544
545	ENDDO
546	!! 2.2.3 dead leaves, for soil cover.
547	dead_leaves(:,j,k) = &
548	dead_leaves(:,j,k) + &
549	litterfrac(:,j,ileaf,k) * ( bm_to_litter(:,j,ileaf,icarbon) + turnover(:,j,ileaf,icarbon) )
550	ENDDO
551
552
553	!! 2.2.4 lignin increase in structural litter
554	lignin_struc_inc(:,j,iabove) = &
555	LC(j,ileaf) * bm_to_litter(:,j,ileaf,icarbon) + &
556	LC(j,ifruit) * bm_to_litter(:,j,ifruit,icarbon) + &
557	LC(j,icarbres) * bm_to_litter(:,j,icarbres,icarbon) + &
558	LC(j,ilabile) * bm_to_litter(:,j,ilabile,icarbon) + &
559	LC(j,ileaf) * turnover(:,j,ileaf,icarbon) + &
560	LC(j,ifruit) * turnover(:,j,ifruit,icarbon) + &
561	LC(j,icarbres) * turnover(:,j,icarbres,icarbon) + &
562	LC(j,ilabile) * turnover(:,j,ilabile,icarbon)
563
564	lignin_struc_inc(:,j,ibelow) = &
565	LC(j,iroot) * bm_to_litter(:,j,iroot,icarbon) + &
566	LC(j,iroot)*turnover(:,j,iroot,icarbon)
567
568
569	IF(is_tree(j)) THEN
570	!! 2.2.4 lignin increase in woody litter
571
572	lignin_wood_inc(:,j,iabove) = &
573	LC(j,isapabove) * turnover(:,j,isapabove,icarbon) + &
574	LC(j,iheartabove) * turnover(:,j,iheartabove,icarbon)
575
576	lignin_wood_inc(:,j,ibelow) = &
577	LC(j,isapbelow)*turnover(:,j,isapbelow,icarbon) + &
578	LC(j,iheartbelow)*turnover(:,j,iheartbelow,icarbon)
579	ELSE
580	!! 2.2.4 lignin increase in structural litter
581	lignin_struc_inc(:,j,iabove) = lignin_struc_inc(:,j,iabove) + &
582	LC(j,isapabove) * turnover(:,j,isapabove,icarbon) + &
583	LC(j,iheartabove) * turnover(:,j,iheartabove,icarbon)
584
585	lignin_struc_inc(:,j,ibelow) = lignin_struc_inc(:,j,ibelow) + &
586	LC(j,isapbelow)*turnover(:,j,isapbelow,icarbon) + &
587	LC(j,iheartbelow)*turnover(:,j,iheartbelow,icarbon)
588	ENDIF
589
590
591	lignin_wood_inc(:,j,iabove) = lignin_wood_inc(:,j,iabove)+&
592	LC(j,isapabove) * tree_bm_to_litter(:,j,isapabove,icarbon) + &
593	LC(j,iheartabove) * tree_bm_to_litter(:,j,iheartabove,icarbon)
594
595	lignin_wood_inc(:,j,ibelow) = lignin_wood_inc(:,j,ibelow)+ &
596	LC(j,isapbelow) * tree_bm_to_litter(:,j,isapbelow,icarbon) + &
597	LC(j,iheartbelow) * tree_bm_to_litter(:,j,iheartbelow,icarbon)
598
599	!! 2.2.4 lignin increase in structural litter
600	lignin_struc_inc(:,j,iabove) = lignin_struc_inc(:,j,iabove) + &
601	LC(j,isapabove) * (bm_to_litter(:,j,isapabove,icarbon)-tree_bm_to_litter(:,j,isapabove,icarbon)) + &
602	LC(j,iheartabove) * (bm_to_litter(:,j,iheartabove,icarbon)-tree_bm_to_litter(:,j,iheartabove,icarbon))
603
604	lignin_struc_inc(:,j,ibelow) = lignin_struc_inc(:,j,ibelow) + &
605	LC(j,isapbelow) * (bm_to_litter(:,j,isapbelow,icarbon)-tree_bm_to_litter(:,j,isapbelow,icarbon)) + &
606	LC(j,iheartbelow) * (bm_to_litter(:,j,iheartbelow,icarbon)-bm_to_litter(:,j,iheartbelow,icarbon))
607
608	ENDDO
609
610	!! 2.2.5 add new litter (struct/met, above/below)
611	litter(:,:,:,:,:) = litter(:,:,:,:,:) + litter_inc(:,:,:,:,:)
612
613	!! 2.2.6 for security: can't add more lignin than structural litter (above/below)
614	DO l = 1, nlevs !Loop over litter levels (above and below ground)
615	DO m = 1,nvm !Lopp over PFTs
616
617	lignin_struc_inc(:,m,l) = &
618	MIN( lignin_struc_inc(:,m,l), litter_inc(:,istructural,m,l,icarbon) )
619	lignin_wood_inc(:,m,l) = &
620	MIN( lignin_wood_inc(:,m,l), litter_inc(:,iwoody,m,l,icarbon) )
621
622	ENDDO
623	ENDDO
624
625
626	!! 2.2.7 new lignin content: add old lignin and lignin increase, divide by
627	!! total structural litter (above/below)
628
629	WHERE ( litter(:,istructural,:,:,icarbon) .GT. min_stomate )
630	lignin_struc(:,:,:) = &
631	( lignin_struc(:,:,:)*old_struc(:,:,:) + lignin_struc_inc(:,:,:) ) / &
632	litter(:,istructural,:,:,icarbon)
633	ELSEWHERE
634	lignin_struc(:,:,:) = zero
635	ENDWHERE
636
637	WHERE ( litter(:,iwoody,:,:,icarbon) .GT. min_stomate )
638	lignin_wood(:,:,:) = &
639	( lignin_wood(:,:,:)*old_woody(:,:,:) + lignin_wood_inc(:,:,:) ) / &
640	litter(:,iwoody,:,:,icarbon)
641	ELSEWHERE
642	lignin_wood(:,:,:) = zero
643	ENDWHERE
644
645	!! 2.2.8 Add the manure into the metabolic above ground pool if enough C is
646	!! harvested over the pixel. If not enough C is harvested a part goes as
647	!! litter and the other part goes at mineral N. Here we assume that N is not
648	!! limiting, only the C harvested is limiting
649
650	DO m=1,nvm
651
652	WHERE((veget_cov_max(:,m).GT.min_stomate) .AND. (harvest_above(:,icarbon)dt .GE. input(:,m,imanure)cn_ratio_manure*dt))
653
654	litter(:,imetabolic,m,iabove,icarbon) = litter(:,imetabolic,m,iabove,icarbon) + &
655	input(:,m,imanure)cn_ratio_manuredt
656	litter(:,imetabolic,m,iabove,initrogen) = litter(:,imetabolic,m,iabove,initrogen) + &
657	input(:,m,imanure)*dt
658
659	harvest_above(:,icarbon) = harvest_above(:,icarbon) - input(:,m,imanure)*cn_ratio_manure
660
661	ELSEWHERE((veget_cov_max(:,m).GT.min_stomate) .AND. (harvest_above(:,icarbon)dt .LT. input(:,m,imanure)cn_ratio_manure*dt))
662
663	litter(:,imetabolic,m,iabove,icarbon) = litter(:,imetabolic,m,iabove,icarbon) + &
664	harvest_above(:,icarbon)*dt
665	litter(:,imetabolic,m,iabove,initrogen) = litter(:,imetabolic,m,iabove,initrogen) + &
666	harvest_above(:,icarbon)/cn_ratio_manure*dt
667
668	soil_n_min(:,m,iammonium) = soil_n_min(:,m,iammonium) &
669	+ (input(:,m,imanure)-harvest_above(:,icarbon)/cn_ratio_manure)ratio_nh4_fertdt
670	soil_n_min(:,m,initrate) = soil_n_min(:,m,initrate) &
671	+ (input(:,m,imanure)-harvest_above(:,icarbon)/cn_ratio_manure)(1.-ratio_nh4_fert)dt
672
673	harvest_above(:,icarbon) = zero
674
675	ENDWHERE
676	ENDDO
677
678	!! 3. Temperature control on decay: Factor between 0 and 1
679
680	!! 3.1 above: surface temperature
681	control_temp(:,iabove) = control_temp_func (npts, tsurf)
682
683
684	!! 3.2 below: convolution of temperature and decomposer profiles
685	!! (exponential decomposer profile supposed)
686
687	!! 3.2.1 rpc is an integration constant such that the integral of the root profile is 1.
688	rpc(:) = un / ( un - EXP( -z_soil(nslm) / z_decomp ) )
689
690	!! 3.2.2 integrate over the nslm levels
691	tsoil_decomp(:) = zero
692
693	DO l = 1, nslm
694
695	tsoil_decomp(:) = &
696	tsoil_decomp(:) + tsoil(:,l) * rpc(:) * &
697	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
698
699	ENDDO
700
701	control_temp(:,ibelow) = control_temp_func (npts, tsoil_decomp)
702
703	!! 4. Moisture control. Factor between 0 and 1
704
705	!! 4.1 above the ground: litter humidity
706	control_moist(:,iabove) = control_moist_func (npts, litterhum)
707
708	!
709	!! 4.2 below: convolution of humidity and decomposer profiles
710	! (exponential decomposer profile supposed)
711
712	!! 4.2.1 rpc is an integration constant such that the integral of the root profile is 1.
713	rpc(:) = un / ( un - EXP( -z_soil(nslm) / z_decomp ) )
714
715	!! 4.2.2 integrate over the nslm levels
716	soilhum_decomp(:) = zero
717
718	DO l = 1, nslm !Loop over soil levels
719
720	soilhum_decomp(:) = &
721	soilhum_decomp(:) + soilhum(:,l) * rpc(:) * &
722	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
723
724	ENDDO
725
726	control_moist(:,ibelow) = control_moist_func (npts, soilhum_decomp)
727
728	!! 5. fluxes from litter to carbon pools and respiration
729
730	DO l = 1, nlevs !Loop over litter levels (above and below ground)
731	DO m = 1,nvm !Loop over PFTs
732
733	IF ( ok_soil_carbon_discretization ) THEN
734	itarget=iactive
735	ELSE
736	IF (l.EQ.iabove)THEN
737	itarget=isurface
738	ELSE
739	itarget=iactive
740	ENDIF
741	ENDIF
742	!! 5.1 structural litter: goes into active and slow carbon pools + respiration
743
744	!! 5.1.1 total quantity of structural litter which is decomposed
745	fd(:) = dtlitter_turn(istructural) &
746	control_temp(:,l) * control_moist(:,l) * exp( -litter_struct_coef * lignin_struc(:,m,l) )
747
748	DO k = 1,nelements
749
750	qd(:,k) = litter(:,istructural,m,l,k) * fd(:)
751
752	END DO
753
754	litter(:,istructural,m,l,:) = litter(:,istructural,m,l,:) - qd(:,:)
755	n_mineralisation(:,m) = n_mineralisation(:,m) + qd(:,initrogen)
756
757	!! 5.1.2 decompose same fraction of structural part of dead leaves. Not exact
758	!! as lignine content is not the same as that of the total structural litter.
759	! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
760	! we do this test to do this calcul only ones in 1,nlev loop
761	if (l == iabove) dead_leaves(:,m,istructural) = dead_leaves(:,m,istructural) * ( un - fd(:) )
762
763	!! 5.1.3 non-lignin fraction of structural litter goes into
764	!! active (or surface) carbon pool + respiration
765	som_input(:,itarget,m,icarbon) = som_input(:,itarget,m,icarbon) + &
766	frac_soil(istructural,itarget,l) * qd(:,icarbon) * ( 1. - lignin_struc(:,m,l) ) / dt
767
768	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
769	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
770	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
771	! time step and avoid some constructions that could create bug during future developments.
772	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
773	( 1. - frac_soil(istructural,itarget,l) ) * qd(:,icarbon) * &
774	( 1. - lignin_struc(:,m,l) ) / dt
775
776	!! 5.1.4 lignin fraction of structural litter goes into
777	!! slow carbon pool + respiration
778	som_input(:,islow,m,icarbon) = som_input(:,islow,m,icarbon) + &
779	frac_soil(istructural,islow,l) * qd(:,icarbon) * lignin_struc(:,m,l) / dt
780
781	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
782	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
783	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dt_sechiba/one_day
784	! time step and avoid some constructions that could create bug during future developments.
785	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
786	( 1. - frac_soil(istructural,islow,l) ) * qd(:,icarbon) * lignin_struc(:,m,l) / dt
787
788
789	!! 5.2 metabolic litter goes into active carbon pool + respiration
790
791	!! 5.2.1 total quantity of metabolic litter that is decomposed
792	fd(:) = dtlitter_turn(imetabolic) control_temp(:,l) * control_moist(:,l)
793
794	DO k = 1,nelements
795
796	qd(:,k) = litter(:,imetabolic,m,l,k) * fd(:)
797
798	END DO
799
800	litter(:,imetabolic,m,l,:) = litter(:,imetabolic,m,l,:) - qd(:,:)
801	n_mineralisation(:,m) = n_mineralisation(:,m) + qd(:,initrogen)
802	!! 5.2.2 decompose same fraction of metabolic part of dead leaves.
803	! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
804	! we do this test to do this calcul only ones in 1,nlev loop
805	if (l == iabove) dead_leaves(:,m,imetabolic) = dead_leaves(:,m,imetabolic) * ( 1. - fd(:) )
806
807	!! 5.2.3 put decomposed litter into active (or surface) pool + respiration
808	som_input(:,itarget,m,icarbon) = som_input(:,itarget,m,icarbon) + &
809	frac_soil(imetabolic,itarget,l) * qd(:,icarbon) / dt
810
811	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
812	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
813	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
814	! time step and avoid some constructions that could create bug during future developments.
815	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
816	( 1. - frac_soil(imetabolic,itarget,l) ) * qd(:,icarbon) / dt
817
818	!! 5.3 woody litter: goes into active and slow carbon pools + respiration
819
820	!! 5.3.1 total quantity of woody litter which is decomposed
821
822	fd(:) = dtlitter_turn(iwoody) &
823	control_temp(:,l) * control_moist(:,l) * EXP( -3. * lignin_wood(:,m,l) )
824
825	DO k = 1,nelements
826
827	qd(:,k) = litter(:,iwoody,m,l,k) * fd(:)
828
829	END DO
830
831	litter(:,iwoody,m,l,:) = litter(:,iwoody,m,l,:) - qd(:,:)
832	n_mineralisation(:,m) = n_mineralisation(:,m) + qd(:,initrogen)
833
834	!! 5.3.2 non-lignin fraction of woody litter goes into
835	!! active/structural carbon pool + respiration (per time unit)
836
837	som_input(:,itarget,m,icarbon) = som_input(:,itarget,m,icarbon) + &
838	frac_soil(istructural,itarget,l) * qd(:,icarbon) * ( 1. - lignin_wood(:,m,l) ) / dt
839
840	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
841	( 1. - frac_soil(istructural,itarget,l) ) * qd(:,icarbon) * &
842	( 1. - lignin_wood(:,m,l) ) / dt
843
844	!! 5.3.3 lignin fraction of woody litter goes into
845	!! slow carbon pool + respiration (per time unit)
846
847	som_input(:,islow,m,icarbon) = som_input(:,islow,m,icarbon) + &
848	frac_soil(istructural,islow,l) * qd(:,icarbon) * lignin_wood(:,m,l) / dt
849
850	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
851	( 1. - frac_soil(istructural,islow,l) ) * qd(:,icarbon) * lignin_wood(:,m,l) / dt
852
853
854	ENDDO
855	ENDDO
856
857	som_input(:,iactive,:,initrogen) = som_input(:,iactive,:,icarbon)/CN_target(:,:,iactive)
858	som_input(:,islow,:,initrogen) = som_input(:,islow,:,icarbon)/CN_target(:,:,islow)
859	som_input(:,isurface,:,initrogen) = som_input(:,isurface,:,icarbon)/CN_target(:,:,isurface)
860
861	n_mineralisation(:,:) = n_mineralisation(:,:) - &
862	( som_input(:,iactive,:,initrogen) + &
863	som_input(:,islow,:,initrogen) + &
864	som_input(:,isurface,:,initrogen))*dt !! multiply by dt !!
865
866
867	DO m=1,nvm
868	summ(:,:)=zero
869	DO l = 1, nlevs
870	DO k = 1,nlitt
871	summ(:,:)=summ(:,:)+litter(:,k,m,l,:)
872	ENDDO
873	ENDDO
874	WHERE(summ(:,icarbon).LE.min_stomate)
875	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
876	summ(:,icarbon)
877	n_mineralisation(:,m) = n_mineralisation(:,m) + &
878	summ(:,initrogen)
879	litter(:,imetabolic,m,iabove,icarbon) = zero
880	litter(:,istructural,m,iabove,icarbon) = zero
881	litter(:,iwoody,m,iabove,icarbon) = zero
882	litter(:,imetabolic,m,ibelow,icarbon) = zero
883	litter(:,istructural,m,ibelow,icarbon) = zero
884	litter(:,iwoody,m,ibelow,icarbon) = zero
885	litter(:,imetabolic,m,iabove,initrogen) = zero
886	litter(:,istructural,m,iabove,initrogen) = zero
887	litter(:,iwoody,m,iabove,initrogen) = zero
888	litter(:,imetabolic,m,ibelow,initrogen) = zero
889	litter(:,istructural,m,ibelow,initrogen) = zero
890	litter(:,iwoody,m,ibelow,initrogen) = zero
891	ENDWHERE
892	ENDDO
893
894	!! 6. calculate fraction of total soil covered by dead leaves
895
896	CALL deadleaf (npts, veget_cov_max, dead_leaves, deadleaf_cover)
897
898	!! 7. (Quasi-)Analytical Spin-up : Start filling MatrixA
899
900	IF (spinup_analytic) THEN
901
902	MatrixA(:,:,:,:) = zero
903	VectorB(:,:,:) = zero
904
905
906	DO m = 1,nvm
907
908	!- MatrixA : carbon fluxes leaving the litter
909
910	MatrixA(:,m,istructural_above,istructural_above)= - dtlitter_turn(istructural) &
911	control_temp(:,iabove) * control_moist(:,iabove) * exp( -litter_struct_coef * lignin_struc(:,m,iabove) )
912
913	MatrixA(:,m,istructural_below,istructural_below) = - dtlitter_turn(istructural) &
914	control_temp(:,ibelow) * control_moist(:,ibelow) * exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )
915
916	MatrixA(:,m,imetabolic_above,imetabolic_above) = - dtlitter_turn(imetabolic) &
917	control_temp(:,iabove) * control_moist(:,iabove)
918
919	MatrixA(:,m,imetabolic_below,imetabolic_below) = - dtlitter_turn(imetabolic) &
920	control_temp(:,ibelow) * control_moist(:,ibelow)
921
922	! Flux leaving the woody above litter pool :
923	MatrixA(:, m, iwoody_above, iwoody_above) = - dt * litter_turn(iwoody) * control_temp(:,iabove) * &
924	control_moist(:,iabove) * exp( -3. * lignin_wood(:,m,iabove) )
925
926	! Flux leaving the woody below litter pool :
927	MatrixA(:, m, iwoody_below, iwoody_below) = - dt * litter_turn(iwoody) * control_temp(:,ibelow) * &
928	control_moist(:,ibelow) * exp( -3. * lignin_wood(:,m,ibelow))
929
930	! Flux received by the carbon surface from the woody above litter pool :
931	MatrixA(:, m, isurface_pool, iwoody_above) = frac_soil(istructural, isurface, iabove) * &
932	dt litter_turn(iwoody) &
933	control_temp(:,iabove) * &
934	control_moist(:,iabove) * &
935	exp( -3. * lignin_wood(:,m,iabove) ) * ( 1. - lignin_wood(:,m,iabove) )
936
937	! Flux received by the carbon active from the woody below litter pool :
938	MatrixA(:, m, iactive_pool, iwoody_below) = frac_soil(istructural, iactive, ibelow) * &
939	dt litter_turn(iwoody) &
940	control_temp(:,ibelow) * &
941	control_moist(:,ibelow) * &
942	exp( -3. * lignin_wood(:,m,ibelow) ) * ( 1. - lignin_wood(:,m,ibelow) )
943
944	! Flux received by the carbon slow from the woody above litter pool :
945	MatrixA(:, m, islow_pool, iwoody_above) = frac_soil(istructural, islow, iabove) * &
946	dt litter_turn(iwoody) &
947	control_temp(:,iabove) * &
948	control_moist(:,iabove) * &
949	exp( -3. * lignin_wood(:,m,iabove) ) * lignin_wood(:,m,iabove)
950
951	! Flux received by the carbon slow from the woody below litter pool :
952	MatrixA(:, m, islow_pool, iwoody_below) = frac_soil(istructural, islow, ibelow) * &
953	dt litter_turn(iwoody) &
954	control_temp(:,ibelow) * &
955	control_moist(:,ibelow) * &
956	exp( -3. * lignin_wood(:,m,ibelow) ) * lignin_wood(:,m,ibelow)
957
958
959	!- MatrixA : carbon fluxes between the litter and the pools (the rest of the matrix is filled in stomate_soilcarbon.f90)
960	MatrixA(:,m,isurface_pool,istructural_above) = frac_soil(istructural,isurface,iabove) * &
961	dtlitter_turn(istructural) &
962	control_temp(:,iabove) * control_moist(:,iabove) * &
963	exp( -litter_struct_coef * lignin_struc(:,m,iabove) ) * &
964	( 1. - lignin_struc(:,m,iabove) )
965
966
967	MatrixA(:,m,iactive_pool,istructural_below) = frac_soil(istructural,iactive,ibelow) * &
968	dtlitter_turn(istructural) &
969	control_temp(:,ibelow) * control_moist(:,ibelow) * &
970	exp( -litter_struct_coef * lignin_struc(:,m,ibelow) ) * &
971	( 1. - lignin_struc(:,m,ibelow) )
972
973	MatrixA(:,m,isurface_pool,imetabolic_above) = frac_soil(imetabolic,isurface,iabove) * &
974	dtlitter_turn(imetabolic) control_temp(:,iabove) * control_moist(:,iabove)
975
976	MatrixA(:,m,iactive_pool,imetabolic_below) = frac_soil(imetabolic,iactive,ibelow) * &
977	dtlitter_turn(imetabolic) control_temp(:,ibelow) * control_moist(:,ibelow)
978
979	MatrixA(:,m,islow_pool,istructural_above) = frac_soil(istructural,islow,iabove) * &
980	dtlitter_turn(istructural) &
981	control_temp(:,iabove) * control_moist(:,iabove) * &
982	exp( -litter_struct_coef * lignin_struc(:,m,iabove) )* &
983	lignin_struc(:,m,iabove)
984
985
986	MatrixA(:,m,islow_pool,istructural_below) = frac_soil(istructural,islow,ibelow) * &
987	dtlitter_turn(istructural) &
988	control_temp(:,ibelow) * control_moist(:,ibelow) * &
989	exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )* &
990	lignin_struc(:,m,ibelow)
991
992
993	!- VectorB : carbon input -
994
995	VectorB(:,m,istructural_above) = litter_inc(:,istructural,m,iabove,icarbon)
996	VectorB(:,m,istructural_below) = litter_inc(:,istructural,m,ibelow,icarbon)
997	VectorB(:,m,imetabolic_above) = litter_inc(:,imetabolic,m,iabove,icarbon)
998	VectorB(:,m,imetabolic_below) = litter_inc(:,imetabolic,m,ibelow,icarbon)
999	VectorB(:,m,iwoody_above) = litter_inc(:,iwoody,m,iabove,icarbon)
1000	VectorB(:,m,iwoody_below) = litter_inc(:,iwoody,m,ibelow,icarbon)
1001
1002	IF (printlev>=4) WRITE(numout,*) 'We filled MatrixA and VectorB'
1003
1004	WHERE(litter(:,istructural,m,iabove,initrogen) .GT. min_stomate)
1005	CN_som_litter_longterm(:,m,istructural_above) = &
1006	( CN_som_litter_longterm(:,m,istructural_above) * (tau_CN_longterm-dt) &
1007	+ litter(:,istructural,m,iabove,icarbon)/litter(:,istructural,m,iabove,initrogen) * dt)/ (tau_CN_longterm)
1008	ENDWHERE
1009
1010	WHERE(litter(:,istructural,m,ibelow,initrogen) .GT. min_stomate)
1011	CN_som_litter_longterm(:,m,istructural_below) = &
1012	( CN_som_litter_longterm(:,m,istructural_below) * (tau_CN_longterm-dt) &
1013	+ litter(:,istructural,m,ibelow,icarbon)/litter(:,istructural,m,ibelow,initrogen) * dt)/ (tau_CN_longterm)
1014	ENDWHERE
1015
1016	WHERE(litter(:,imetabolic,m,iabove,initrogen) .GT. min_stomate)
1017	CN_som_litter_longterm(:,m,imetabolic_above) = &
1018	( CN_som_litter_longterm(:,m,imetabolic_above) * (tau_CN_longterm-dt) &
1019	+ litter(:,imetabolic,m,iabove,icarbon)/litter(:,imetabolic,m,iabove,initrogen) * dt)/ (tau_CN_longterm)
1020	ENDWHERE
1021
1022	WHERE(litter(:,imetabolic,m,ibelow,initrogen) .GT. min_stomate)
1023	CN_som_litter_longterm(:,m,imetabolic_below) = &
1024	( CN_som_litter_longterm(:,m,imetabolic_below) * (tau_CN_longterm-dt) &
1025	+ litter(:,imetabolic,m,ibelow,icarbon)/litter(:,imetabolic,m,ibelow,initrogen) * dt)/ (tau_CN_longterm)
1026	ENDWHERE
1027
1028	WHERE(litter(:,iwoody,m,iabove,initrogen) .GT. min_stomate)
1029	CN_som_litter_longterm(:,m,iwoody_above) = ( CN_som_litter_longterm(:,m,iwoody_above) * (tau_CN_longterm-dt) &
1030	+ litter(:,iwoody,m,iabove,icarbon)/litter(:,iwoody,m,iabove,initrogen) * dt)/ (tau_CN_longterm)
1031	ENDWHERE
1032
1033	WHERE(litter(:,iwoody,m,ibelow,initrogen) .GT. min_stomate)
1034	CN_som_litter_longterm(:,m,iwoody_below) = ( CN_som_litter_longterm(:,m,iwoody_below) * (tau_CN_longterm-dt) &
1035	+ litter(:,iwoody,m,ibelow,icarbon)/litter(:,iwoody,m,ibelow,initrogen) * dt)/ (tau_CN_longterm)
1036	ENDWHERE
1037
1038
1039
1040
1041	ENDDO ! Loop over # PFTs
1042
1043
1044	ENDIF ! spinup analytic
1045
1046	IF (printlev>=4) WRITE(numout,*) 'Leaving littercalc'
1047
1048	END SUBROUTINE littercalc
1049
1050
1051	!! ==============================================================================================================================\n
1052	!! SUBROUTINE : deadleaf
1053	!!
1054	!>\BRIEF This routine calculates the deadleafcover.
1055	!!
1056	!! DESCRIPTION : It first calculates the lai corresponding to the dead leaves (LAI) using
1057	!! the dead leaves carbon content (DL) the specific leaf area (sla) and the
1058	!! maximal coverage fraction of a PFT (vegetmax) using the following equations:
1059	!! \latexonly
1060	!! \input{deadleaf1.tex}
1061	!! \endlatexonly
1062	!! \n
1063	!! Then, the dead leaf cover (DLC) is calculated as following:\n
1064	!! \latexonly
1065	!! \input{deadleaf2.tex}
1066	!! \endlatexonly
1067	!! \n
1068	!!
1069	!! RECENT CHANGE(S) : None
1070	!!
1071	!! MAIN OUTPUT VARIABLE: ::deadleaf_cover
1072	!!
1073	!! REFERENCE(S) : None
1074	!!
1075	!! FLOWCHART : None
1076	!! \n
1077	!_ ================================================================================================================================
1078
1079	SUBROUTINE deadleaf (npts, veget_cov_max, dead_leaves, deadleaf_cover)
1080
1081	!! 0. Variable and parameter declaration
1082
1083	!! 0.1 Input variables
1084
1085	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels (unitless)
1086	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
1087	!! metabolic and structural
1088	!! @tex $(gC m^{-2})$ @endtex
1089	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_cov_max !! PFT "Maximal" coverage fraction of a PFT defined in
1090	!! the input vegetation map
1091	!! @tex $(m^2 m^{-2})$ @endtex
1092
1093	!! 0.2 Output variables
1094
1095	REAL(r_std), DIMENSION(npts), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves over all PFTs
1096	!! (0-1, unitless)
1097
1098	!! 0.3 Modified variables
1099
1100	!! 0.4 Local variables
1101
1102	REAL(r_std), DIMENSION(npts) :: dead_lai !! LAI of dead leaves @tex $(m^2 m^{-2})$ @endtex
1103	INTEGER(i_std) :: j !! Index (unitless)
1104	REAL(r_std), DIMENSION(npts) :: total_dead_leaves !! imetabolic + istructural
1105	!_ ================================================================================================================================
1106
1107	!! 1. LAI of dead leaves
1108
1109	dead_lai(:) = zero
1110
1111	DO j = 1,nvm !Loop over PFTs
1112	total_dead_leaves(:) = dead_leaves(:,j,imetabolic) + dead_leaves(:,j,istructural)
1113	dead_lai(:) = dead_lai(:) + biomass_to_lai( total_dead_leaves, npts, j) &
1114	* veget_cov_max(:,j)
1115	ENDDO
1116
1117	!! 2. fraction of soil covered by dead leaves
1118
1119	deadleaf_cover(:) = un - exp( - 0.5 * dead_lai(:) )
1120
1121	IF (printlev>=4) WRITE(numout,*) 'Leaving deadleaf'
1122
1123	END SUBROUTINE deadleaf
1124
1125
1126	!! ================================================================================================================================
1127	!! FUNCTION : control_moist_func
1128	!!
1129	!>\BRIEF Calculate moisture control for litter and soil C decomposition
1130	!!
1131	!! DESCRIPTION : Calculate moisture control factor applied
1132	!! to litter decomposition and to soil carbon decomposition in
1133	!! stomate_soilcarbon.f90 using the following equation: \n
1134	!! \latexonly
1135	!! \input{control_moist_func1.tex}
1136	!! \endlatexonly
1137	!! \n
1138	!! with M the moisture control factor and soilmoisutre, the soil moisture
1139	!! calculated in sechiba.
1140	!! Then, the function is ranged between Moistcont_min and 1:\n
1141	!! \latexonly
1142	!! \input{control_moist_func2.tex}
1143	!! \endlatexonly
1144	!! \n
1145	!! RECENT CHANGE(S) : None
1146	!!
1147	!! RETURN VALUE : ::moistfunc_result
1148	!!
1149	!! REFERENCE(S) : None
1150	!!
1151	!! FLOWCHART : None
1152	!! \n
1153	!_ ================================================================================================================================
1154
1155	FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)
1156
1157	!! 0. Variable and parameter declaration
1158
1159	!! 0.1 Input variables
1160
1161	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixel (unitless)
1162	REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in !! relative humidity (unitless)
1163
1164	!! 0.2 Output variables
1165
1166	REAL(r_std), DIMENSION(npts) :: moistfunc_result !! Moisture control factor (0.25-1, unitless)
1167
1168	!! 0.3 Modified variables
1169
1170	!! 0.4 Local variables
1171
1172	!_ ================================================================================================================================
1173
1174	moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
1175	moistfunc_result(:) = MAX( moistcont_min, MIN( un, moistfunc_result(:) ) )
1176
1177	END FUNCTION control_moist_func
1178
1179
1180	!! ================================================================================================================================
1181	!! FUNCTION : control_temp_func
1182	!!
1183	!>\BRIEF Calculate temperature control for litter and soild C decomposition
1184	!!
1185	!! DESCRIPTION : Calculate temperature control factor applied
1186	!! to litter decomposition and to soil carbon decomposition in
1187	!! stomate_soilcarbon.f90 using the following equation: \n
1188	!! \latexonly
1189	!! \input{control_temp_func1.tex}
1190	!! \endlatexonly
1191	!! \n
1192	!! with T the temperature control factor, temp the temperature in Kelvin of
1193	!! the air (for aboveground litter) or of the soil (for belowground litter
1194	!! and soil)
1195	!! Then, the function is limited in its maximal range to 1:\n
1196	!! \latexonly
1197	!! \input{control_temp_func2.tex}
1198	!! \endlatexonly
1199	!! \n
1200	!! RECENT CHANGE(S) : None
1201	!!
1202	!! RETURN VALUE: ::tempfunc_result
1203	!!
1204	!! REFERENCE(S) : None
1205	!!
1206	!! FLOWCHART : None
1207	!! \n
1208	!_ ================================================================================================================================
1209
1210	FUNCTION control_temp_func (npts, temp_in) RESULT (tempfunc_result)
1211
1212	!! 0. Variable and parameter declaration
1213
1214	!! 0.1 Input variables
1215	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of land pixels (unitless)
1216	REAL(r_std), DIMENSION(npts), INTENT(in) :: temp_in !! Temperature (K)
1217
1218	!! 0.2 Output variables
1219	REAL(r_std), DIMENSION(npts) :: tempfunc_result !! Temperature control factor (0-1, unitless)
1220
1221	!! 0.3 Modified variables
1222
1223	!! 0.4 Local variables
1224
1225	!_ ================================================================================================================================
1226
1227	tempfunc_result(:) = exp( soil_Q10 * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
1228	tempfunc_result(:) = MIN( un, tempfunc_result(:) )
1229
1230	END FUNCTION control_temp_func
1231
1232	END MODULE stomate_litter

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

Download in other formats: