Context Navigation

stomate_soilcarbon.f90 @ 5491

Last change on this file since 5491 was 5323, checked in by chunjing.qiu, 6 years ago
Add oldpeat
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 24.2 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_soilcarbon
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Calculate soil dynamics largely following the Century model
10	!!
11	!!\n DESCRIPTION: None
12	!!
13	!! RECENT CHANGE(S): None
14	!!
15	!! SVN :
16	!! $HeadURL$
17	!! $Date$
18	!! $Revision$
19	!! \n
20	!_ ================================================================================================================================
21
22	MODULE stomate_soilcarbon
23
24	! modules used:
25
26	USE ioipsl_para
27	USE stomate_data
28	USE constantes
29
30	IMPLICIT NONE
31
32	! private & public routines
33
34	PRIVATE
35	PUBLIC soilcarbon,soilcarbon_clear
36
37	! Variables shared by all subroutines in this module
38
39	LOGICAL, SAVE :: firstcall_soilcarbon = .TRUE. !! Is this the first call? (true/false)
40	!$OMP THREADPRIVATE(firstcall_soilcarbon)
41
42	CONTAINS
43
44
45	!! ================================================================================================================================
46	!! SUBROUTINE : soilcarbon_clear
47	!!
48	!>\BRIEF Set the flag ::firstcall_soilcarbon to .TRUE. and as such activate sections 1.1.2 and 1.2 of the subroutine soilcarbon
49	!! (see below).
50	!!
51	!_ ================================================================================================================================
52
53	SUBROUTINE soilcarbon_clear
54	firstcall_soilcarbon=.TRUE.
55	ENDSUBROUTINE soilcarbon_clear
56
57
58	!! ================================================================================================================================
59	!! SUBROUTINE : soilcarbon
60	!!
61	!>\BRIEF Computes the soil respiration and carbon stocks, essentially
62	!! following Parton et al. (1987).
63	!!
64	!! DESCRIPTION : The soil is divided into 3 carbon pools, with different
65	!! characteristic turnover times : active (1-5 years), slow (20-40 years)
66	!! and passive (200-1500 years).\n
67	!! There are three types of carbon transferred in the soil:\n
68	!! - carbon input in active and slow pools from litter decomposition,\n
69	!! - carbon fluxes between the three pools,\n
70	!! - carbon losses from the pools to the atmosphere, i.e., soil respiration.\n
71	!!
72	!! The subroutine performs the following tasks:\n
73	!!
74	!! Section 1.\n
75	!! The flux fractions (f) between carbon pools are defined based on Parton et
76	!! al. (1987). The fractions are constants, except for the flux fraction from
77	!! the active pool to the slow pool, which depends on the clay content,\n
78	!! \latexonly
79	!! \input{soilcarbon_eq1.tex}
80	!! \endlatexonly\n
81	!! In addition, to each pool is assigned a constant turnover time.\n
82	!!
83	!! Section 2.\n
84	!! The carbon input, calculated in the stomate_litter module, is added to the
85	!! carbon stock of the different pools.\n
86	!!
87	!! Section 3.\n
88	!! First, the outgoing carbon flux of each pool is calculated. It is
89	!! proportional to the product of the carbon stock and the ratio between the
90	!! iteration time step and the residence time:\n
91	!! \latexonly
92	!! \input{soilcarbon_eq2.tex}
93	!! \endlatexonly
94	!! ,\n
95	!! Note that in the case of crops, the additional multiplicative factor
96	!! integrates the faster decomposition due to tillage (following Gervois et
97	!! al. (2008)).
98	!! In addition, the flux from the active pool depends on the clay content:\n
99	!! \latexonly
100	!! \input{soilcarbon_eq3.tex}
101	!! \endlatexonly
102	!! ,\n
103	!! Each pool is then cut from the carbon amount corresponding to each outgoing
104	!! flux:\n
105	!! \latexonly
106	!! \input{soilcarbon_eq4.tex}
107	!! \endlatexonly\n
108	!! Second, the flux fractions lost to the atmosphere is calculated in each pool
109	!! by subtracting from 1 the pool-to-pool flux fractions. The soil respiration
110	!! is then the summed contribution of all the pools,\n
111	!! \latexonly
112	!! \input{soilcarbon_eq5.tex}
113	!! \endlatexonly\n
114	!! Finally, each carbon pool accumulates the contribution of the other pools:
115	!! \latexonly
116	!! \input{soilcarbon_eq6.tex}
117	!! \endlatexonly
118	!!
119	!! Section 4.\n
120	!! If the flag SPINUP_ANALYTIC is set to true, the matrix A is updated following
121	!! Lardy (2011).
122	!!
123	!! RECENT CHANGE(S): None
124	!!
125	!! MAIN OUTPUTS VARIABLE(S): carbon, resp_hetero_soil
126	!!
127	!! REFERENCE(S) :
128	!! - Parton, W.J., D.S. Schimel, C.V. Cole, and D.S. Ojima. 1987. Analysis of
129	!! factors controlling soil organic matter levels in Great Plains grasslands.
130	!! Soil Sci. Soc. Am. J., 51, 1173-1179.
131	!! - Gervois, S., P. Ciais, N. de Noblet-Ducoudre, N. Brisson, N. Vuichard,
132	!! and N. Viovy (2008), Carbon and water balance of European croplands
133	!! throughout the 20th century, Global Biogeochem. Cycles, 22, GB2022,
134	!! doi:10.1029/2007GB003018.
135	!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
136	!! Environmental Modelling & Software (2011), doi:10.1016\|j.envsoft.2011.05.016
137	!!
138	!! FLOWCHART :
139	!! \latexonly
140	!! \includegraphics[scale=0.5]{soilcarbon_flowchart.jpg}
141	!! \endlatexonly
142	!! \n
143	!_ ================================================================================================================================
144
145	SUBROUTINE soilcarbon (npts, dt, clay, &
146	soilcarbon_input, control_temp, control_moist, &
147	carbon, &
148	resp_hetero_soil, &
149	MatrixA, &
150	!!!qcj++ peatland
151	height_acro,height_cato,carbon_acro,carbon_cato,tcarbon_acro,tcarbon_cato,resp_acro_oxic, &
152	resp_acro_anoxic,resp_cato,acro_to_cato,litter_to_acro, wtp_peat)
153
154	!gmjc
155	! resp_hetero_soil_part)
156	!end gmjc
157
158	!! 0. Variable and parameter declaration
159
160	!! 0.1 Input variables
161
162	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
163	REAL(r_std), INTENT(in) :: dt !! Time step
164	REAL(r_std), DIMENSION(npts), INTENT(in) :: clay !! Clay fraction (unitless, 0-1)
165	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in) :: soilcarbon_input !! Amount of carbon going into the carbon pools from litter decomposition \f$(gC m^{-2} day^{-1})$\f
166	REAL(r_std), DIMENSION(npts,nlevs), INTENT(in) :: control_temp !! Temperature control of heterotrophic respiration (unitless: 0->1)
167	REAL(r_std), DIMENSION(npts,nlevs), INTENT(in) :: control_moist !! Moisture control of heterotrophic respiration (unitless: 0.25->1)
168
169	!! 0.2 Output variables
170
171	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_hetero_soil !! Soil heterotrophic respiration \f$(gC m^{-2} dt^{-1})$\f
172
173	!! 0.3 Modified variables
174
175	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f
176	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA !! Matrix containing the fluxes between the carbon pools
177	!! per sechiba time step
178	!! @tex $(gC.m^2.day^{-1})$ @endtex
179	!gmjc
180	!! 0.4 Local variables
181	REAL(r_std), SAVE, DIMENSION(ncarb) :: carbon_tau !! Residence time in carbon pools (days)
182	!$OMP THREADPRIVATE(carbon_tau)
183	REAL(r_std), DIMENSION(npts,ncarb,ncarb) :: frac_carb !! Flux fractions between carbon pools
184	!! (second index=origin, third index=destination)
185	!! (unitless, 0-1)
186	REAL(r_std), DIMENSION(npts,ncarb) :: frac_resp !! Flux fractions from carbon pools to the atmosphere (respiration) (unitless, 0-1)
187	REAL(r_std), DIMENSION(npts,ncarb,nelements) :: fluxtot !! Total flux out of carbon pools \f$(gC m^{2})$\f
188	REAL(r_std), DIMENSION(npts,ncarb,ncarb,nelements) :: flux !! Fluxes between carbon pools \f$(gC m^{2})$\f
189	CHARACTER(LEN=7), DIMENSION(ncarb) :: carbon_str !! Name of the carbon pools for informative outputs (unitless)
190	INTEGER(i_std) :: k,kk,m,j !! Indices (unitless)
191
192	!!!qcj++ peatland
193
194	REAL(r_std),DIMENSION(npts),INTENT(in) :: wtp_peat
195	!! 0.2 Output variables
196	REAL(r_std), DIMENSION(npts),INTENT(out) :: tcarbon_acro
197	REAL(r_std), DIMENSION(npts),INTENT(out) :: tcarbon_cato
198	REAL(r_std), DIMENSION(npts,nvm),INTENT(out) :: resp_acro_oxic !!respiration of acrotelm( oxic )
199	REAL(r_std), DIMENSION(npts,nvm),INTENT(out) :: resp_acro_anoxic !!respiration of acrotelm( anoxic )
200	REAL(r_std), DIMENSION(npts,nvm),INTENT(out) :: resp_cato !!respiration of catotelm
201	REAL(r_std), DIMENSION(npts,nvm),INTENT(out) :: litter_to_acro !!total carbon input from litter to acrotelm
202	REAL(r_std), DIMENSION(npts,nvm),INTENT(out) :: acro_to_cato !!carbon flux from acrotelm to catotelm
203	REAL(r_std), DIMENSION(npts), INTENT(out) :: height_cato
204
205	!! 0.3 Modified variables
206	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) ::carbon_acro
207	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) ::carbon_cato
208	REAL(r_std), DIMENSION(npts), INTENT(inout) ::height_acro
209
210	!! 0.4 Local variables
211	REAL(r_std), DIMENSION(npts) ::B !acrotelm oxic respiration
212	REAL(r_std), DIMENSION(npts) ::KA !acrotelm decomposition rate
213	REAL(r_std), DIMENSION(npts) ::KP !catotelm formation rate
214	REAL(r_std), DIMENSION(npts) ::KC !catotelm decomposition rate
215	INTEGER(i_std) ::pft_peat
216	REAL(r_std) ::wtd_min !interface between acrotelm and catotelm
217	REAL(r_std), DIMENSION(npts) ::wtd
218
219
220	!_ ================================================================================================================================
221
222	!! printlev is the level of diagnostic information, 0 (none) to 4 (full)
223	IF (printlev>=3) WRITE(numout,*) 'Entering soilcarbon'
224
225	!! 1. Initializations
226
227	!! 1.1 Get soil "constants"
228	!! 1.1.1 Flux fractions between carbon pools: depend on clay content, recalculated each time
229	! From active pool: depends on clay content
230	frac_carb(:,iactive,iactive) = zero
231	frac_carb(:,iactive,ipassive) = frac_carb_ap
232	frac_carb(:,iactive,islow) = un - (metabolic_ref_frac - active_to_pass_clay_frac*clay(:)) - frac_carb(:,iactive,ipassive)
233
234	! 1.1.1.2 from slow pool
235
236	frac_carb(:,islow,islow) = zero
237	frac_carb(:,islow,iactive) = frac_carb_sa
238	frac_carb(:,islow,ipassive) = frac_carb_sp
239
240	! From passive pool
241	frac_carb(:,ipassive,ipassive) = zero
242	frac_carb(:,ipassive,iactive) = frac_carb_pa
243	frac_carb(:,ipassive,islow) = frac_carb_ps
244
245	IF ( firstcall_soilcarbon ) THEN
246
247	!!!qcj++ peatland
248	! ok_peat = .FALSE.
249	! CALL getin_p('OK_PEAT', ok_peat)
250
251	!! 1.1.2 Residence times in carbon pools (days)
252	carbon_tau(iactive) = carbon_tau_iactive * one_year ! 1.5 years. This is same as CENTURY. But, in Parton et al. (1987), it's weighted by moisture and temperature dependences.
253	carbon_tau(islow) = carbon_tau_islow * one_year ! 25 years. This is same as CENTURY. But, in Parton et al. (1987), it's weighted by moisture and temperature dependences.
254	carbon_tau(ipassive) = carbon_tau_ipassive * one_year ! 1000 years. This is same as CENTURY. But, in Parton et al. (1987), it's weighted by moisture and temperature dependences.
255
256	!! 1.2 Messages : display the residence times
257	carbon_str(iactive) = 'active'
258	carbon_str(islow) = 'slow'
259	carbon_str(ipassive) = 'passive'
260
261	WRITE(numout,*) 'soilcarbon:'
262
263	WRITE(numout,*) ' > minimal carbon residence time in carbon pools (d):'
264	DO k = 1, ncarb ! Loop over carbon pools
265	WRITE(numout,*) '(1, ::carbon_str(k)):',carbon_str(k),' : (1, ::carbon_tau(k)):',carbon_tau(k)
266	ENDDO
267
268	WRITE(numout,*) ' > flux fractions between carbon pools: depend on clay content'
269
270	firstcall_soilcarbon = .FALSE.
271
272	ENDIF
273
274	!! 1.3 Set soil respiration to zero
275	resp_hetero_soil(:,:) = zero
276	! gmjc
277	! resp_hetero_soil_part(:,:,:) = zero
278	! end gmjc
279	!! 2. Update the carbon stocks with the different soil carbon input
280
281	carbon(:,:,:) = carbon(:,:,:) + soilcarbon_input(:,:,:) * dt
282
283	!! 3. Fluxes between carbon reservoirs, and to the atmosphere (respiration) \n
284
285	!! 3.1. Determine the respiration fraction : what's left after
286	! subtracting all the 'pool-to-pool' flux fractions
287	! Diagonal elements of frac_carb are zero
288	! VPP killer:
289	! frac_resp(:,:) = 1. - SUM( frac_carb(:,:,:), DIM=3 )
290	frac_resp(:,:) = un - frac_carb(:,:,iactive) - frac_carb(:,:,islow) - &
291	frac_carb(:,:,ipassive)
292	!! 3.2. Calculate fluxes
293	DO m = 2,nvm ! Loop over # PFTs
294	!! 3.2.1. Flux out of pools
295	!!!qcj++ peatland
296	IF (.NOT. is_peat(m)) THEN
297	DO k = 1, ncarb ! Loop over carbon pools from which the flux comes
298
299	! Determine total flux out of pool
300	! S.L. Piao 2006/05/05 - for crop multiply tillage factor of decomposition
301	! Not crop
302
303	IF ( natural(m) ) THEN
304	fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * &
305	control_moist(:,ibelow) * control_temp(:,ibelow)
306	! C3 crop
307	ELSEIF ( (.NOT. natural(m)) .AND. (.NOT. is_c4(m)) ) THEN
308	fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * &
309	control_moist(:,ibelow) * control_temp(:,ibelow) * flux_tot_coeff(1)
310	! C4 Crop
311	ELSEIF ( (.NOT. natural(m)) .AND. is_c4(m) ) THEN
312	fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * &
313	control_moist(:,ibelow) * control_temp(:,ibelow) * flux_tot_coeff(2)
314	ENDIF
315	! END - S.L. Piao 2006/05/05 - for crop multiply tillage factor of decomposition
316
317	! Carbon flux from active pools depends on clay content
318	IF ( k .EQ. iactive ) THEN
319	fluxtot(:,k,icarbon) = fluxtot(:,k,icarbon) * ( un - flux_tot_coeff(3) * clay(:) )
320	ENDIF
321
322	! Update the loss in each carbon pool
323	carbon(:,k,m) = carbon(:,k,m) - fluxtot(:,k,icarbon)
324
325	! Fluxes towards the other pools (k -> kk)
326	DO kk = 1, ncarb ! Loop over the carbon pools where the flux goes
327	flux(:,k,kk,icarbon) = frac_carb(:,k,kk) * fluxtot(:,k,icarbon)
328	ENDDO
329
330	ENDDO ! End of loop over carbon pools
331
332	!! 3.2.2 respiration
333	! VPP killer:
334	! resp_hetero_soil(:,m) = SUM( frac_resp(:,:) * fluxtot(:,:), DIM=2 ) / dt
335
336	resp_hetero_soil(:,m) = &
337	( frac_resp(:,iactive) * fluxtot(:,iactive,icarbon) + &
338	frac_resp(:,islow) * fluxtot(:,islow,icarbon) + &
339	frac_resp(:,ipassive) * fluxtot(:,ipassive,icarbon) ) / dt
340
341	!gmjc
342	! resp_hetero_soil_part(:,iactive,m) = &
343	! frac_resp(:,iactive) * fluxtot(:,iactive,icarbon)/dt
344	! resp_hetero_soil_part(:,islow,m) = &
345	! frac_resp(:,islow) * fluxtot(:,islow,icarbon)/dt
346	! resp_hetero_soil_part(:,ipassive,m) = &
347	! frac_resp(:,ipassive) * fluxtot(:,ipassive,icarbon)/dt
348	!end gmjc
349	!! 3.2.3 add fluxes to active, slow, and passive pools
350	! VPP killer:
351	! carbon(:,:,m) = carbon(:,:,m) + SUM( flux(:,:,:), DIM=2 )
352
353	DO k = 1, ncarb ! Loop over carbon pools
354	carbon(:,k,m) = carbon(:,k,m) + &
355	flux(:,iactive,k,icarbon) + flux(:,ipassive,k,icarbon) + flux(:,islow,k,icarbon)
356	ENDDO ! Loop over carbon pools
357
358	IF (ok_peat) THEN
359	carbon_acro(:,m)=zero
360	carbon_cato(:,m)=zero
361	acro_to_cato(:,m)=zero
362	litter_to_acro(:,m)=zero
363	resp_cato(:,m)=zero
364	resp_acro_anoxic(:,m)=zero
365	resp_acro_oxic(:,m)=zero
366	ENDIF
367	ENDIF ! is_peat
368
369	ENDDO ! End loop over PFTs
370
371	!! 4. (Quasi-)Analytical Spin-up
372
373	!! 4.1.1 Finish to fill MatrixA with fluxes between soil pools
374
375	IF (spinup_analytic) THEN
376
377	DO m = 2,nvm
378
379	! flux leaving the active pool
380	MatrixA(:,m,iactive_pool,iactive_pool) = moins_un * &
381	dt/carbon_tau(iactive) * &
382	control_moist(:,ibelow) * control_temp(:,ibelow) * &
383	( 1. - flux_tot_coeff(3) * clay(:))
384
385	! flux received by the active pool from the slow pool
386	MatrixA(:,m,iactive_pool,islow_pool) = frac_carb(:,islow,iactive)dt/carbon_tau(islow) &
387	control_moist(:,ibelow) * control_temp(:,ibelow)
388
389	! flux received by the active pool from the passive pool
390	MatrixA(:,m,iactive_pool,ipassive_pool) = frac_carb(:,ipassive,iactive)dt/carbon_tau(ipassive) &
391	control_moist(:,ibelow) * control_temp(:,ibelow)
392
393	! flux received by the slow pool from the active pool
394	MatrixA(:,m,islow_pool,iactive_pool) = frac_carb(:,iactive,islow) *&
395	dt/carbon_tau(iactive) * &
396	control_moist(:,ibelow) * control_temp(:,ibelow) * &
397	( 1. - flux_tot_coeff(3) * clay(:) )
398
399	! flux leaving the slow pool
400	MatrixA(:,m,islow_pool,islow_pool) = moins_un * &
401	dt/carbon_tau(islow) * &
402	control_moist(:,ibelow) * control_temp(:,ibelow)
403
404	! flux received by the passive pool from the active pool
405	MatrixA(:,m,ipassive_pool,iactive_pool) = frac_carb(:,iactive,ipassive)* &
406	dt/carbon_tau(iactive) * &
407	control_moist(:,ibelow) * control_temp(:,ibelow) *&
408	( 1. - flux_tot_coeff(3) * clay(:) )
409
410	! flux received by the passive pool from the slow pool
411	MatrixA(:,m,ipassive_pool,islow_pool) = frac_carb(:,islow,ipassive) * &
412	dt/carbon_tau(islow) * &
413	control_moist(:,ibelow) * control_temp(:,ibelow)
414
415	! flux leaving the passive pool
416	MatrixA(:,m,ipassive_pool,ipassive_pool) = moins_un * &
417	dt/carbon_tau(ipassive) * &
418	control_moist(:,ibelow) * control_temp(:,ibelow)
419
420
421	IF ( (.NOT. natural(m)) .AND. (.NOT. is_c4(m)) ) THEN ! C3crop
422
423	! flux leaving the active pool
424	MatrixA(:,m,iactive_pool,iactive_pool) = MatrixA(:,m,iactive_pool,iactive_pool) * &
425	flux_tot_coeff(1)
426
427	! flux received by the active pool from the slow pool
428	MatrixA(:,m,iactive_pool,islow_pool)= MatrixA(:,m,iactive_pool,islow_pool) * &
429	flux_tot_coeff(1)
430
431	! flux received by the active pool from the passive pool
432	MatrixA(:,m,iactive_pool,ipassive_pool) = MatrixA(:,m,iactive_pool,ipassive_pool) * &
433	flux_tot_coeff(1)
434
435	! flux received by the slow pool from the active pool
436	MatrixA(:,m,islow_pool,iactive_pool) = MatrixA(:,m,islow_pool,iactive_pool) * &
437	flux_tot_coeff(1)
438
439	! flux leaving the slow pool
440	MatrixA(:,m,islow_pool,islow_pool) = MatrixA(:,m,islow_pool,islow_pool) * &
441	flux_tot_coeff(1)
442
443	! flux received by the passive pool from the active pool
444	MatrixA(:,m,ipassive_pool,iactive_pool) = MatrixA(:,m,ipassive_pool,iactive_pool) * &
445	flux_tot_coeff(1)
446
447	! flux received by the passive pool from the slow pool
448	MatrixA(:,m,ipassive_pool,islow_pool) = MatrixA(:,m,ipassive_pool,islow_pool) * &
449	flux_tot_coeff(1)
450
451	! flux leaving the passive pool
452	MatrixA(:,m,ipassive_pool,ipassive_pool) = MatrixA(:,m,ipassive_pool,ipassive_pool) *&
453	flux_tot_coeff(1)
454
455	ENDIF ! (.NOT. natural(m)) .AND. (.NOT. is_c4(m))
456
457
458	IF ( (.NOT. natural(m)) .AND. is_c4(m) ) THEN ! C4crop
459
460	! flux leaving the active pool
461	MatrixA(:,m,iactive_pool,iactive_pool) = MatrixA(:,m,iactive_pool,iactive_pool) * &
462	flux_tot_coeff(2)
463
464	! flux received by the active pool from the slow pool
465	MatrixA(:,m,iactive_pool,islow_pool)= MatrixA(:,m,iactive_pool,islow_pool) * &
466	flux_tot_coeff(2)
467
468	! flux received by the active pool from the passive pool
469	MatrixA(:,m,iactive_pool,ipassive_pool) = MatrixA(:,m,iactive_pool,ipassive_pool) * &
470	flux_tot_coeff(2)
471
472	! flux received by the slow pool from the active pool
473	MatrixA(:,m,islow_pool,iactive_pool) = MatrixA(:,m,islow_pool,iactive_pool) * &
474	flux_tot_coeff(2)
475
476	! flux leaving the slow pool
477	MatrixA(:,m,islow_pool,islow_pool) = MatrixA(:,m,islow_pool,islow_pool) * &
478	flux_tot_coeff(2)
479
480	! flux received by the passive pool from the active pool
481	MatrixA(:,m,ipassive_pool,iactive_pool) = MatrixA(:,m,ipassive_pool,iactive_pool) * &
482	flux_tot_coeff(2)
483
484	! flux received by the passive pool from the slow pool
485	MatrixA(:,m,ipassive_pool,islow_pool) = MatrixA(:,m,ipassive_pool,islow_pool) * &
486	flux_tot_coeff(2)
487
488	! flux leaving the passive pool
489	MatrixA(:,m,ipassive_pool,ipassive_pool) = MatrixA(:,m,ipassive_pool,ipassive_pool) * &
490	flux_tot_coeff(2)
491
492	ENDIF ! (.NOT. natural(m)) .AND. is_c4(m)
493
494	IF (printlev>=4) WRITE(numout,*)'Finish to fill MatrixA'
495
496	ENDDO ! Loop over # PFTS
497
498
499	! 4.2 Add Identity for each submatrix(7,7)
500
501	DO j = 1,nbpools
502	MatrixA(:,:,j,j) = MatrixA(:,:,j,j) + un
503	ENDDO
504
505	ENDIF ! (spinup_analytic)
506
507	IF (printlev>=4) WRITE(numout,*) 'Leaving soilcarbon'
508
509	IF (ok_peat) THEN
510	!!! decomposition rates and catotelm formation rates are constrained by temperature
511	KA(:)=control_temp(:,ibelow)KA_inidt/one_year !control_temp(:,ibelow): integrated temperature of of 11layers
512	KP(:)=control_temp(:,ibelow)KP_inidt/one_year
513	KC(:)=KC_inidt/one_year !control_temp(:,ibelow)?
514
515	height_cato(:)=zero
516	tcarbon_acro(:) = zero
517	tcarbon_cato(:) = zero
518	resp_acro_oxic(:,:)=zero
519	resp_acro_anoxic(:,:)=zero
520	resp_cato(:,:)=zero
521	litter_to_acro(:,:)=zero
522	acro_to_cato(:,:)=zero
523	wtd_min=300. !in mm
524	CALL getin_p('wtd_min', wtd_min)
525
526	DO k=1,npts
527	wtd(k)=wtd_min-wtp_peat(k) !in mm
528	wtd(k)=wtd(k)*0.001 ! in m
529	IF (wtd(k) .LE. 0.0) THEN
530	B(k)=1.0
531	ELSE IF (wtd(k) .LT. height_acro(k)) THEN
532	B(k)=(height_acro(k)-wtd(k))/height_acro(k)
533	ELSE IF (wtd(k) .GE. height_acro(k)) THEN
534	B(k)=0.0
535	ENDIF
536	ENDDO
537
538	DO k=1,npts
539	DO m=2,nvm
540	IF (is_peat(m)) THEN
541
542	carbon(k,:,m) =zero
543
544	resp_acro_oxic(k,m)=resp_acro_oxic(k,m)+B(k)KA(k)carbon_acro(k,m)
545	resp_acro_anoxic(k,m)=resp_acro_anoxic(k,m)+(1.-B(k))v_ratioKA(k)*carbon_acro(k,m)
546	acro_to_cato(k,m)=acro_to_cato(k,m)+KP(k)*carbon_acro(k,m)
547	resp_cato(k,m)=resp_cato(k,m)+KC(k)*carbon_cato(k,m)
548	!!!litter added to default active and slow soil carbon pool are added into acrotelm
549	resp_hetero_soil(k,m) = (resp_acro_oxic(k,m)+resp_acro_anoxic(k,m)+resp_cato(k,m))/dt
550	litter_to_acro(k,m)=litter_to_acro(k,m)+soilcarbon_input(k,iactive,m)*dt+ &
551	& soilcarbon_input(k,islow,m)*dt
552	carbon_acro(k,m)=carbon_acro(k,m)+litter_to_acro(k,m)- &
553	& acro_to_cato(k,m)-resp_acro_oxic(k,m)-resp_acro_anoxic(k,m)
554	carbon_cato(k,m)=carbon_cato(k,m)+acro_to_cato(k,m)-resp_cato(k,m)
555
556	tcarbon_acro(k)=tcarbon_acro(k)+carbon_acro(k,m)
557	tcarbon_cato(k)=tcarbon_cato(k)+carbon_cato(k,m)
558	height_acro(k)=tcarbon_acro(k)/(p_A*cf_A)
559	height_cato(k)=tcarbon_cato(k)/(p_C*cf_C)
560	ENDIF
561	ENDDO
562	ENDDO
563	ENDIF
564
565	END SUBROUTINE soilcarbon
566
567	END MODULE stomate_soilcarbon

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

Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_stomate/stomate_soilcarbon.f90 @ 5491

Download in other formats: