source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_stomate/stomate_soilcarbon.f90 @ 8616

Last change on this file since 8616 was 2738, checked in by josefine.ghattas, 9 years ago

Small modifications for TAF :

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