Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-182/src_stomate/stomate_litter.f90

Last change on this file was 1170, checked in by josefine.ghattas, 11 years ago

Merge revision 1137 from branches/MERGE-OCN into the trunk.

Added supplementary dimension for carbon to following variables : bm_sapl, litter, bm_to_litter, turnover_daily, turnover_littercalc, turnover_longterm, bm_to_litttercalc and biomass. No change in results.

by Didier

Property svn:keywords set to HeadURL Date Author Revision

File size: 44.6 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_litter
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Update litter and lignine content after litter fall and
10	!! calculating litter decomposition.
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL$
20	!! $Date$
21	!! $Revision$
22	!! \n
23	!_ ================================================================================================================================
24
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
35	IMPLICIT NONE
36
37	! private & public routines
38
39	PRIVATE
40	PUBLIC littercalc,littercalc_clear, deadleaf
41
42	LOGICAL, SAVE :: firstcall = .TRUE. !! first call
43	!$OMP THREADPRIVATE(firstcall)
44
45	CONTAINS
46
47	!! ================================================================================================================================
48	!! SUBROUTINE : littercalc_calc
49	!!
50	!!\BRIEF Set the flag ::firstcall to .TRUE. and as such activate section
51	!! 1.1 of the subroutine littercalc (see below).
52	!!
53	!! DESCRIPTION : None
54	!!
55	!! RECENT CHANGE(S) : None
56	!!
57	!! MAIN OUTPUT VARIABLE(S) : None
58	!!
59	!! REFERENCE(S) : None
60	!!
61	!! FLOWCHART : None
62	!! \n
63	!_ ================================================================================================================================
64
65	SUBROUTINE littercalc_clear
66	firstcall =.TRUE.
67	END SUBROUTINE littercalc_clear
68
69
70	!! ================================================================================================================================
71	!! SUBROUTINE : littercalc
72	!!
73	!!\BRIEF Calculation of the litter decomposition and therefore of the
74	!! heterotrophic respiration from litter following Parton et al. (1987).
75	!!
76	!! DESCRIPTION : The littercal routine splits the litter in 4 pools:
77	!! aboveground metaboblic, aboveground structural, belowground metabolic and
78	!! belowground structural. the fraction (F) of plant material going to metabolic
79	!! and structural is defined following Parton et al. (1987)
80	!! \latexonly
81	!! \input{littercalc1.tex}
82	!! \endlatexonly
83	!! \n
84	!! where L is the lignin content of the plant carbon pools considered and CN
85	!! its CN ratio. L and CN are fixed parameters for each plant carbon pools,
86	!! therefore it is the ratio between each plant carbon pool within a PFT, which
87	!! controlled the part of the total litter, that will be considered as
88	!! recalcitrant (i.e. structural litter) or labile (i.e. metabolic litter).\n
89	!!
90	!! The routine calculates the fraction of aboveground litter which is metabolic
91	!! or structural (the litterpart variable) which is then used in lpj_fire.f90.\n
92	!!
93	!! In the section 2, the routine calculate the new plant material entering the
94	!! litter pools by phenological death of plants organs (corresponding to the
95	!! variable turnover) and by fire, herbivory and others non phenological causes
96	!! (variable bm_to_litter). This calculation is first done for each PFT and then
97	!! the values calculated for each PFT are added up. Following the same approach
98	!! the lignin content of the total structural litter is calculated and will be
99	!! then used as a factor control of the decomposition of the structural litter
100	!! (lignin_struc) in the section 5.1.2. A test is performed to avoid that we add
101	!! more lignin than structural litter. Finally, the variable litterpart is
102	!! updated.\n
103	!!
104	!! In the section 3 and 4 the temperature and the moisture controlling the
105	!! decomposition are calculated for above and belowground. For aboveground
106	!! litter, air temperature and litter moisture are calculated in sechiba and used
107	!! directly. For belowground, soil temperature and moisture are also calculated
108	!! in sechiba but are modulated as a function of the soil depth. The modulation
109	!! is a multiplying factor exponentially distributed between 0 (in depth) and 1
110	!! in surface.\n
111	!!
112	!! Then, in the section 5, the routine calculates the structural litter decomposition
113	!! (C) following first order kinetics following Parton et al. (1987).
114	!! \latexonly
115	!! \input{littercalc2.tex}
116	!! \endlatexonly
117	!! \n
118	!! with k the decomposition rate of the structural litter.
119	!! k corresponds to
120	!! \latexonly
121	!! \input{littercalc3.tex}
122	!! \endlatexonly
123	!! \n
124	!! with littertau the turnover rate, T a function of the temperature and M a function of
125	!! the moisture described below.\n
126	!!
127	!! Then, the fraction of dead leaves (DL) composed by aboveground structural litter is
128	!! calculated as following
129	!! \latexonly
130	!! \input{littercalc4.tex}
131	!! \endlatexonly
132	!! \n
133	!! with k the decomposition rate of the structural litter previously
134	!! described.\n
135	!!
136	!! In the section 5.1, the fraction of decomposed structural litter
137	!! incorporated to the soil (Input) and its associated heterotrophic respiration are
138	!! calculated. For structural litter, the C decomposed could go in the active
139	!! soil carbon pool or in the slow carbon, as described in
140	!! stomate_soilcarbon.f90.\n
141	!! \latexonly
142	!! \input{littercalc5.tex}
143	!! \endlatexonly
144	!! \n
145	!! with f a parameter describing the fraction of structural litter incorporated
146	!! into the considered soil carbon pool, C the amount of litter decomposed and L
147	!! the amount of lignin in the litter. The litter decomposed which is not
148	!! incorporated into the soil is respired.\n
149	!!
150	!! In the section 5.2, the fraction of decomposed metabolic litter
151	!! incorporated to the soil and its associated heterotrophic respiration are
152	!! calculated with the same approaches presented for 5.1 but no control factor
153	!! depending on the lignin content are used.\n
154	!!
155	!! In the section 6 the dead leaf cover is calculated through a call to the
156	!! deadleaf subroutine presented below.\n
157	!!
158	!! In the section 7, if the flag SPINUP_ANALYTIC is set to true, we fill MatrixA
159	!! and VectorB following Lardy(2011).
160	!!
161	!! MAIN OUTPUT VARIABLES: ::deadleaf_cover, ::resp_hetero_litter, ::soilcarbon_input,
162	!! ::control_temp, ::control_moist
163	!!
164	!! REFERENCES:
165	!! - Parton, WJ, Schimel, DS, Cole, CV, and Ojima, DS. 1987. Analysis
166	!! of factors controlling soil organic matter levels in Great Plains
167	!! grasslands. Soil Science Society of America journal (USA)
168	!! (51):1173-1179.
169	!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
170	!! Environmental Modelling & Software (2011), doi:10.1016\|j.envsoft.2011.05.016
171	!!
172	!! FLOWCHART :
173	!! \latexonly
174	!! \includegraphics(scale=0.5){littercalcflow.jpg}
175	!! \endlatexonly
176	!! \n
177	!_ ================================================================================================================================
178
179	SUBROUTINE littercalc (npts, dt, &
180	turnover, bm_to_litter, &
181	veget_max, tsurf, tsoil, soilhum, litterhum, &
182	litterpart, litter, dead_leaves, lignin_struc, &
183	deadleaf_cover, resp_hetero_litter, &
184	soilcarbon_input, control_temp, control_moist, &
185	MatrixA, VectorB)
186
187	!! 0. Variable and parameter declaration
188
189	!! 0.1 Input variables
190
191	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels
192	REAL(r_std), INTENT(in) :: dt !! Time step of the simulations for stomate
193	!! @tex $(dtradia one\_day^{-1})$ @endtex
194	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: turnover !! Turnover rates of plant biomass
195	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
196	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: bm_to_litter !! Conversion of biomass to litter
197	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
198	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
199	!! defined in the input vegetation map
200	!! @tex $(m^2 m^{-2})$ @endtex
201	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf !! Temperature (K) at the surface
202	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: tsoil !! Soil temperature (K)
203	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: soilhum !! Daily soil humidity of each soil layer
204	!! (unitless)
205	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhum !! Daily litter humidity (unitless)
206
207	!! 0.2 Output variables
208
209	REAL(r_std), DIMENSION(npts), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves
210	!! over all PFTs (0-1, unitless)
211	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_hetero_litter !! Litter heterotrophic respiration. The unit
212	!! is given by m^2 of ground.
213	!! @tex $(gC dtradia one\_day^{-1}) m^{-2})$ @endtex
214	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(out) :: soilcarbon_input !! Quantity of carbon going into carbon pools
215	!! from litter decomposition. The unit is
216	!! given by m^2 of ground
217	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
218	REAL(r_std), DIMENSION(npts,nlevs), INTENT(out) :: control_temp !! Temperature control of heterotrophic
219	!! respiration, above and below (0-1,
220	!! unitless)
221	REAL(r_std), DIMENSION(npts,nlevs), INTENT(out) :: control_moist !! Moisture control of heterotrophic
222	!! respiration (0.25-1, unitless)
223	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(out) :: MatrixA !! Matrix containing the fluxes between the
224	!! carbon pools per sechiba time step
225	!! @tex $(gC.m^2.day^{-1})$ @endtex
226	REAL(r_std), DIMENSION(npts,nvm,nbpools), INTENT(out) :: VectorB !! Vector containing the litter increase per
227	!! sechiba time step
228	!! @tex $(gC m^{-2})$ @endtex
229
230	!! 0.3 Modified variables
231
232	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: litterpart !! Fraction of litter above the ground
233	!! belonging to the different PFTs (0-1,
234	!! unitless)
235	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter !! Metabolic and structural litter,above and
236	!! below ground. The unit is given by m^2 of
237	!! ground @tex $(gC m^{-2})$ @endtex
238	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
239	!! metabolic and structural in
240	!! @tex $(gC m^{-2})$ @endtex
241	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_struc !! Ratio Lignin content in structural litter,
242	!! above and below ground, (0-1, unitless)
243
244	!! 0.4 Local variables
245
246	REAL(r_std), SAVE, DIMENSION(nparts,nlitt) :: litterfrac !! The fraction of leaves, wood, etc. that
247	!! goes into metabolic and structural
248	!! litterpools (0-1, unitless)
249	!$OMP THREADPRIVATE(litterfrac)
250	REAL(r_std), SAVE, DIMENSION(0:nbdl) :: z_soil !! Soil levels (m)
251	!$OMP THREADPRIVATE(z_soil)
252	REAL(r_std), DIMENSION(npts) :: rpc !! Integration constant for vertical root
253	!! profiles (unitless)
254	REAL(r_std), SAVE, DIMENSION(nlitt) :: litter_tau !! Turnover time in litter pools (days)
255	!$OMP THREADPRIVATE(litter_tau)
256	REAL(r_std), SAVE, DIMENSION(nlitt,ncarb,nlevs) :: frac_soil !! Fraction of litter that goes into soil
257	!! (litter -> carbon, above and below). The
258	!! remaining part goes to the atmosphere
259	!$OMP THREADPRIVATE(frac_soil)
260	REAL(r_std), DIMENSION(npts) :: tsoil_decomp !! Temperature used for decompostition in
261	!! soil (K)
262	REAL(r_std), DIMENSION(npts) :: soilhum_decomp !! Humidity used for decompostition in soil
263	!! (unitless)
264	REAL(r_std), DIMENSION(npts) :: fd !! Fraction of structural or metabolic litter
265	!! decomposed (unitless)
266	REAL(r_std), DIMENSION(npts,nelements) :: qd !! Quantity of structural or metabolic litter
267	!! decomposed @tex $(gC m^{-2})$ @endtex
268	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: old_struc !! Old structural litter, above and below
269	!! @tex $(gC m^{-2})$ @endtex
270	REAL(r_std), DIMENSION(npts,nvm,nlitt,nlevs,nelements) :: litter_inc_PFT !! Increase of litter, per PFT, metabolic and
271	!! structural, above and below ground. The
272	!! unit is given by m^2 of ground.
273	!! @tex $(gC m^{-2})$ @endtex
274	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements) :: litter_inc !! Increase of metabolic and structural
275	!! litter, above and below ground. The unit
276	!! is given by m^2 of ground.
277	!! @tex $(gC m^{-2})$ @endtex
278	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: lignin_struc_inc !! Lignin increase in structural litter,
279	!! above and below ground. The unit is given
280	!! by m^2 of ground.
281	!! @tex $(gC m^{-2})$ @endtex
282	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements) :: litter_pft !! Metabolic and structural litter above the
283	!! ground per PFT
284	REAL(r_std), DIMENSION(npts) :: zdiff_min !! Intermediate field for looking for minimum
285	!! of what? this is not used in the code.
286	!! [??CHECK] could we delete it?
287	CHARACTER(LEN=10), DIMENSION(nlitt) :: litter_str !! Messages to write output information about
288	!! the litter
289	CHARACTER(LEN=22), DIMENSION(nparts) :: part_str !! Messages to write output information about
290	!! the plant
291	CHARACTER(LEN=7), DIMENSION(ncarb) :: carbon_str !! Messages to write output information about
292	!! the soil carbon
293	CHARACTER(LEN=5), DIMENSION(nlevs) :: level_str !! Messages to write output information about
294	!! the level (aboveground or belowground litter)
295	INTEGER(i_std) :: i,j,k,l,m !! Indices (unitless)
296	!_ ================================================================================================================================
297
298	IF (bavard.GE.3) WRITE(numout,*) 'Entering littercalc'
299
300	!! 1. Initialisations of the different fields during the first call of the routine
301
302	IF ( firstcall ) THEN
303
304	!! 1.1.3 litter fractions:
305	!! what fraction of leaves, wood, etc. goes into metabolic and structural litterpools
306	DO k = 1, nparts
307
308	litterfrac(k,imetabolic) = metabolic_ref_frac - metabolic_LN_ratio * LC(k) * CN(k)
309	litterfrac(k,istructural) = un - litterfrac(k,imetabolic)
310
311	ENDDO
312
313	!! 1.1.4 residence times in litter pools (days)
314	litter_tau(imetabolic) = tau_metabolic * one_year ! .5 years
315	litter_tau(istructural) = tau_struct * one_year ! 3 years
316
317	!! 1.1.5 decomposition flux fraction that goes into soil
318	! (litter -> carbon, above and below)
319	! 1-frac_soil goes into atmosphere
320	frac_soil(:,:,:) = zero
321
322	! structural litter: lignin fraction goes into slow pool + respiration,
323	! rest into active pool + respiration
324	frac_soil(istructural,iactive,iabove) = frac_soil_struct_aa
325	frac_soil(istructural,iactive,ibelow) = frac_soil_struct_ab
326	frac_soil(istructural,islow,iabove) = frac_soil_struct_sa
327	frac_soil(istructural,islow,ibelow) = frac_soil_struct_sb
328
329	! metabolic litter: all goes into active pool + respiration.
330	! Nothing into slow or passive pool.
331	frac_soil(imetabolic,iactive,iabove) = frac_soil_metab_aa
332	frac_soil(imetabolic,iactive,ibelow) = frac_soil_metab_ab
333
334
335	!! 1.2 soil levels
336	z_soil(0) = zero
337	z_soil(1:nbdl) = diaglev(1:nbdl)
338
339
340	!! 1.3 messages
341	litter_str(imetabolic) = 'metabolic'
342	litter_str(istructural) = 'structural'
343
344	carbon_str(iactive) = 'active'
345	carbon_str(islow) = 'slow'
346	carbon_str(ipassive) = 'passive'
347
348	level_str(iabove) = 'above'
349	level_str(ibelow) = 'below'
350
351	part_str(ileaf) = 'leaves'
352	part_str(isapabove) = 'sap above ground'
353	part_str(isapbelow) = 'sap below ground'
354	part_str(iheartabove) = 'heartwood above ground'
355	part_str(iheartbelow) = 'heartwood below ground'
356	part_str(iroot) = 'roots'
357	part_str(ifruit) = 'fruits'
358	part_str(icarbres) = 'carbohydrate reserve'
359
360	WRITE(numout,*) 'litter:'
361
362	WRITE(numout,*) ' > C/N ratios: '
363	DO k = 1, nparts
364	WRITE(numout,*) ' ', part_str(k), ': ',CN(k)
365	ENDDO
366
367	WRITE(numout,*) ' > Lignine/C ratios: '
368	DO k = 1, nparts
369	WRITE(numout,*) ' ', part_str(k), ': ',LC(k)
370	ENDDO
371
372	WRITE(numout,*) ' > fraction of compartment that goes into litter: '
373	DO k = 1, nparts
374	DO m = 1, nlitt
375	WRITE(numout,*) ' ', part_str(k), '-> ',litter_str(m), ':',litterfrac(k,m)
376	ENDDO
377	ENDDO
378
379	WRITE(numout,*) ' > scaling depth for decomposition (m): ',z_decomp
380
381	WRITE(numout,*) ' > minimal carbon residence time in litter pools (d):'
382	DO m = 1, nlitt
383	WRITE(numout,*) ' ',litter_str(m),':',litter_tau(m)
384	ENDDO
385
386	WRITE(numout,*) ' > litter decomposition flux fraction that really goes '
387	WRITE(numout,*) ' into carbon pools (rest into the atmosphere):'
388	DO m = 1, nlitt
389	DO l = 1, nlevs
390	DO k = 1, ncarb
391	WRITE(numout,*) ' ',litter_str(m),' ',level_str(l),' -> ',&
392	carbon_str(k),':', frac_soil(m,k,l)
393	ENDDO
394	ENDDO
395	ENDDO
396
397	firstcall = .FALSE.
398
399	ENDIF
400
401
402	!! 1.3 litter above the ground per PFT.
403	DO j = 2, nvm ! Loop over # PFTs
404
405	DO k = 1, nlitt !Loop over litter pool
406	litter_pft(:,j,k,icarbon) = litterpart(:,j,k) * litter(:,k,j,iabove,icarbon)
407	ENDDO
408
409	ENDDO
410
411
412	!! 1.4 set output to zero
413	deadleaf_cover(:) = zero
414	resp_hetero_litter(:,:) = zero
415	soilcarbon_input(:,:,:) = zero
416
417
418	!! 2. Add biomass to different litterpools (per m^2 of ground)
419
420	!! 2.1 first, save old structural litter (needed for lignin fractions).
421	! above/below
422	DO l = 1, nlevs !Loop over litter levels (above and below ground)
423	DO m = 2,nvm !Loop over PFTs
424
425	old_struc(:,m,l) = litter(:,istructural,m,l,icarbon)
426
427	ENDDO
428	ENDDO
429
430
431	!! 2.2 update litter, dead leaves, and lignin content in structural litter
432	litter_inc(:,:,:,:,:) = zero
433	lignin_struc_inc(:,:,:) = zero
434
435	DO j = 2,nvm !Loop over PFTs
436
437	!! 2.2.1 litter
438	DO k = 1, nlitt !Loop over litter pools (metabolic and structural)
439
440	!! 2.2.2 calculate litter increase (per m^2 of ground).
441	! Only a given fracion of fruit turnover is directly coverted into litter.
442	! Litter increase for each PFT, structural and metabolic, above/below
443	litter_inc_PFT(:,j,k,iabove,icarbon) = &
444	litterfrac(ileaf,k) * bm_to_litter(:,j,ileaf,icarbon) + &
445	litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
446	litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
447	litterfrac(ifruit,k) * bm_to_litter(:,j,ifruit,icarbon) + &
448	litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
449	litterfrac(ileaf,k) * turnover(:,j,ileaf,icarbon) + &
450	litterfrac(isapabove,k) * turnover(:,j,isapabove,icarbon) + &
451	litterfrac(iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
452	litterfrac(ifruit,k) * turnover(:,j,ifruit,icarbon) + &
453	litterfrac(icarbres,k) * turnover(:,j,icarbres,icarbon)
454
455	litter_inc_PFT(:,j,k,ibelow,icarbon) = &
456	litterfrac(isapbelow,k) * bm_to_litter(:,j,isapbelow,icarbon) + &
457	litterfrac(iheartbelow,k) * bm_to_litter(:,j,iheartbelow,icarbon) + &
458	litterfrac(iroot,k) * bm_to_litter(:,j,iroot,icarbon) + &
459	litterfrac(isapbelow,k) * turnover(:,j,isapbelow,icarbon) + &
460	litterfrac(iheartbelow,k) * turnover(:,j,iheartbelow,icarbon) + &
461	litterfrac(iroot,k) * turnover(:,j,iroot,icarbon)
462
463	! litter increase, met/struct, above/below
464	litter_inc(:,k,j,iabove,icarbon) = litter_inc(:,k,j,iabove,icarbon) + litter_inc_PFT(:,j,k,iabove,icarbon)
465	litter_inc(:,k,j,ibelow,icarbon) = litter_inc(:,k,j,ibelow,icarbon) + litter_inc_PFT(:,j,k,ibelow,icarbon)
466
467	!! 2.2.3 dead leaves, for soil cover.
468	dead_leaves(:,j,k) = &
469	dead_leaves(:,j,k) + &
470	litterfrac(ileaf,k) * ( bm_to_litter(:,j,ileaf,icarbon) + turnover(:,j,ileaf,icarbon) )
471
472	!! 2.2.4 lignin increase in structural litter
473	IF ( k .EQ. istructural ) THEN
474
475	lignin_struc_inc(:,j,iabove) = &
476	lignin_struc_inc(:,j,iabove) + &
477	LC(ileaf) * bm_to_litter(:,j,ileaf,icarbon) + &
478	LC(isapabove) * bm_to_litter(:,j,isapabove,icarbon) + &
479	LC(iheartabove) * bm_to_litter(:,j,iheartabove,icarbon) + &
480	LC(ifruit) * bm_to_litter(:,j,ifruit,icarbon) + &
481	LC(icarbres) * bm_to_litter(:,j,icarbres,icarbon) + &
482	LC(ileaf) * turnover(:,j,ileaf,icarbon) + &
483	LC(isapabove) * turnover(:,j,isapabove,icarbon) + &
484	LC(iheartabove) * turnover(:,j,iheartabove,icarbon) + &
485	LC(ifruit) * turnover(:,j,ifruit,icarbon) + &
486	LC(icarbres) * turnover(:,j,icarbres,icarbon)
487
488	lignin_struc_inc(:,j,ibelow) = &
489	lignin_struc_inc(:,j,ibelow) + &
490	LC(isapbelow) * bm_to_litter(:,j,isapbelow,icarbon) + &
491	LC(iheartbelow) * bm_to_litter(:,j,iheartbelow,icarbon) + &
492	LC(iroot) * bm_to_litter(:,j,iroot,icarbon) + &
493	LC(isapbelow)*turnover(:,j,isapbelow,icarbon) + &
494	LC(iheartbelow)*turnover(:,j,iheartbelow,icarbon) + &
495	LC(iroot)*turnover(:,j,iroot,icarbon)
496
497	ENDIF
498
499	ENDDO
500	ENDDO
501
502	!! 2.2.5 add new litter (struct/met, above/below)
503	litter(:,:,:,:,:) = litter(:,:,:,:,:) + litter_inc(:,:,:,:,:)
504
505	!! 2.2.6 for security: can't add more lignin than structural litter (above/below)
506	DO l = 1, nlevs !Loop over litter levels (above and below ground)
507	DO m = 2,nvm !Lopp over PFTs
508
509	lignin_struc_inc(:,m,l) = &
510	MIN( lignin_struc_inc(:,m,l), litter_inc(:,istructural,m,l,icarbon) )
511
512	ENDDO
513	ENDDO
514
515	!! 2.2.7 new lignin content: add old lignin and lignin increase, divide by
516	!! total structural litter (above/below)
517	DO l = 1, nlevs !Loop over litter levels (above and below ground)
518	DO m = 2,nvm !Loop over PFTs
519	WHERE( litter(:,istructural,m,l,icarbon) .GT. min_stomate )
520
521	!MM : Soenke modif
522	! Best vectorization ?
523	!!$ lignin_struc(:,:,:) = &
524	!!$ ( lignin_struc(:,:,:)*old_struc(:,:,:) + lignin_struc_inc(:,:,:) ) / &
525	!!$ litter(:,istructural,:,:,icarbon)
526
527	lignin_struc(:,m,l) = lignin_struc(:,m,l) * old_struc(:,m,l)
528	lignin_struc(:,m,l) = lignin_struc(:,m,l) + lignin_struc_inc(:,m,l)
529	lignin_struc(:,m,l) = lignin_struc(:,m,l) / litter(:,istructural,m,l,icarbon)
530	ELSEWHERE
531	lignin_struc(:,m,l) = zero
532	ENDWHERE
533	ENDDO
534	ENDDO
535
536
537	!! 2.3 new litter fraction per PFT (for structural and metabolic litter, above
538	!! the ground).
539	DO j = 2,nvm !Loop over PFTs
540
541	WHERE ( litter(:,:,j,iabove,icarbon) .GT. min_stomate )
542
543	litterpart(:,j,:) = &
544	( litter_pft(:,j,:,icarbon) + litter_inc_PFT(:,j,:,iabove,icarbon) ) / litter(:,:,j,iabove,icarbon)
545
546	ELSEWHERE
547
548	litterpart(:,j,:) = zero
549
550	ENDWHERE
551
552	ENDDO
553
554
555	!! 3. Temperature control on decay: Factor between 0 and 1
556
557	!! 3.1 above: surface temperature
558	control_temp(:,iabove) = control_temp_func (npts, tsurf)
559
560
561	!! 3.2 below: convolution of temperature and decomposer profiles
562	!! (exponential decomposer profile supposed)
563
564	!! 3.2.1 rpc is an integration constant such that the integral of the root profile is 1.
565	rpc(:) = un / ( un - EXP( -z_soil(nbdl) / z_decomp ) )
566
567	!! 3.2.2 integrate over the nbdl levels
568	tsoil_decomp(:) = zero
569
570	DO l = 1, nbdl
571
572	tsoil_decomp(:) = &
573	tsoil_decomp(:) + tsoil(:,l) * rpc(:) * &
574	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
575
576	ENDDO
577
578	control_temp(:,ibelow) = control_temp_func (npts, tsoil_decomp)
579
580	!! 4. Moisture control. Factor between 0 and 1
581
582	!! 4.1 above the ground: litter humidity
583	control_moist(:,iabove) = control_moist_func (npts, litterhum)
584
585	!
586	!! 4.2 below: convolution of humidity and decomposer profiles
587	! (exponential decomposer profile supposed)
588
589	!! 4.2.1 rpc is an integration constant such that the integral of the root profile is 1.
590	rpc(:) = un / ( un - EXP( -z_soil(nbdl) / z_decomp ) )
591
592	!! 4.2.2 integrate over the nbdl levels
593	soilhum_decomp(:) = zero
594
595	DO l = 1, nbdl !Loop over soil levels
596
597	soilhum_decomp(:) = &
598	soilhum_decomp(:) + soilhum(:,l) * rpc(:) * &
599	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
600
601	ENDDO
602
603	control_moist(:,ibelow) = control_moist_func (npts, soilhum_decomp)
604
605	!! 5. fluxes from litter to carbon pools and respiration
606
607	DO l = 1, nlevs !Loop over litter levels (above and below ground)
608	DO m = 2,nvm !Loop over PFTs
609
610	!! 5.1 structural litter: goes into active and slow carbon pools + respiration
611
612	!! 5.1.1 total quantity of structural litter which is decomposed
613	fd(:) = dt/litter_tau(istructural) * &
614	control_temp(:,l) * control_moist(:,l) * exp( -litter_struct_coef * lignin_struc(:,m,l) )
615
616	DO k = 1,nelements
617
618	qd(:,k) = litter(:,istructural,m,l,k) * fd(:)
619
620	END DO
621
622	litter(:,istructural,m,l,:) = litter(:,istructural,m,l,:) - qd(:,:)
623
624	!! 5.1.2 decompose same fraction of structural part of dead leaves. Not exact
625	!! as lignine content is not the same as that of the total structural litter.
626	! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
627	! we do this test to do this calcul only ones in 1,nlev loop
628	if (l == iabove) dead_leaves(:,m,istructural) = dead_leaves(:,m,istructural) * ( un - fd(:) )
629
630	!! 5.1.3 non-lignin fraction of structural litter goes into
631	!! active carbon pool + respiration
632	soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
633	frac_soil(istructural,iactive,l) * qd(:,icarbon) * ( 1. - lignin_struc(:,m,l) ) / dt
634
635	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
636	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
637	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
638	! time step and avoid some constructions that could create bug during future developments.
639	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
640	( 1. - frac_soil(istructural,iactive,l) ) * qd(:,icarbon) * &
641	( 1. - lignin_struc(:,m,l) ) / dt
642
643	!! 5.1.4 lignin fraction of structural litter goes into
644	!! slow carbon pool + respiration
645	soilcarbon_input(:,islow,m) = soilcarbon_input(:,islow,m) + &
646	frac_soil(istructural,islow,l) * qd(:,icarbon) * lignin_struc(:,m,l) / dt
647
648	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
649	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
650	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
651	! time step and avoid some constructions that could create bug during future developments.
652	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
653	( 1. - frac_soil(istructural,islow,l) ) * qd(:,icarbon) * lignin_struc(:,m,l) / dt
654
655
656	!! 5.2 metabolic litter goes into active carbon pool + respiration
657
658	!! 5.2.1 total quantity of metabolic litter that is decomposed
659	fd(:) = dt/litter_tau(imetabolic) * control_temp(:,l) * control_moist(:,l)
660
661	DO k = 1,nelements
662
663	qd(:,k) = litter(:,imetabolic,m,l,k) * fd(:)
664
665	END DO
666
667	litter(:,imetabolic,m,l,:) = litter(:,imetabolic,m,l,:) - qd(:,:)
668
669	!! 5.2.2 decompose same fraction of metabolic part of dead leaves.
670	! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
671	! we do this test to do this calcul only ones in 1,nlev loop
672	if (l == iabove) dead_leaves(:,m,imetabolic) = dead_leaves(:,m,imetabolic) * ( 1. - fd(:) )
673
674
675	!! 5.2.3 put decomposed litter into carbon pool + respiration
676	soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
677	frac_soil(imetabolic,iactive,l) * qd(:,icarbon) / dt
678
679	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
680	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
681	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
682	! time step and avoid some constructions that could create bug during future developments.
683	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
684	( 1. - frac_soil(imetabolic,iactive,l) ) * qd(:,icarbon) / dt
685
686	ENDDO
687	ENDDO
688
689
690	!! 6. calculate fraction of total soil covered by dead leaves
691
692	CALL deadleaf (npts, veget_max, dead_leaves, deadleaf_cover)
693
694	!! 7. (Quasi-)Analytical Spin-up : Start filling MatrixA
695
696	IF (spinup_analytic) THEN
697
698	MatrixA(:,:,:,:) = zero
699	VectorB(:,:,:) = zero
700
701
702	DO m = 2,nvm
703
704	!- MatrixA : carbon fluxes leaving the litter
705
706	MatrixA(:,m,istructural_above,istructural_above)= - dt/litter_tau(istructural) * &
707	control_temp(:,iabove) * control_moist(:,iabove) * exp( -litter_struct_coef * lignin_struc(:,m,iabove) )
708
709	MatrixA(:,m,istructural_below,istructural_below) = - dt/litter_tau(istructural) * &
710	control_temp(:,ibelow) * control_moist(:,ibelow) * exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )
711
712	MatrixA(:,m,imetabolic_above,imetabolic_above) = - dt/litter_tau(imetabolic) * &
713	control_temp(:,iabove) * control_moist(:,iabove)
714
715	MatrixA(:,m,imetabolic_below,imetabolic_below) = - dt/litter_tau(imetabolic) * &
716	control_temp(:,ibelow) * control_moist(:,ibelow)
717
718
719	!- MatrixA : carbon fluxes between the litter and the pools (the rest of the matrix is filled in stomate_soilcarbon.f90)
720
721	MatrixA(:,m,iactive_pool,istructural_above) = frac_soil(istructural,iactive,iabove) * &
722	dt/litter_tau(istructural) * &
723	control_temp(:,iabove) * control_moist(:,iabove) * &
724	exp( -litter_struct_coef * lignin_struc(:,m,iabove) ) * &
725	( 1. - lignin_struc(:,m,iabove) )
726
727
728	MatrixA(:,m,iactive_pool,istructural_below) = frac_soil(istructural,iactive,ibelow) * &
729	dt/litter_tau(istructural) * &
730	control_temp(:,ibelow) * control_moist(:,ibelow) * &
731	exp( -litter_struct_coef * lignin_struc(:,m,ibelow) ) * &
732	( 1. - lignin_struc(:,m,ibelow) )
733
734
735	MatrixA(:,m,iactive_pool,imetabolic_above) = frac_soil(imetabolic,iactive,iabove) * &
736	dt/litter_tau(imetabolic) * control_temp(:,iabove) * control_moist(:,iabove)
737
738
739	MatrixA(:,m,iactive_pool,imetabolic_below) = frac_soil(imetabolic,iactive,ibelow) * &
740	dt/litter_tau(imetabolic) * control_temp(:,ibelow) * control_moist(:,ibelow)
741
742	MatrixA(:,m,islow_pool,istructural_above) = frac_soil(istructural,islow,iabove) * &
743	dt/litter_tau(istructural) * &
744	control_temp(:,iabove) * control_moist(:,iabove) * &
745	exp( -litter_struct_coef * lignin_struc(:,m,iabove) )* &
746	lignin_struc(:,m,iabove)
747
748
749	MatrixA(:,m,islow_pool,istructural_below) = frac_soil(istructural,islow,ibelow) * &
750	dt/litter_tau(istructural) * &
751	control_temp(:,ibelow) * control_moist(:,ibelow) * &
752	exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )* &
753	lignin_struc(:,m,ibelow)
754
755
756	!- VectorB : carbon input -
757
758	VectorB(:,m,istructural_above) = litter_inc_PFT(:,m,istructural,iabove,icarbon)
759	VectorB(:,m,istructural_below) = litter_inc_PFT(:,m,istructural,ibelow,icarbon)
760	VectorB(:,m,imetabolic_above) = litter_inc_PFT(:,m,imetabolic,iabove,icarbon)
761	VectorB(:,m,imetabolic_below) = litter_inc_PFT(:,m,imetabolic,ibelow,icarbon)
762
763	IF (bavard.GE.4) WRITE(numout,*) 'We filled MatrixA and VectorB'
764
765	ENDDO ! Loop over # PFTs
766
767	ENDIF ! spinup analytic
768
769	IF (bavard.GE.4) WRITE(numout,*) 'Leaving littercalc'
770
771	END SUBROUTINE littercalc
772
773
774	!! ==============================================================================================================================\n
775	!! SUBROUTINE : deadleaf
776	!!
777	!>\BRIEF This routine calculates the deadleafcover.
778	!!
779	!! DESCRIPTION : It first calculates the lai corresponding to the dead leaves (LAI) using
780	!! the dead leaves carbon content (DL) the specific leaf area (sla) and the
781	!! maximal coverage fraction of a PFT (vegetmax) using the following equations:
782	!! \latexonly
783	!! \input{deadleaf1.tex}
784	!! \endlatexonly
785	!! \n
786	!! Then, the dead leaf cover (DLC) is calculated as following:\n
787	!! \latexonly
788	!! \input{deadleaf2.tex}
789	!! \endlatexonly
790	!! \n
791	!!
792	!! RECENT CHANGE(S) : None
793	!!
794	!! MAIN OUTPUT VARIABLE: ::deadleaf_cover
795	!!
796	!! REFERENCE(S) : None
797	!!
798	!! FLOWCHART : None
799	!! \n
800	!_ ================================================================================================================================
801
802	SUBROUTINE deadleaf (npts, veget_max, dead_leaves, deadleaf_cover)
803
804	!! 0. Variable and parameter declaration
805
806	!! 0.1 Input variables
807
808	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels (unitless)
809	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
810	!! metabolic and structural
811	!! @tex $(gC m^{-2})$ @endtex
812	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT defined in
813	!! the input vegetation map
814	!! @tex $(m^2 m^{-2})$ @endtex
815
816	!! 0.2 Output variables
817
818	REAL(r_std), DIMENSION(npts), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves over all PFTs
819	!! (0-1, unitless)
820
821	!! 0.3 Modified variables
822
823	!! 0.4 Local variables
824
825	REAL(r_std), DIMENSION(npts) :: dead_lai !! LAI of dead leaves @tex $(m^2 m^{-2})$ @endtex
826	INTEGER(i_std) :: j !! Index (unitless)
827	!_ ================================================================================================================================
828
829	!! 1. LAI of dead leaves
830
831	dead_lai(:) = zero
832
833	DO j = 2,nvm !Loop over PFTs
834	dead_lai(:) = dead_lai(:) + ( dead_leaves(:,j,imetabolic) + dead_leaves(:,j,istructural) ) * sla(j) &
835	* veget_max(:,j)
836	ENDDO
837
838	!! 2. fraction of soil covered by dead leaves
839
840	deadleaf_cover(:) = un - exp( - 0.5 * dead_lai(:) )
841
842	IF (bavard.GE.4) WRITE(numout,*) 'Leaving deadleaf'
843
844	END SUBROUTINE deadleaf
845
846
847	!! ================================================================================================================================
848	!! FUNCTION : control_moist_func
849	!!
850	!>\BRIEF Calculate moisture control for litter and soild C decomposition
851	!!
852	!! DESCRIPTION : Calculate moisture control factor applied
853	!! to litter decomposition and to soil carbon decomposition in
854	!! stomate_soilcarbon.f90 using the following equation: \n
855	!! \latexonly
856	!! \input{control_moist_func1.tex}
857	!! \endlatexonly
858	!! \n
859	!! with M the moisture control factor and soilmoisutre, the soil moisture
860	!! calculated in sechiba.
861	!! Then, the function is ranged between 0.25 and 1:\n
862	!! \latexonly
863	!! \input{control_moist_func2.tex}
864	!! \endlatexonly
865	!! \n
866	!! RECENT CHANGE(S) : None
867	!!
868	!! RETURN VALUE : ::moistfunc_result
869	!!
870	!! REFERENCE(S) : None
871	!!
872	!! FLOWCHART : None
873	!! \n
874	!_ ================================================================================================================================
875
876	FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)
877
878	!! 0. Variable and parameter declaration
879
880	!! 0.1 Input variables
881
882	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixel (unitless)
883	REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in !! relative humidity (unitless)
884
885	!! 0.2 Output variables
886
887	REAL(r_std), DIMENSION(npts) :: moistfunc_result !! Moisture control factor (0.25-1, unitless)
888
889	!! 0.3 Modified variables
890
891	!! 0.4 Local variables
892
893	!_ ================================================================================================================================
894
895	moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
896	moistfunc_result(:) = MAX( 0.25_r_std, MIN( un, moistfunc_result(:) ) )
897
898	END FUNCTION control_moist_func
899
900
901	!! ================================================================================================================================
902	!! FUNCTION : control_temp_func
903	!!
904	!>\BRIEF Calculate temperature control for litter and soild C decomposition
905	!!
906	!! DESCRIPTION : Calculate temperature control factor applied
907	!! to litter decomposition and to soil carbon decomposition in
908	!! stomate_soilcarbon.f90 using the following equation: \n
909	!! \latexonly
910	!! \input{control_temp_func1.tex}
911	!! \endlatexonly
912	!! \n
913	!! with T the temperature control factor, temp the temperature in Kelvin of
914	!! the air (for aboveground litter) or of the soil (for belowground litter
915	!! and soil)
916	!! Then, the function is limited in its maximal range to 1:\n
917	!! \latexonly
918	!! \input{control_temp_func2.tex}
919	!! \endlatexonly
920	!! \n
921	!! RECENT CHANGE(S) : None
922	!!
923	!! RETURN VALUE: ::tempfunc_result
924	!!
925	!! REFERENCE(S) : None
926	!!
927	!! FLOWCHART : None
928	!! \n
929	!_ ================================================================================================================================
930
931	FUNCTION control_temp_func (npts, temp_in) RESULT (tempfunc_result)
932
933	!! 0. Variable and parameter declaration
934
935	!! 0.1 Input variables
936	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of land pixels (unitless)
937	REAL(r_std), DIMENSION(npts), INTENT(in) :: temp_in !! Temperature (K)
938
939	!! 0.2 Output variables
940	REAL(r_std), DIMENSION(npts) :: tempfunc_result !! Temperature control factor (0-1, unitless)
941
942	!! 0.3 Modified variables
943
944	!! 0.4 Local variables
945
946	!_ ================================================================================================================================
947
948	tempfunc_result(:) = exp( soil_Q10 * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
949	tempfunc_result(:) = MIN( un, tempfunc_result(:) )
950
951	END FUNCTION control_temp_func
952
953	END MODULE stomate_litter

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

Download in other formats: