Context Navigation

stomate_litter.f90 @ 8579

Last change on this file since 8579 was 8462, checked in by josefine.ghattas, 4 months ago
The Moyano function describing the soil moisture effect on OM decomposition is added. It has been developed by Elodie Salmon in another branch and integrated in ORCHIDEE_2_2 by Bertrad Guenet. This commit corresponds to a corrected version of [8418].
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 54.9 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_litter
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Update litter and lignine content after litter fall and
10	!! calculating litter decomposition.
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL$
20	!! $Date$
21	!! $Revision$
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE stomate_litter
26
27	! modules used:
28	USE xios_orchidee
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_litter = .TRUE. !! first call
43	!$OMP THREADPRIVATE(firstcall_litter)
44
45	CONTAINS
46
47	!! ================================================================================================================================
48	!! SUBROUTINE : littercalc_calc
49	!!
50	!!\BRIEF Set the flag ::firstcall_litter 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_litter =.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, &
180	turnover, bm_to_litter, &
181	veget_cov_max, tsurf, tsoil, soilhum, litterhum, &
182	soilhumSAT, litterhumSAT, &
183	litterpart, litter, dead_leaves, lignin_struc, &
184	deadleaf_cover, resp_hetero_litter, &
185	soilcarbon_input, control_temp, control_moist, &
186	MatrixA, VectorB, bulk, clay, carbon)
187
188	!! 0. Variable and parameter declaration
189
190	!! 0.1 Input variables
191
192	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels
193	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: turnover !! Turnover rates of plant biomass
194	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
195	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: bm_to_litter !! Conversion of biomass to litter
196	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
197	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_cov_max !! PFT "Maximal" coverage fraction of a PFT
198	!! defined in the input vegetation map
199	!! @tex $(m^2 m^{-2})$ @endtex
200	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf !! Temperature (K) at the surface
201	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: tsoil !! Soil temperature (K)
202	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: soilhum !! Daily soil humidity of each soil layer
203	!! (unitless)
204	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: soilhumSAT !! Daily soil humidity of each soil layer
205	!! (unitless)
206	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhum !! Daily litter humidity (unitless)
207	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhumSAT !! Daily litter humidity (unitless)
208	REAL(r_std), DIMENSION(npts), INTENT(in) :: clay !! Clay fraction (unitless, 0-1)
209	REAL(r_std), DIMENSION(npts),INTENT(in) :: bulk !! Bulk density (kg/m**3)
210	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in) :: carbon !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f
211
212	!! 0.2 Output variables
213
214	REAL(r_std), DIMENSION(npts), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves
215	!! over all PFTs (0-1, unitless)
216	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_hetero_litter !! Litter heterotrophic respiration. The unit
217	!! is given by m^2 of ground.
218	!! @tex $(gC dt_sechiba one\_day^{-1}) m^{-2})$ @endtex
219	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(out) :: soilcarbon_input !! Quantity of carbon going into carbon pools
220	!! from litter decomposition. The unit is
221	!! given by m^2 of ground
222	!! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex
223	REAL(r_std), DIMENSION(npts,nlevs), INTENT(out) :: control_temp !! Temperature control of heterotrophic
224	!! respiration, above and below (0-1,
225	!! unitless)
226	REAL(r_std), DIMENSION(npts,nlevs), INTENT(out) :: control_moist !! Moisture control of heterotrophic
227	!! respiration (0.25-1, unitless)
228	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(out) :: MatrixA !! Matrix containing the fluxes between the
229	!! carbon pools per sechiba time step
230	!! @tex $(gC.m^2.day^{-1})$ @endtex
231	REAL(r_std), DIMENSION(npts,nvm,nbpools), INTENT(out) :: VectorB !! Vector containing the litter increase per
232	!! sechiba time step
233	!! @tex $(gC m^{-2})$ @endtex
234
235	!! 0.3 Modified variables
236
237	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: litterpart !! Fraction of litter above the ground
238	!! belonging to the different PFTs (0-1,
239	!! unitless)
240	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter !! Metabolic and structural litter,above and
241	!! below ground. The unit is given by m^2 of
242	!! ground @tex $(gC m^{-2})$ @endtex
243	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
244	!! metabolic and structural in
245	!! @tex $(gC m^{-2})$ @endtex
246	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_struc !! Ratio Lignin content in structural litter,
247	!! above and below ground, (0-1, unitless)
248
249	!! 0.4 Local variables
250
251	REAL(r_std) :: dt !! Number of sechiba(fast processes) time-step per day
252	REAL(r_std), SAVE, DIMENSION(nparts,nlitt) :: litterfrac !! The fraction of leaves, wood, etc. that
253	!! goes into metabolic and structural
254	!! litterpools (0-1, unitless)
255	!$OMP THREADPRIVATE(litterfrac)
256	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: z_soil !! Soil levels (m)
257	!$OMP THREADPRIVATE(z_soil)
258	REAL(r_std), DIMENSION(npts) :: rpc !! Integration constant for vertical root
259	!! profiles (unitless)
260	REAL(r_std), SAVE, DIMENSION(nlitt) :: litter_tau !! Turnover time in litter pools (days)
261	!$OMP THREADPRIVATE(litter_tau)
262	REAL(r_std), SAVE, DIMENSION(nlitt,ncarb,nlevs) :: frac_soil !! Fraction of litter that goes into soil
263	!! (litter -> carbon, above and below). The
264	!! remaining part goes to the atmosphere
265	!$OMP THREADPRIVATE(frac_soil)
266	REAL(r_std), DIMENSION(npts) :: tsoil_decomp !! Temperature used for decompostition in
267	!! soil (K)
268	REAL(r_std), DIMENSION(npts) :: soilhum_decomp !! Humidity used for decompostition in soil
269	!! (unitless)
270	REAL(r_std), DIMENSION(npts) :: soilhumSAT_decomp !! Humidity used for decompostition in soil
271	!! (unitless) for Moyano et al 2012
272	REAL(r_std), DIMENSION(npts) :: fd !! Fraction of structural or metabolic litter
273	!! decomposed (unitless)
274	REAL(r_std), DIMENSION(npts,nelements) :: qd !! Quantity of structural or metabolic litter
275	!! decomposed @tex $(gC m^{-2})$ @endtex
276	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: old_struc !! Old structural litter, above and below
277	!! @tex $(gC m^{-2})$ @endtex
278	REAL(r_std), DIMENSION(npts,nvm,nlitt,nlevs,nelements) :: litter_inc_PFT !! Increase of litter, per PFT, metabolic and
279	!! structural, above and below ground. The
280	!! unit is given by m^2 of ground.
281	!! @tex $(gC m^{-2})$ @endtex
282	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements) :: litter_inc !! Increase of metabolic and structural
283	!! litter, above and below ground. The unit
284	!! is given by m^2 of ground.
285	!! @tex $(gC m^{-2})$ @endtex
286	REAL(r_std), DIMENSION(npts,nvm,nlevs) :: lignin_struc_inc !! Lignin increase in structural litter,
287	!! above and below ground. The unit is given
288	!! by m^2 of ground.
289	!! @tex $(gC m^{-2})$ @endtex
290	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements) :: litter_pft !! Metabolic and structural litter above the
291	!! ground per PFT
292	REAL(r_std), DIMENSION(npts) :: zdiff_min !! Intermediate field for looking for minimum
293	!! of what? this is not used in the code.
294	!! [??CHECK] could we delete it?
295	CHARACTER(LEN=10), DIMENSION(nlitt) :: litter_str !! Messages to write output information about
296	!! the litter
297	CHARACTER(LEN=22), DIMENSION(nparts) :: part_str !! Messages to write output information about
298	!! the plant
299	CHARACTER(LEN=7), DIMENSION(ncarb) :: carbon_str !! Messages to write output information about
300	!! the soil carbon
301	CHARACTER(LEN=5), DIMENSION(nlevs) :: level_str !! Messages to write output information about
302	!! the level (aboveground or belowground litter)
303	INTEGER(i_std) :: i,j,k,l,m !! Indices (unitless)
304	INTEGER(i_std) :: ier !! Error handling
305	!_ ================================================================================================================================
306
307	IF (printlev>=3) WRITE(numout,*) 'Entering littercalc'
308
309	!! 1. Initialisations of the different fields during the first call of the routine
310	dt = dt_sechiba/one_day
311	IF ( firstcall_litter ) THEN
312
313	!! 1.1.3 litter fractions:
314	!! what fraction of leaves, wood, etc. goes into metabolic and structural litterpools
315	DO k = 1, nparts
316
317	litterfrac(k,imetabolic) = metabolic_ref_frac - metabolic_LN_ratio * LC(k) * CN(k)
318	litterfrac(k,istructural) = un - litterfrac(k,imetabolic)
319
320	ENDDO
321
322	!! 1.1.4 residence times in litter pools (days)
323	litter_tau(imetabolic) = tau_metabolic * one_year ! .5 years
324	litter_tau(istructural) = tau_struct * one_year ! 3 years
325
326	!! 1.1.5 decomposition flux fraction that goes into soil
327	! (litter -> carbon, above and below)
328	! 1-frac_soil goes into atmosphere
329	frac_soil(:,:,:) = zero
330
331	! structural litter: lignin fraction goes into slow pool + respiration,
332	! rest into active pool + respiration
333	frac_soil(istructural,iactive,iabove) = frac_soil_struct_aa
334	frac_soil(istructural,iactive,ibelow) = frac_soil_struct_ab
335	frac_soil(istructural,islow,iabove) = frac_soil_struct_sa
336	frac_soil(istructural,islow,ibelow) = frac_soil_struct_sb
337
338	! metabolic litter: all goes into active pool + respiration.
339	! Nothing into slow or passive pool.
340	frac_soil(imetabolic,iactive,iabove) = frac_soil_metab_aa
341	frac_soil(imetabolic,iactive,ibelow) = frac_soil_metab_ab
342
343
344	!! 1.2 soil levels
345	ALLOCATE(z_soil(0:nslm), stat=ier)
346	IF ( ier /= 0 ) CALL ipslerr_p(3,'littercalc','Pb in allocate of z_soil','','')
347
348	z_soil(0) = zero
349	z_soil(1:nslm) = diaglev(1:nslm)
350
351
352	!! 1.3 messages
353	litter_str(imetabolic) = 'metabolic'
354	litter_str(istructural) = 'structural'
355
356	carbon_str(iactive) = 'active'
357	carbon_str(islow) = 'slow'
358	carbon_str(ipassive) = 'passive'
359
360	level_str(iabove) = 'above'
361	level_str(ibelow) = 'below'
362
363	part_str(ileaf) = 'leaves'
364	part_str(isapabove) = 'sap above ground'
365	part_str(isapbelow) = 'sap below ground'
366	part_str(iheartabove) = 'heartwood above ground'
367	part_str(iheartbelow) = 'heartwood below ground'
368	part_str(iroot) = 'roots'
369	part_str(ifruit) = 'fruits'
370	part_str(icarbres) = 'carbohydrate reserve'
371
372	IF (printlev >=2) THEN
373	WRITE(numout,*) 'litter:'
374
375	WRITE(numout,*) ' > C/N ratios: '
376	DO k = 1, nparts
377	WRITE(numout,*) ' ', part_str(k), ': ',CN(k)
378	ENDDO
379
380	WRITE(numout,*) ' > Lignine/C ratios: '
381	DO k = 1, nparts
382	WRITE(numout,*) ' ', part_str(k), ': ',LC(k)
383	ENDDO
384
385	WRITE(numout,*) ' > fraction of compartment that goes into litter: '
386	DO k = 1, nparts
387	DO m = 1, nlitt
388	WRITE(numout,*) ' ', part_str(k), '-> ',litter_str(m), ':',litterfrac(k,m)
389	ENDDO
390	ENDDO
391
392	WRITE(numout,*) ' > scaling depth for decomposition (m): ',z_decomp
393
394	WRITE(numout,*) ' > minimal carbon residence time in litter pools (d):'
395	DO m = 1, nlitt
396	WRITE(numout,*) ' ',litter_str(m),':',litter_tau(m)
397	ENDDO
398
399	WRITE(numout,*) ' > litter decomposition flux fraction that really goes '
400	WRITE(numout,*) ' into carbon pools (rest into the atmosphere):'
401	DO m = 1, nlitt
402	DO l = 1, nlevs
403	DO k = 1, ncarb
404	WRITE(numout,*) ' ',litter_str(m),' ',level_str(l),' -> ',&
405	carbon_str(k),':', frac_soil(m,k,l)
406	ENDDO
407	ENDDO
408	ENDDO
409	END IF ! printlev >=2
410	firstcall_litter = .FALSE.
411
412	ENDIF
413
414
415	!! 1.3 litter above the ground per PFT.
416	DO j = 1, nvm ! Loop over # PFTs
417
418	DO k = 1, nlitt !Loop over litter pool
419	litter_pft(:,j,k,icarbon) = litterpart(:,j,k) * litter(:,k,j,iabove,icarbon)
420	ENDDO
421
422	ENDDO
423
424
425	!! 1.4 set output to zero
426	deadleaf_cover(:) = zero
427	resp_hetero_litter(:,:) = zero
428	soilcarbon_input(:,:,:) = zero
429
430
431	!! 2. Add biomass to different litterpools (per m^2 of ground)
432
433	!! 2.1 first, save old structural litter (needed for lignin fractions).
434	! above/below
435	DO l = 1, nlevs !Loop over litter levels (above and below ground)
436	DO m = 1,nvm !Loop over PFTs
437
438	old_struc(:,m,l) = litter(:,istructural,m,l,icarbon)
439
440	ENDDO
441	ENDDO
442
443
444	!! 2.2 update litter, dead leaves, and lignin content in structural litter
445	litter_inc(:,:,:,:,:) = zero
446	lignin_struc_inc(:,:,:) = zero
447
448	DO j = 1,nvm !Loop over PFTs
449
450	!! 2.2.1 litter
451	DO k = 1, nlitt !Loop over litter pools (metabolic and structural)
452
453	!! 2.2.2 calculate litter increase (per m^2 of ground).
454	! Only a given fracion of fruit turnover is directly coverted into litter.
455	! Litter increase for each PFT, structural and metabolic, above/below
456	litter_inc_PFT(:,j,k,iabove,icarbon) = &
457	litterfrac(ileaf,k) * bm_to_litter(:,j,ileaf,icarbon) + &
458	litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
459	litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
460	litterfrac(ifruit,k) * bm_to_litter(:,j,ifruit,icarbon) + &
461	litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
462	litterfrac(ileaf,k) * turnover(:,j,ileaf,icarbon) + &
463	litterfrac(isapabove,k) * turnover(:,j,isapabove,icarbon) + &
464	litterfrac(iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
465	litterfrac(ifruit,k) * turnover(:,j,ifruit,icarbon) + &
466	litterfrac(icarbres,k) * turnover(:,j,icarbres,icarbon)
467
468	litter_inc_PFT(:,j,k,ibelow,icarbon) = &
469	litterfrac(isapbelow,k) * bm_to_litter(:,j,isapbelow,icarbon) + &
470	litterfrac(iheartbelow,k) * bm_to_litter(:,j,iheartbelow,icarbon) + &
471	litterfrac(iroot,k) * bm_to_litter(:,j,iroot,icarbon) + &
472	litterfrac(isapbelow,k) * turnover(:,j,isapbelow,icarbon) + &
473	litterfrac(iheartbelow,k) * turnover(:,j,iheartbelow,icarbon) + &
474	litterfrac(iroot,k) * turnover(:,j,iroot,icarbon)
475
476	! litter increase, met/struct, above/below
477	litter_inc(:,k,j,iabove,icarbon) = litter_inc(:,k,j,iabove,icarbon) + litter_inc_PFT(:,j,k,iabove,icarbon)
478	litter_inc(:,k,j,ibelow,icarbon) = litter_inc(:,k,j,ibelow,icarbon) + litter_inc_PFT(:,j,k,ibelow,icarbon)
479
480	!! 2.2.3 dead leaves, for soil cover.
481	dead_leaves(:,j,k) = &
482	dead_leaves(:,j,k) + &
483	litterfrac(ileaf,k) * ( bm_to_litter(:,j,ileaf,icarbon) + turnover(:,j,ileaf,icarbon) )
484
485	!! 2.2.4 lignin increase in structural litter
486	IF ( k .EQ. istructural ) THEN
487
488	lignin_struc_inc(:,j,iabove) = &
489	lignin_struc_inc(:,j,iabove) + &
490	LC(ileaf) * bm_to_litter(:,j,ileaf,icarbon) + &
491	LC(isapabove) * bm_to_litter(:,j,isapabove,icarbon) + &
492	LC(iheartabove) * bm_to_litter(:,j,iheartabove,icarbon) + &
493	LC(ifruit) * bm_to_litter(:,j,ifruit,icarbon) + &
494	LC(icarbres) * bm_to_litter(:,j,icarbres,icarbon) + &
495	LC(ileaf) * turnover(:,j,ileaf,icarbon) + &
496	LC(isapabove) * turnover(:,j,isapabove,icarbon) + &
497	LC(iheartabove) * turnover(:,j,iheartabove,icarbon) + &
498	LC(ifruit) * turnover(:,j,ifruit,icarbon) + &
499	LC(icarbres) * turnover(:,j,icarbres,icarbon)
500
501	lignin_struc_inc(:,j,ibelow) = &
502	lignin_struc_inc(:,j,ibelow) + &
503	LC(isapbelow) * bm_to_litter(:,j,isapbelow,icarbon) + &
504	LC(iheartbelow) * bm_to_litter(:,j,iheartbelow,icarbon) + &
505	LC(iroot) * bm_to_litter(:,j,iroot,icarbon) + &
506	LC(isapbelow)*turnover(:,j,isapbelow,icarbon) + &
507	LC(iheartbelow)*turnover(:,j,iheartbelow,icarbon) + &
508	LC(iroot)*turnover(:,j,iroot,icarbon)
509
510	ENDIF
511
512	ENDDO
513	ENDDO
514
515	!! 2.2.5 add new litter (struct/met, above/below)
516	litter(:,:,:,:,:) = litter(:,:,:,:,:) + litter_inc(:,:,:,:,:)
517
518	!! 2.2.6 for security: can't add more lignin than structural litter (above/below)
519	DO l = 1, nlevs !Loop over litter levels (above and below ground)
520	DO m = 1,nvm !Lopp over PFTs
521
522	lignin_struc_inc(:,m,l) = &
523	MIN( lignin_struc_inc(:,m,l), litter_inc(:,istructural,m,l,icarbon) )
524
525	ENDDO
526	ENDDO
527
528	!! 2.2.7 new lignin content: add old lignin and lignin increase, divide by
529	!! total structural litter (above/below)
530	DO l = 1, nlevs !Loop over litter levels (above and below ground)
531	DO m = 1,nvm !Loop over PFTs
532	WHERE( litter(:,istructural,m,l,icarbon) .GT. min_stomate )
533
534	!MM : Soenke modif
535	! Best vectorization ?
536	!!$ lignin_struc(:,:,:) = &
537	!!$ ( lignin_struc(:,:,:)*old_struc(:,:,:) + lignin_struc_inc(:,:,:) ) / &
538	!!$ litter(:,istructural,:,:,icarbon)
539
540	lignin_struc(:,m,l) = lignin_struc(:,m,l) * old_struc(:,m,l)
541	lignin_struc(:,m,l) = lignin_struc(:,m,l) + lignin_struc_inc(:,m,l)
542	lignin_struc(:,m,l) = lignin_struc(:,m,l) / litter(:,istructural,m,l,icarbon)
543	ELSEWHERE
544	lignin_struc(:,m,l) = zero
545	ENDWHERE
546	ENDDO
547	ENDDO
548
549
550	!! 2.3 new litter fraction per PFT (for structural and metabolic litter, above
551	!! the ground).
552	DO j = 1,nvm !Loop over PFTs
553
554	WHERE ( litter(:,:,j,iabove,icarbon) .GT. min_stomate )
555
556	litterpart(:,j,:) = &
557	( litter_pft(:,j,:,icarbon) + litter_inc_PFT(:,j,:,iabove,icarbon) ) / litter(:,:,j,iabove,icarbon)
558
559	ELSEWHERE
560
561	litterpart(:,j,:) = zero
562
563	ENDWHERE
564
565	ENDDO
566
567
568	!! 3. Temperature control on decay: Factor between 0 and 1
569
570	!! 3.1 above: surface temperature
571	control_temp(:,iabove) = control_temp_func (npts, tsurf)
572
573
574	!! 3.2 below: convolution of temperature and decomposer profiles
575	!! (exponential decomposer profile supposed)
576
577	!! 3.2.1 rpc is an integration constant such that the integral of the root profile is 1.
578	rpc(:) = un / ( un - EXP( -z_soil(nslm) / z_decomp ) )
579
580	!! 3.2.2 integrate over the nslm levels
581	tsoil_decomp(:) = zero
582
583	DO l = 1, nslm
584
585	tsoil_decomp(:) = &
586	tsoil_decomp(:) + tsoil(:,l) * rpc(:) * &
587	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
588
589	ENDDO
590
591	control_temp(:,ibelow) = control_temp_func (npts, tsoil_decomp)
592
593	!! 4. Moisture control. Factor between 0 and 1
594
595	!! 4.1 above the ground: litter humidity
596
597	IF (ok_moyano_soilhumsat)THEN !!ok_moyano_soilhumsat is true
598	control_moist(:,iabove) = control_moist_func_moyano (npts,clay,bulk,carbon, litter ,litterhumSAT, veget_cov_max)
599	ELSE !ok_moyano_soilhumsat is false by default
600	control_moist(:,iabove) = control_moist_func (npts, litterhum)
601	ENDIF
602
603	CALL xios_orchidee_send_field("ControlMoistLitter",control_moist(:,iabove))
604	CALL xios_orchidee_send_field("litterhumSAT",litterhumSAT(:))
605	!
606	!! 4.2 below: convolution of humidity and decomposer profiles
607	! (exponential decomposer profile supposed)
608
609	!! 4.2.1 rpc is an integration constant such that the integral of the root profile is 1.
610	rpc(:) = un / ( un - EXP( -z_soil(nslm) / z_decomp ) )
611
612	!! 4.2.2 integrate over the nslm levels
613	soilhum_decomp(:) = zero
614	soilhumSAT_decomp(:) = zero
615
616	IF (ok_moyano_soilhumsat)THEN !!ok_moyano_soilhumsat is true
617	DO l = 1, nslm !Loop over soil levels
618	soilhumSAT_decomp(:) = &
619	soilhumSAT_decomp(:) + soilhumSAT(:,l) * rpc(:) * &
620	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
621	ENDDO
622	control_moist(:,ibelow) = control_moist_func_moyano (npts,clay,bulk,carbon, litter, soilhumSAT_decomp,veget_cov_max)
623
624	ELSE !ok_moyano_soilhumsat is false by default
625	DO l = 1, nslm !Loop over soil levels
626
627	soilhum_decomp(:) = &
628	soilhum_decomp(:) + soilhum(:,l) * rpc(:) * &
629	( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
630	ENDDO
631	control_moist(:,ibelow) = control_moist_func (npts, soilhum_decomp)
632	ENDIF
633
634	CALL xios_orchidee_send_field("ControlMoistSoil",control_moist(:,ibelow))
635	CALL xios_orchidee_send_field("soilhumSAT",soilhumSAT_decomp(:))
636
637	!! 5. fluxes from litter to carbon pools and respiration
638
639	DO l = 1, nlevs !Loop over litter levels (above and below ground)
640	DO m = 1,nvm !Loop over PFTs
641
642	!! 5.1 structural litter: goes into active and slow carbon pools + respiration
643
644	!! 5.1.1 total quantity of structural litter which is decomposed
645	fd(:) = dt/litter_tau(istructural) * &
646	control_temp(:,l) * control_moist(:,l) * exp( -litter_struct_coef * lignin_struc(:,m,l) )
647
648	DO k = 1,nelements
649
650	qd(:,k) = litter(:,istructural,m,l,k) * fd(:)
651
652	END DO
653
654	litter(:,istructural,m,l,:) = litter(:,istructural,m,l,:) - qd(:,:)
655
656	!! 5.1.2 decompose same fraction of structural part of dead leaves. Not exact
657	!! as lignine content is not the same as that of the total structural litter.
658	! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
659	! we do this test to do this calcul only ones in 1,nlev loop
660	if (l == iabove) dead_leaves(:,m,istructural) = dead_leaves(:,m,istructural) * ( un - fd(:) )
661
662	!! 5.1.3 non-lignin fraction of structural litter goes into
663	!! active carbon pool + respiration
664	soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
665	frac_soil(istructural,iactive,l) * qd(:,icarbon) * ( 1. - lignin_struc(:,m,l) ) / dt
666
667	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
668	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
669	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dt_sechiba/one_day
670	! time step and avoid some constructions that could create bug during future developments.
671	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
672	( 1. - frac_soil(istructural,iactive,l) ) * qd(:,icarbon) * &
673	( 1. - lignin_struc(:,m,l) ) / dt
674
675	!! 5.1.4 lignin fraction of structural litter goes into
676	!! slow carbon pool + respiration
677	soilcarbon_input(:,islow,m) = soilcarbon_input(:,islow,m) + &
678	frac_soil(istructural,islow,l) * qd(:,icarbon) * lignin_struc(:,m,l) / dt
679
680	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
681	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
682	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dt_sechiba/one_day
683	! time step and avoid some constructions that could create bug during future developments.
684	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
685	( 1. - frac_soil(istructural,islow,l) ) * qd(:,icarbon) * lignin_struc(:,m,l) / dt
686
687
688	!! 5.2 metabolic litter goes into active carbon pool + respiration
689
690	!! 5.2.1 total quantity of metabolic litter that is decomposed
691	fd(:) = dt/litter_tau(imetabolic) * control_temp(:,l) * control_moist(:,l)
692
693	DO k = 1,nelements
694
695	qd(:,k) = litter(:,imetabolic,m,l,k) * fd(:)
696
697	END DO
698
699	litter(:,imetabolic,m,l,:) = litter(:,imetabolic,m,l,:) - qd(:,:)
700
701	!! 5.2.2 decompose same fraction of metabolic part of dead leaves.
702	! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
703	! we do this test to do this calcul only ones in 1,nlev loop
704	if (l == iabove) dead_leaves(:,m,imetabolic) = dead_leaves(:,m,imetabolic) * ( 1. - fd(:) )
705
706
707	!! 5.2.3 put decomposed litter into carbon pool + respiration
708	soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
709	frac_soil(imetabolic,iactive,l) * qd(:,icarbon) / dt
710
711	!BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
712	! the sechiba time step but then in stomate.f90 resp_hetero_litter is multiplied by dt.
713	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dt_sechiba/one_day
714	! time step and avoid some constructions that could create bug during future developments.
715	resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
716	( 1. - frac_soil(imetabolic,iactive,l) ) * qd(:,icarbon) / dt
717
718	ENDDO
719	ENDDO
720
721
722	!! 6. calculate fraction of total soil covered by dead leaves
723
724	CALL deadleaf (npts, veget_cov_max, dead_leaves, deadleaf_cover)
725
726	!! 7. (Quasi-)Analytical Spin-up : Start filling MatrixA
727
728	IF (spinup_analytic) THEN
729
730	MatrixA(:,:,:,:) = zero
731	VectorB(:,:,:) = zero
732
733
734	DO m = 1,nvm
735
736	!- MatrixA : carbon fluxes leaving the litter
737
738	MatrixA(:,m,istructural_above,istructural_above)= - dt/litter_tau(istructural) * &
739	control_temp(:,iabove) * control_moist(:,iabove) * exp( -litter_struct_coef * lignin_struc(:,m,iabove) )
740
741	MatrixA(:,m,istructural_below,istructural_below) = - dt/litter_tau(istructural) * &
742	control_temp(:,ibelow) * control_moist(:,ibelow) * exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )
743
744	MatrixA(:,m,imetabolic_above,imetabolic_above) = - dt/litter_tau(imetabolic) * &
745	control_temp(:,iabove) * control_moist(:,iabove)
746
747	MatrixA(:,m,imetabolic_below,imetabolic_below) = - dt/litter_tau(imetabolic) * &
748	control_temp(:,ibelow) * control_moist(:,ibelow)
749
750
751	!- MatrixA : carbon fluxes between the litter and the pools (the rest of the matrix is filled in stomate_soilcarbon.f90)
752
753	MatrixA(:,m,iactive_pool,istructural_above) = frac_soil(istructural,iactive,iabove) * &
754	dt/litter_tau(istructural) * &
755	control_temp(:,iabove) * control_moist(:,iabove) * &
756	exp( -litter_struct_coef * lignin_struc(:,m,iabove) ) * &
757	( 1. - lignin_struc(:,m,iabove) )
758
759
760	MatrixA(:,m,iactive_pool,istructural_below) = frac_soil(istructural,iactive,ibelow) * &
761	dt/litter_tau(istructural) * &
762	control_temp(:,ibelow) * control_moist(:,ibelow) * &
763	exp( -litter_struct_coef * lignin_struc(:,m,ibelow) ) * &
764	( 1. - lignin_struc(:,m,ibelow) )
765
766
767	MatrixA(:,m,iactive_pool,imetabolic_above) = frac_soil(imetabolic,iactive,iabove) * &
768	dt/litter_tau(imetabolic) * control_temp(:,iabove) * control_moist(:,iabove)
769
770
771	MatrixA(:,m,iactive_pool,imetabolic_below) = frac_soil(imetabolic,iactive,ibelow) * &
772	dt/litter_tau(imetabolic) * control_temp(:,ibelow) * control_moist(:,ibelow)
773
774	MatrixA(:,m,islow_pool,istructural_above) = frac_soil(istructural,islow,iabove) * &
775	dt/litter_tau(istructural) * &
776	control_temp(:,iabove) * control_moist(:,iabove) * &
777	exp( -litter_struct_coef * lignin_struc(:,m,iabove) )* &
778	lignin_struc(:,m,iabove)
779
780
781	MatrixA(:,m,islow_pool,istructural_below) = frac_soil(istructural,islow,ibelow) * &
782	dt/litter_tau(istructural) * &
783	control_temp(:,ibelow) * control_moist(:,ibelow) * &
784	exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )* &
785	lignin_struc(:,m,ibelow)
786
787
788	!- VectorB : carbon input -
789
790	VectorB(:,m,istructural_above) = litter_inc_PFT(:,m,istructural,iabove,icarbon)
791	VectorB(:,m,istructural_below) = litter_inc_PFT(:,m,istructural,ibelow,icarbon)
792	VectorB(:,m,imetabolic_above) = litter_inc_PFT(:,m,imetabolic,iabove,icarbon)
793	VectorB(:,m,imetabolic_below) = litter_inc_PFT(:,m,imetabolic,ibelow,icarbon)
794
795	IF (printlev>=4) WRITE(numout,*) 'We filled MatrixA and VectorB'
796
797	ENDDO ! Loop over # PFTs
798
799	ENDIF ! spinup analytic
800
801	IF (printlev>=4) WRITE(numout,*) 'Leaving littercalc'
802
803	END SUBROUTINE littercalc
804
805
806	!! ==============================================================================================================================\n
807	!! SUBROUTINE : deadleaf
808	!!
809	!>\BRIEF This routine calculates the deadleafcover.
810	!!
811	!! DESCRIPTION : It first calculates the lai corresponding to the dead leaves (LAI) using
812	!! the dead leaves carbon content (DL) the specific leaf area (sla) and the
813	!! maximal coverage fraction of a PFT (vegetmax) using the following equations:
814	!! \latexonly
815	!! \input{deadleaf1.tex}
816	!! \endlatexonly
817	!! \n
818	!! Then, the dead leaf cover (DLC) is calculated as following:\n
819	!! \latexonly
820	!! \input{deadleaf2.tex}
821	!! \endlatexonly
822	!! \n
823	!!
824	!! RECENT CHANGE(S) : None
825	!!
826	!! MAIN OUTPUT VARIABLE: ::deadleaf_cover
827	!!
828	!! REFERENCE(S) : None
829	!!
830	!! FLOWCHART : None
831	!! \n
832	!_ ================================================================================================================================
833
834	SUBROUTINE deadleaf (npts, veget_cov_max, dead_leaves, deadleaf_cover)
835
836	!! 0. Variable and parameter declaration
837
838	!! 0.1 Input variables
839
840	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixels (unitless)
841	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: dead_leaves !! Dead leaves per ground unit area, per PFT,
842	!! metabolic and structural
843	!! @tex $(gC m^{-2})$ @endtex
844	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_cov_max !! PFT "Maximal" coverage fraction of a PFT defined in
845	!! the input vegetation map
846	!! @tex $(m^2 m^{-2})$ @endtex
847
848	!! 0.2 Output variables
849
850	REAL(r_std), DIMENSION(npts), INTENT(out) :: deadleaf_cover !! Fraction of soil covered by dead leaves over all PFTs
851	!! (0-1, unitless)
852
853	!! 0.3 Modified variables
854
855	!! 0.4 Local variables
856
857	REAL(r_std), DIMENSION(npts) :: dead_lai !! LAI of dead leaves @tex $(m^2 m^{-2})$ @endtex
858	INTEGER(i_std) :: j !! Index (unitless)
859	!_ ================================================================================================================================
860
861	!! 1. LAI of dead leaves
862
863	dead_lai(:) = zero
864
865	DO j = 1,nvm !Loop over PFTs
866	dead_lai(:) = dead_lai(:) + ( dead_leaves(:,j,imetabolic) + dead_leaves(:,j,istructural) ) * sla(j) &
867	* veget_cov_max(:,j)
868	ENDDO
869
870	!! 2. fraction of soil covered by dead leaves
871
872	deadleaf_cover(:) = un - exp( - 0.5 * dead_lai(:) )
873
874	IF (printlev>=4) WRITE(numout,*) 'Leaving deadleaf'
875
876	END SUBROUTINE deadleaf
877
878
879	!! ================================================================================================================================
880	!! FUNCTION : control_moist_func_moyano
881	!!
882	!>\BRIEF Calculate moisture control for litter and soil C decomposition
883	!!
884	!! DESCRIPTION : Calculate moisture control factor applied
885	!! to litter decomposition and to soil carbon decomposition in
886	!! stomate_soilcarbon.f90 using the following equation based on the Moyano et al., 2012 BG paper: \n
887	!! \latexonly
888	!! \input{control_moist_func1.tex}
889	!! \endlatexonly
890	!! \n
891	!! with M the moisture control factor and soilmoisutre, the soil moisture
892	!! calculated in sechiba.
893	!! Then, the function is ranged between Moistcont_min and 1:\n
894	!! \latexonly
895	!! \input{control_moist_func2.tex}
896	!! \endlatexonly
897	!! \n
898	!! RECENT CHANGE(S) : None
899	!!
900	!! RETURN VALUE : ::moistfunc_result
901	!!
902	!! REFERENCE(S) : None
903	!!
904	!! FLOWCHART : None
905	!! \n
906	!_ ================================================================================================================================
907
908	FUNCTION control_moist_func_moyano (npts,clay,bulk,carbon,litter,moist_in,veget_cov_max) RESULT (moistfunc_result)
909
910	!! 0. Variable and parameter declaration
911
912	!! 0.1 Input variables
913	INTEGER(i_std) :: k,m,n,i,l,e !!indices
914	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixel (unitless)
915	REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in !! relative humidity (unitless)
916	REAL(r_std), DIMENSION(npts), INTENT(in) :: clay !! Clay fraction (unitless, 0-1)
917	REAL(r_std), DIMENSION(npts), INTENT(in) :: bulk !! Bulk density (kg/m**3)
918	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in) :: carbon !! Soil carbon pools: active, s low, or passive, \f$(gC m^{2})$\f
919	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) ::litter !! Metabolic and structural litter,above and
920	!! below ground. The unit is given by m^2 of ground
921	!! @tex $(gC m^{-2})$ @endtex
922	REAL(r_std),DIMENSION(npts,nvm),INTENT(in):: veget_cov_max !! PFT "Maximal" coverage fraction of a PFT
923
924	!! 0.2 Output variables
925
926	REAL(r_std), DIMENSION(npts) :: moistfunc_result !! Moisture control factor (0.25-1, unitless)
927
928	!! 0.3 Modified variables
929
930	!! 0.4 Local variables
931	REAL(r_std), DIMENSION(101) :: Mint !! range [0,0.01,1] to compute PRSRmax
932	REAL(r_std), DIMENSION(npts,101) :: PRSR !! Proportional Response of Soil Respiration
933	REAL(r_std), DIMENSION(npts,101) :: SR !! Soil respiration compute using Moyano et al 2012 procedure
934	REAL(r_std) :: SRmax !! maximum value of SR
935	REAL(r_std), DIMENSION(npts,101) :: SRnorm !! SR normalized by the SRmax
936	REAL(r_std) :: ind !! i index of maximum SRnorm value
937	REAL(r_std) :: SRmin !! minimum value of SRnorm
938	REAL(r_std), DIMENSION(npts,101) :: SRsc !! Soil respiration rescale between 0 and 1 following Moyano et al 2012 procedure
939	REAL(r_std) :: SRscmax !! maximum value of SRsc
940	REAL(r_std) :: indmc !! i index of maximum mois_round values
941	REAL(r_std) :: moist_round !! moist_in round to 2decimal points
942	REAL(r_std), DIMENSION(npts) :: carbon_gCgs !! total Soil carbon pools: active +slow+ passive, (gC/gsoil)
943	REAL(r_std), DIMENSION(npts) :: carbon_gCm2s !! total Soil carbon pools: active +slow+ passive , (gC/gsoil)
944	REAL(r_std), DIMENSION(npts) :: clay_Moyano !! Clay fraction borned within the values used by Moyano et al., (unitless, 0-1)
945
946	!_ ===============================================================================================================================
947
948	!! In ok_moyano_soilhumsat total soil carbon per grid cells need to be computed and
949	!! converted from gC/m2soil to gC/gsoil
950	carbon_gCm2s(:)=zero
951	carbon_gCgs(:)=zero
952	DO n = 1, npts ! loop over grids
953	DO m = 1, nvm ! Loop over # PFTs
954	DO k = 1, ncarb ! Loop over carbon pools
955	carbon_gCm2s(n)=carbon_gCm2s(n)+carbon(n,k,m)*veget_cov_max(n,m)
956	ENDDO
957	ENDDO
958	ENDDO
959
960	IF (ok_orga) THEN
961	DO n = 1, npts ! loop over grids
962	carbon_gCgs(n)=carbon_gCm2s(n) / (bulk(n) * 1E3 * soilheight)
963	carbon_gCgs(n)= MAX(min_carbon_moyano,MIN(max_carbon_moyano,carbon_gCgs(n)))
964	clay_Moyano(n) = MAX(min_clay_moyano,MIN(max_clay_moyano,clay(n)))
965	ENDDO
966
967	!!1. compute Max(Prsr(M[0,0.01,1])); Prsr:Proportional Response of Soil
968	!!Respiration (PRSR) and Soil respiration(SR): SR(M)=SRo x Prsr(M)/Max(Prsr(M[0,0.01,1]))
969	!! M:soilhumSAT, SRo=1 (arbitrary)
970	DO n = 1,npts
971	Mint(:)= zero
972	SRmax = zero
973	DO i = 1,101
974	Mint(i)=0.01*(i-1)
975	IF (carbon_gCgs(n) .LT. limit_carbon_orga) THEN
976	PRSR(n,i) = beta1 * Mint(i) + beta2 * Mint(i)*2.0 + beta3 Mint(i)**3.0 &
977	& + beta4 * clay_Moyano(n) + beta5 * clay_Moyano(n) * Mint(i) &
978	& + beta6 * carbon_gCgs(n) + intercept
979	ELSE
980	PRSR(n,i) = beta1_orga * Mint(i) + beta2_orga * Mint(i)**2.0 + &
981	& beta3_orga * Mint(i)**3.0 + intercept_orga
982	ENDIF
983
984	IF (i.LT.2) THEN
985	SR(n,i) = PRSR(n,i)
986	SRmax = MAX(SR(n,i),SRmax)
987	ELSE
988	SR(n,i) = PRSR(n,i) * SR(n,i-1)
989	SRmax = MAX(SR(n,i),SRmax)
990	ENDIF
991	ENDDO
992	SRnorm(n,:)= SRo * SR(n,:)/SRmax
993	ENDDO
994	ELSE
995	DO n = 1, npts ! loop over grids
996	carbon_gCgs(n)=carbon_gCm2s(n) / (bulk(n) * 1E3 * soilheight)
997	carbon_gCgs(n)= MAX(min_carbon_moyano,MIN(limit_carbon_orga,carbon_gCgs(n)))
998	clay_Moyano(n) = MAX(min_clay_moyano,MIN(max_clay_moyano,clay(n)))
999	ENDDO
1000
1001	!!1. compute Max(Prsr(M[0,0.01,1])); Prsr:Proportional Response of Soil
1002	!!Respiration (PRSR) and Soil respiration(SR): SR(M)=SRo x
1003	!Prsr(M)/Max(Prsr(M[0,0.01,1]))
1004	!! M:soilhumSAT, SRo=1 (arbitrary)
1005	DO n = 1,npts
1006	Mint(:)= zero
1007	SRmax = zero
1008	DO i = 1,101
1009	Mint(i)=0.01*(i-1)
1010	PRSR(n,i) = beta1 * Mint(i) + beta2 * Mint(i)*2.0 + beta3 Mint(i)**3.0 &
1011	& + beta4 * clay_Moyano(n) + beta5 * clay_Moyano(n) * Mint(i) &
1012	& + beta6 * carbon_gCgs(n) + intercept
1013
1014	IF (i.LT.2) THEN
1015	SR(n,i) = PRSR(n,i)
1016	SRmax = MAX(SR(n,i),SRmax)
1017	ELSE
1018	SR(n,i) = PRSR(n,i) * SR(n,i-1)
1019	SRmax = MAX(SR(n,i),SRmax)
1020	ENDIF
1021	ENDDO
1022	SRnorm(n,:)= SRo * SR(n,:)/SRmax
1023	ENDDO
1024	ENDIF
1025
1026	!!2. Rescale SR values between 0 and 1 and defined SR for moist_in
1027	DO n = 1,npts
1028	ind = 1.0
1029	indmc = zero
1030	moist_round = zero
1031	SRscmax = zero
1032	SRsc(n,:) = zero
1033	SRmin = 1.0
1034
1035	ind = MAXLOC(SRnorm(n,:),1)
1036	DO i = 1,ind
1037	SRmin = MIN(SRnorm(n,i),SRmin)
1038	ENDDO
1039	DO i = 1,101
1040	SRsc(n,i) = SRnorm(n,i)- SRmin
1041	IF (i.LE.ind)THEN
1042	SRscmax = MAX(SRsc(n,i),SRscmax)
1043	ENDIF
1044	ENDDO
1045	SRsc(n,:)=SRsc(n,:)/SRscmax
1046
1047	moist_round=(NINT(moist_in(n)*100.))/100.
1048	indmc=INT(moist_round/0.01) + un
1049
1050	moistfunc_result(n) = SRsc(n,indmc)
1051	moistfunc_result(n) = MAX( moistcontSAT_min, MIN( un, moistfunc_result(n) ))
1052	ENDDO
1053
1054	END FUNCTION control_moist_func_moyano
1055
1056	!! ================================================================================================================================
1057	!! FUNCTION : control_moist_func
1058	!!
1059	!>\BRIEF Calculate moisture control for litter and soil C decomposition
1060	!!
1061	!! DESCRIPTION : Calculate moisture control factor applied
1062	!! to litter decomposition and to soil carbon decomposition in
1063	!! stomate_soilcarbon.f90 using the following equation: \n
1064	!! \latexonly
1065	!! \input{control_moist_func1.tex}
1066	!! \endlatexonly
1067	!! \n
1068	!! with M the moisture control factor and soilmoisutre, the soil moisture
1069	!! calculated in sechiba.
1070	!! Then, the function is ranged between Moistcont_min and 1:\n
1071	!! \latexonly
1072	!! \input{control_moist_func2.tex}
1073	!! \endlatexonly
1074	!! \n
1075	!! RECENT CHANGE(S) : None
1076	!!
1077	!! RETURN VALUE : ::moistfunc_result
1078	!!
1079	!! REFERENCE(S) : None
1080	!!
1081	!! FLOWCHART : None
1082	!! \n
1083	!_ ================================================================================================================================
1084
1085	FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)
1086
1087	!! 0. Variable and parameter declaration
1088
1089	!! 0.1 Input variables
1090
1091	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid pixel (unitless)
1092	REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in !! relative humidity (unitless)
1093
1094	!! 0.2 Output variables
1095
1096	REAL(r_std), DIMENSION(npts) :: moistfunc_result !! Moisture control factor (0.25-1, unitless)
1097
1098	!! 0.3 Modified variables
1099
1100	!! 0.4 Local variables
1101
1102	!_ ================================================================================================================================
1103
1104	moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
1105	moistfunc_result(:) = MAX( moistcont_min, MIN( un, moistfunc_result(:) ) )
1106
1107	END FUNCTION control_moist_func
1108
1109	!! ================================================================================================================================
1110	!! FUNCTION : control_temp_func
1111	!!
1112	!>\BRIEF Calculate temperature control for litter and soild C decomposition
1113	!!
1114	!! DESCRIPTION : Calculate temperature control factor applied
1115	!! to litter decomposition and to soil carbon decomposition in
1116	!! stomate_soilcarbon.f90 using the following equation: \n
1117	!! \latexonly
1118	!! \input{control_temp_func1.tex}
1119	!! \endlatexonly
1120	!! \n
1121	!! with T the temperature control factor, temp the temperature in Kelvin of
1122	!! the air (for aboveground litter) or of the soil (for belowground litter
1123	!! and soil)
1124	!! Then, the function is limited in its maximal range to 1:\n
1125	!! \latexonly
1126	!! \input{control_temp_func2.tex}
1127	!! \endlatexonly
1128	!! \n
1129	!! RECENT CHANGE(S) : None
1130	!!
1131	!! RETURN VALUE: ::tempfunc_result
1132	!!
1133	!! REFERENCE(S) : None
1134	!!
1135	!! FLOWCHART : None
1136	!! \n
1137	!_ ================================================================================================================================
1138
1139	FUNCTION control_temp_func (npts, temp_in) RESULT (tempfunc_result)
1140
1141	!! 0. Variable and parameter declaration
1142
1143	!! 0.1 Input variables
1144	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of land pixels (unitless)
1145	REAL(r_std), DIMENSION(npts), INTENT(in) :: temp_in !! Temperature (K)
1146
1147	!! 0.2 Output variables
1148	REAL(r_std), DIMENSION(npts) :: tempfunc_result !! Temperature control factor (0-1, unitless)
1149
1150	!! 0.3 Modified variables
1151
1152	!! 0.4 Local variables
1153
1154	!_ ================================================================================================================================
1155
1156	tempfunc_result(:) = exp( soil_Q10 * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
1157	tempfunc_result(:) = MIN( un, tempfunc_result(:) )
1158
1159	END FUNCTION control_temp_func
1160
1161	END MODULE stomate_litter

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

Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_litter.f90 @ 8579

Download in other formats: