source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/stomate_soilcarbon.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 19.6 KB
Line 
1! =================================================================================================================================
2! MODULE       : stomate_soilcarbon
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       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
22MODULE 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
42CONTAINS
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!gmjc
151!       resp_hetero_soil_part)
152!end gmjc
153
154!! 0. Variable and parameter declaration
155
156    !! 0.1 Input variables
157   
158    INTEGER(i_std), INTENT(in)                            :: npts             !! Domain size (unitless)
159    REAL(r_std), INTENT(in)                               :: dt               !! Time step
160    REAL(r_std), DIMENSION(npts), INTENT(in)              :: clay             !! Clay fraction (unitless, 0-1)
161    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
162    REAL(r_std), DIMENSION(npts,nlevs), INTENT(in)        :: control_temp     !! Temperature control of heterotrophic respiration (unitless: 0->1)
163    REAL(r_std), DIMENSION(npts,nlevs), INTENT(in)        :: control_moist    !! Moisture control of heterotrophic respiration (unitless: 0.25->1)
164
165    !! 0.2 Output variables
166   
167    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)         :: resp_hetero_soil !! Soil heterotrophic respiration \f$(gC m^{-2} dt^{-1})$\f
168
169    !! 0.3 Modified variables
170   
171    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon             !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f
172    REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA  !! Matrix containing the fluxes between the carbon pools
173                                                                                !! per sechiba time step
174                                                                                !! @tex $(gC.m^2.day^{-1})$ @endtex
175!gmjc
176    !! 0.4 Local variables
177    REAL(r_std), SAVE, DIMENSION(ncarb)                   :: carbon_tau       !! Residence time in carbon pools (days)
178!$OMP THREADPRIVATE(carbon_tau)
179    REAL(r_std), DIMENSION(npts,ncarb,ncarb)              :: frac_carb        !! Flux fractions between carbon pools
180                                                                              !! (second index=origin, third index=destination)
181                                                                              !! (unitless, 0-1)
182    REAL(r_std), DIMENSION(npts,ncarb)                    :: frac_resp        !! Flux fractions from carbon pools to the atmosphere (respiration) (unitless, 0-1)
183    REAL(r_std), DIMENSION(npts,ncarb,nelements)          :: fluxtot          !! Total flux out of carbon pools \f$(gC m^{2})$\f
184    REAL(r_std), DIMENSION(npts,ncarb,ncarb,nelements)    :: flux             !! Fluxes between carbon pools \f$(gC m^{2})$\f
185    CHARACTER(LEN=7), DIMENSION(ncarb)                    :: carbon_str       !! Name of the carbon pools for informative outputs (unitless)
186    INTEGER(i_std)                                        :: k,kk,m,j         !! Indices (unitless)
187
188!_ ================================================================================================================================
189
190    !! printlev is the level of diagnostic information, 0 (none) to 4 (full)
191    IF (printlev>=3) WRITE(numout,*) 'Entering soilcarbon' 
192
193!! 1. Initializations
194
195    !! 1.1 Get soil "constants"
196    !! 1.1.1 Flux fractions between carbon pools: depend on clay content, recalculated each time
197    ! From active pool: depends on clay content
198    frac_carb(:,iactive,iactive) = zero
199    frac_carb(:,iactive,ipassive) = frac_carb_ap
200    frac_carb(:,iactive,islow) = un - (metabolic_ref_frac - active_to_pass_clay_frac*clay(:)) - frac_carb(:,iactive,ipassive)
201
202    ! 1.1.1.2 from slow pool
203
204    frac_carb(:,islow,islow) = zero
205    frac_carb(:,islow,iactive) = frac_carb_sa
206    frac_carb(:,islow,ipassive) = frac_carb_sp
207
208    ! From passive pool
209    frac_carb(:,ipassive,ipassive) = zero
210    frac_carb(:,ipassive,iactive) = frac_carb_pa
211    frac_carb(:,ipassive,islow) = frac_carb_ps
212
213    IF ( firstcall_soilcarbon ) THEN
214
215        !! 1.1.2 Residence times in carbon pools (days)
216        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.
217        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.
218        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.
219       
220        !! 1.2 Messages : display the residence times 
221        carbon_str(iactive) = 'active'
222        carbon_str(islow) = 'slow'
223        carbon_str(ipassive) = 'passive'
224       
225        IF (printlev >= 2) THEN
226           WRITE(numout,*) 'soilcarbon:'
227           WRITE(numout,*) '   > minimal carbon residence time in carbon pools (d):'
228           DO k = 1, ncarb ! Loop over carbon pools
229              WRITE(numout,*) '(1, ::carbon_str(k)):',carbon_str(k),' : (1, ::carbon_tau(k)):',carbon_tau(k)
230           ENDDO
231           WRITE(numout,*) '   > flux fractions between carbon pools: depend on clay content'
232        END IF
233        firstcall_soilcarbon = .FALSE.
234       
235    ENDIF
236
237    !! 1.3 Set soil respiration to zero
238    resp_hetero_soil(:,:) = zero
239! gmjc
240!    resp_hetero_soil_part(:,:,:) = zero
241! end gmjc
242!! 2. Update the carbon stocks with the different soil carbon input
243
244    carbon(:,:,:) = carbon(:,:,:) + soilcarbon_input(:,:,:) * dt
245
246!! 3. Fluxes between carbon reservoirs, and to the atmosphere (respiration) \n
247
248    !! 3.1. Determine the respiration fraction : what's left after
249    ! subtracting all the 'pool-to-pool' flux fractions
250    ! Diagonal elements of frac_carb are zero
251    !    VPP killer:
252    !     frac_resp(:,:) = 1. - SUM( frac_carb(:,:,:), DIM=3 )
253    frac_resp(:,:) = un - frac_carb(:,:,iactive) - frac_carb(:,:,islow) - &
254         frac_carb(:,:,ipassive) 
255
256    !! 3.2. Calculate fluxes
257
258    DO m = 2,nvm ! Loop over # PFTs
259
260      !! 3.2.1. Flux out of pools
261
262      DO k = 1, ncarb ! Loop over carbon pools from which the flux comes
263       
264        ! Determine total flux out of pool
265        ! S.L. Piao 2006/05/05 - for crop multiply tillage factor of decomposition
266        ! Not crop
267         IF ( natural(m) ) THEN
268            fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * &
269                  control_moist(:,ibelow) * control_temp(:,ibelow)
270         ! C3 crop
271          ELSEIF ( (.NOT. natural(m)) .AND. (.NOT. is_c4(m)) ) THEN
272             fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * &
273                  control_moist(:,ibelow) * control_temp(:,ibelow) * flux_tot_coeff(1)
274          ! C4 Crop   
275          ELSEIF ( (.NOT. natural(m)) .AND. is_c4(m) ) THEN
276             fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * &
277                  control_moist(:,ibelow) * control_temp(:,ibelow) * flux_tot_coeff(2)
278        ENDIF
279        ! END - S.L. Piao 2006/05/05 - for crop multiply tillage factor of decomposition
280
281        ! Carbon flux from active pools depends on clay content
282        IF ( k .EQ. iactive ) THEN
283            fluxtot(:,k,icarbon) = fluxtot(:,k,icarbon) * ( un - flux_tot_coeff(3) * clay(:) )
284        ENDIF
285       
286        ! Update the loss in each carbon pool
287        carbon(:,k,m) = carbon(:,k,m) - fluxtot(:,k,icarbon)
288       
289        ! Fluxes towards the other pools (k -> kk)
290        DO kk = 1, ncarb ! Loop over the carbon pools where the flux goes
291          flux(:,k,kk,icarbon) = frac_carb(:,k,kk) * fluxtot(:,k,icarbon)
292        ENDDO
293       
294      ENDDO ! End of loop over carbon pools
295     
296      !! 3.2.2 respiration
297      !       VPP killer:
298      !       resp_hetero_soil(:,m) = SUM( frac_resp(:,:) * fluxtot(:,:), DIM=2 ) / dt
299     
300      resp_hetero_soil(:,m) = &
301         ( frac_resp(:,iactive) * fluxtot(:,iactive,icarbon) + &
302         frac_resp(:,islow) * fluxtot(:,islow,icarbon) + &
303         frac_resp(:,ipassive) * fluxtot(:,ipassive,icarbon)  ) / dt
304!gmjc
305!       resp_hetero_soil_part(:,iactive,m) = &
306!            frac_resp(:,iactive) * fluxtot(:,iactive,icarbon)/dt
307!       resp_hetero_soil_part(:,islow,m) = &
308!            frac_resp(:,islow) * fluxtot(:,islow,icarbon)/dt
309!       resp_hetero_soil_part(:,ipassive,m) = &
310!            frac_resp(:,ipassive) * fluxtot(:,ipassive,icarbon)/dt
311!end gmjc     
312      !! 3.2.3 add fluxes to active, slow, and passive pools
313      !       VPP killer:
314      !       carbon(:,:,m) = carbon(:,:,m) + SUM( flux(:,:,:), DIM=2 )
315     
316      DO k = 1, ncarb ! Loop over carbon pools
317        carbon(:,k,m) = carbon(:,k,m) + &
318           flux(:,iactive,k,icarbon) + flux(:,ipassive,k,icarbon) + flux(:,islow,k,icarbon)
319      ENDDO ! Loop over carbon pools
320     
321    ENDDO ! End loop over PFTs
322   
323 !! 4. (Quasi-)Analytical Spin-up
324   
325    !! 4.1.1 Finish to fill MatrixA with fluxes between soil pools
326   
327    IF (spinup_analytic) THEN
328
329       DO m = 2,nvm 
330
331          ! flux leaving the active pool
332          MatrixA(:,m,iactive_pool,iactive_pool) = moins_un * &
333               dt/carbon_tau(iactive) * &
334               control_moist(:,ibelow) * control_temp(:,ibelow) * &
335               ( 1. - flux_tot_coeff(3) * clay(:)) 
336
337          ! flux received by the active pool from the slow pool
338          MatrixA(:,m,iactive_pool,islow_pool) =  frac_carb(:,islow,iactive)*dt/carbon_tau(islow) * &
339               control_moist(:,ibelow) * control_temp(:,ibelow)
340
341          ! flux received by the active pool from the passive pool
342          MatrixA(:,m,iactive_pool,ipassive_pool) =  frac_carb(:,ipassive,iactive)*dt/carbon_tau(ipassive) * &
343               control_moist(:,ibelow) * control_temp(:,ibelow) 
344
345          ! flux received by the slow pool from the active pool
346          MatrixA(:,m,islow_pool,iactive_pool) =  frac_carb(:,iactive,islow) *&
347               dt/carbon_tau(iactive) * &
348               control_moist(:,ibelow) * control_temp(:,ibelow) * &
349               ( 1. - flux_tot_coeff(3) * clay(:) ) 
350
351          ! flux leaving the slow pool
352          MatrixA(:,m,islow_pool,islow_pool) = moins_un * &
353               dt/carbon_tau(islow) * &
354               control_moist(:,ibelow) * control_temp(:,ibelow)
355
356          ! flux received by the passive pool from the active pool
357          MatrixA(:,m,ipassive_pool,iactive_pool) =  frac_carb(:,iactive,ipassive)* &
358               dt/carbon_tau(iactive) * &
359               control_moist(:,ibelow) * control_temp(:,ibelow) *&
360               ( 1. - flux_tot_coeff(3) * clay(:) )
361
362          ! flux received by the passive pool from the slow pool
363          MatrixA(:,m,ipassive_pool,islow_pool) =  frac_carb(:,islow,ipassive) * &
364               dt/carbon_tau(islow) * &
365               control_moist(:,ibelow) * control_temp(:,ibelow)
366
367          ! flux leaving the passive pool
368          MatrixA(:,m,ipassive_pool,ipassive_pool) =  moins_un * &
369               dt/carbon_tau(ipassive) * &
370               control_moist(:,ibelow) * control_temp(:,ibelow)     
371
372
373          IF ( (.NOT. natural(m)) .AND. (.NOT. is_c4(m)) ) THEN ! C3crop
374
375             ! flux leaving the active pool
376             MatrixA(:,m,iactive_pool,iactive_pool) = MatrixA(:,m,iactive_pool,iactive_pool) * &
377                  flux_tot_coeff(1) 
378
379             ! flux received by the active pool from the slow pool
380             MatrixA(:,m,iactive_pool,islow_pool)= MatrixA(:,m,iactive_pool,islow_pool) * &
381                  flux_tot_coeff(1) 
382
383             ! flux received by the active pool from the passive pool
384             MatrixA(:,m,iactive_pool,ipassive_pool) = MatrixA(:,m,iactive_pool,ipassive_pool) * &
385                  flux_tot_coeff(1)   
386
387             ! flux received by the slow pool from the active pool
388             MatrixA(:,m,islow_pool,iactive_pool) =  MatrixA(:,m,islow_pool,iactive_pool) * &
389                  flux_tot_coeff(1) 
390
391             ! flux leaving the slow pool
392             MatrixA(:,m,islow_pool,islow_pool) = MatrixA(:,m,islow_pool,islow_pool) * &
393                  flux_tot_coeff(1)     
394
395             ! flux received by the passive pool from the active pool
396             MatrixA(:,m,ipassive_pool,iactive_pool) = MatrixA(:,m,ipassive_pool,iactive_pool) * &
397                  flux_tot_coeff(1)
398
399             ! flux received by the passive pool from the slow pool
400             MatrixA(:,m,ipassive_pool,islow_pool) = MatrixA(:,m,ipassive_pool,islow_pool) * &
401                  flux_tot_coeff(1)
402
403             ! flux leaving the passive pool
404             MatrixA(:,m,ipassive_pool,ipassive_pool) =  MatrixA(:,m,ipassive_pool,ipassive_pool) *&
405                  flux_tot_coeff(1)
406
407          ENDIF ! (.NOT. natural(m)) .AND. (.NOT. is_c4(m))
408
409
410          IF ( (.NOT. natural(m)) .AND. is_c4(m) ) THEN ! C4crop
411
412             ! flux leaving the active pool
413             MatrixA(:,m,iactive_pool,iactive_pool) = MatrixA(:,m,iactive_pool,iactive_pool) * &
414                  flux_tot_coeff(2) 
415
416            ! flux received by the active pool from the slow pool
417             MatrixA(:,m,iactive_pool,islow_pool)= MatrixA(:,m,iactive_pool,islow_pool) * &
418                  flux_tot_coeff(2) 
419
420             ! flux received by the active pool from the passive pool
421             MatrixA(:,m,iactive_pool,ipassive_pool) = MatrixA(:,m,iactive_pool,ipassive_pool) * &
422                  flux_tot_coeff(2)   
423
424             ! flux received by the slow pool from the active pool
425             MatrixA(:,m,islow_pool,iactive_pool) =  MatrixA(:,m,islow_pool,iactive_pool) * &
426                  flux_tot_coeff(2) 
427
428             ! flux leaving the slow pool
429             MatrixA(:,m,islow_pool,islow_pool) = MatrixA(:,m,islow_pool,islow_pool) * &
430                  flux_tot_coeff(2)     
431
432             ! flux received by the passive pool from the active pool
433             MatrixA(:,m,ipassive_pool,iactive_pool) = MatrixA(:,m,ipassive_pool,iactive_pool) * &
434                  flux_tot_coeff(2)
435
436             ! flux received by the passive pool from the slow pool
437             MatrixA(:,m,ipassive_pool,islow_pool) = MatrixA(:,m,ipassive_pool,islow_pool) * &
438                  flux_tot_coeff(2)
439
440             ! flux leaving the passive pool
441             MatrixA(:,m,ipassive_pool,ipassive_pool) =  MatrixA(:,m,ipassive_pool,ipassive_pool) * &
442                  flux_tot_coeff(2)
443
444          ENDIF ! (.NOT. natural(m)) .AND. is_c4(m)
445
446          IF (printlev>=4) WRITE(numout,*)'Finish to fill MatrixA'
447
448       ENDDO ! Loop over # PFTS
449
450
451       ! 4.2 Add Identity for each submatrix(7,7)
452
453       DO j = 1,nbpools
454          MatrixA(:,:,j,j) = MatrixA(:,:,j,j) + un 
455       ENDDO
456
457    ENDIF ! (spinup_analytic)
458
459    IF (printlev>=4) WRITE(numout,*) 'Leaving soilcarbon'
460   
461  END SUBROUTINE soilcarbon
462
463END MODULE stomate_soilcarbon
Note: See TracBrowser for help on using the repository browser.