source: branches/publications/ORCHIDEE_SOM_13C_r4859/src_stomate/stomate_lcchange.f90 @ 6889

Last change on this file since 6889 was 4021, checked in by marta.camino, 8 years ago

13C and delta 13C included without 13C discretization in soil

  • Property svn:keywords set to HeadURL Date Author Revision
File size: 21.7 KB
Line 
1! =================================================================================================================================
2! MODULE       : stomate_lcchange
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       Impact of land cover change on carbon stocks
10!!
11!!\n DESCRIPTION: None
12!!
13!! RECENT CHANGE(S): None
14!!
15!! REFERENCE(S) : None
16!!
17!! SVN          :
18!! $HeadURL$
19!! $Date$
20!! $Revision$
21!! \n
22!_ ================================================================================================================================
23
24
25MODULE stomate_lcchange
26
27  ! modules used:
28 
29  USE ioipsl_para
30  USE stomate_data
31  USE pft_parameters
32  USE constantes
33  USE constantes_soil
34 
35  IMPLICIT NONE
36 
37  PRIVATE
38  PUBLIC lcchange_main
39 
40CONTAINS
41
42
43!! ================================================================================================================================
44!! SUBROUTINE   : lcchange_main
45!!
46!>\BRIEF        Impact of land cover change on carbon stocks
47!!
48!! DESCRIPTION  : This subroutine is always activate if VEGET_UPDATE>0Y in the configuration file, which means that the
49!! vegetation map is updated regulary. lcchange_main is called from stomateLpj the first time step after the vegetation
50!! map has been changed.
51!! The impact of land cover change on carbon stocks is computed in this subroutine. The land cover change is written
52!! by the difference of current and previous "maximal" coverage fraction of a PFT.
53!! On the basis of this difference, the amount of 'new establishment'/'biomass export',
54!! and increase/decrease of each component, are estimated.\n
55!!
56!! Main structure of lpj_establish.f90 is:
57!! 1. Initialization
58!! 2. Calculation of changes in carbon stocks and biomass by land cover change
59!! 3. Update 10 year- and 100 year-turnover pool contents
60!! 4. History
61!!
62!! RECENT CHANGE(S) : None
63!!
64!! MAIN OUTPUT VARIABLE(S) : ::prod10, ::prod100, ::flux10, ::flux100,
65!!   :: cflux_prod10 and :: cflux_prod100
66!!
67!! REFERENCES   : None
68!!
69!! FLOWCHART    :
70!! \latexonly
71!!     \includegraphics[scale=0.5]{lcchange.png}
72!! \endlatexonly
73!! \n
74!_ ================================================================================================================================
75
76 
77  SUBROUTINE lcchange_main ( npts, dt_days, veget_max_old, veget_max_new, &
78       biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &       
79       co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,flux10,flux100, &
80       prod10,prod100,&
81       convflux,&
82       cflux_prod10,cflux_prod100, leaf_frac,&
83       npp_longterm, lm_lastyearmax, litter_above, litter_below,&
84       carbon,DOC,moist_soil)
85   
86    IMPLICIT NONE
87   
88  !! 0. Variable and parameter declaration
89   
90    !! 0.1 Input variables
91   
92    INTEGER, INTENT(in)                                       :: npts             !! Domain size - number of pixels (unitless)
93    REAL(r_std), INTENT(in)                                   :: dt_days          !! Time step of vegetation dynamics for stomate
94                                                                                  !! (days)
95    REAL(r_std), DIMENSION(nvm, nparts,nelements), INTENT(in) :: bm_sapl          !! biomass of sapling
96                                                                                  !! @tex ($gC individual^{-1}$) @endtex
97    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)             :: moist_soil       !! moisture in soil for one pixel (m3/m3)
98    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: veget_max_old    !! Current "maximal" coverage fraction of a PFT (LAI
99                                                                                  !! -> infinity) on ground
100    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: veget_max_new    !! New "maximal" coverage fraction of a PFT (LAI ->
101                                                                                  !! infinity) on ground (unitless)
102 
103    !! 0.2 Output variables
104
105    REAL(r_std), DIMENSION(npts), INTENT(out)                 :: convflux         !! release during first year following land cover
106                                                                                  !! change
107    REAL(r_std), DIMENSION(npts), INTENT(out)                 :: cflux_prod10     !! total annual release from the 10 year-turnover
108                                                                                  !! pool @tex ($gC m^{-2}$) @endtex
109    REAL(r_std), DIMENSION(npts), INTENT(out)                 :: cflux_prod100    !! total annual release from the 100 year-
110                                                                                  !! turnover pool @tex ($gC m^{-2}$) @endtex
111    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: turnover_daily   !! Turnover rates
112                                                                                      !! @tex ($gC m^{-2} day^{-1}$) @endtex
113
114    !! 0.3 Modified variables   
115   
116    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: biomass    !! biomass @tex ($gC m^{-2}$) @endtex
117    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: ind              !! Number of individuals @tex ($m^{-2}$) @endtex
118    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: age              !! mean age (years)
119    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)               :: senescence       !! plant senescent (only for deciduous trees) Set
120                                                                                  !! to .FALSE. if PFT is introduced or killed
121    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)               :: PFTpresent       !! Is pft there (unitless)
122    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: everywhere       !! is the PFT everywhere in the grid box or very
123                                                                                  !! localized (unitless)
124    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: when_growthinit  !! how many days ago was the beginning of the
125                                                                                  !! growing season (days)
126    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: co2_to_bm        !! biomass uptaken
127                                                                                  !! @tex ($gC m^{-2} day^{-1}$) @endtex
128    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! conversion of biomass to litter
129                                                                                  !! @tex ($gC m^{-2} day^{-1}$) @endtex
130    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: cn_ind           !! crown area of individuals
131                                                                                  !! @tex ($m^{2}$) @endtex
132    REAL(r_std), DIMENSION(npts,0:10), INTENT(inout)          :: prod10           !! products remaining in the 10 year-turnover
133                                                                                  !! pool after the annual release for each
134                                                                                  !! compartment (10 + 1 : input from year of land
135                                                                                  !! cover change)
136    REAL(r_std), DIMENSION(npts,0:100), INTENT(inout)         :: prod100          !! products remaining in the 100 year-turnover
137                                                                                  !! pool after the annual release for each
138                                                                                  !! compartment (100 + 1 : input from year of land
139                                                                                  !! cover change)
140    REAL(r_std), DIMENSION(npts,10), INTENT(inout)            :: flux10           !! annual release from the 10/100 year-turnover
141                                                                                  !! pool compartments
142    REAL(r_std), DIMENSION(npts,100), INTENT(inout)           :: flux100          !! annual release from the 10/100 year-turnover
143                                                                                  !! pool compartments
144    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac        !! fraction of leaves in leaf age class
145                                                                                  !! (unitless)
146    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: lm_lastyearmax   !! last year's maximum leaf mass for each PFT
147                                                                                  !! @tex ($gC m^{-2}$) @endtex
148    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: npp_longterm     !! "long term" net primary productivity
149                                                                                  !! @tex ($gC m^{-2} year^{-1}$) @endtex
150    REAL(r_std),DIMENSION(npts,nlitt,nvm, nelements), INTENT(inout):: litter_above !! metabolic and structural litter, above
151                                                                                  !! ground @tex ($gC m^{-2}$) @endtex
152    REAL(r_std),DIMENSION(npts,nlitt,nvm,nbdl,nelements), INTENT(inout):: litter_below !! metabolic and structural litter, below
153                                                                                  !! ground @tex ($gC m^{-2}$) @endtex
154                                                                                  !@endtex
155    REAL(r_std),DIMENSION(npts,ncarb,nvm,nbdl,nelements), INTENT(inout) :: carbon           !! carbon pool: active, slow, or passive
156                                                                                  !! @tex ($gC m^{-2}$) @endtex
157    REAL(r_std), DIMENSION(npts,nvm,nbdl,ndoc,npool,nelements), INTENT(inout) :: DOC     !! Dissolved Organic Carbon in soil
158                                                                                       !! The unit is guven by m^3 of water
159                                                                                       !! @tex $(gC m{-3} of water)$ @endtex
160    !! 0.4 Local variables
161
162    INTEGER(i_std)                                            :: i, j, k, l, m    !! indices (unitless)
163    REAL(r_std),DIMENSION(npts,nelements)                     :: bm_new           !! biomass increase @tex ($gC m^{-2}$) @endtex
164    REAL(r_std),DIMENSION(npts,nparts,nelements)              :: biomass_loss     !! biomass loss @tex ($gC m^{-2}$) @endtex
165    REAL(r_std)                                               :: above            !! aboveground biomass @tex ($gC m^{-2}$) @endtex
166    REAL(r_std),DIMENSION(npts,nlitt,nelements)               :: dilu_lit_above   !! Above ground litter dilution @tex ($gC m^{-2}$) @endtex
167    REAL(r_std),DIMENSION(npts,nlitt,nbdl,nelements)          :: dilu_lit_below   !! Below ground litter dilution @tex ($gC m^{-2}$) @endtex
168    REAL(r_std),DIMENSION(npts,ncarb,nbdl,nelements)          :: dilu_soil_carbon !! Soil Carbondilution @tex ($gC m^{-2}$) @endtex
169    REAL(r_std),DIMENSION(nvm)                                :: delta_veg        !! changes in "maximal" coverage fraction of PFT
170    REAL(r_std)                                               :: delta_veg_sum    !! sum of delta_veg
171    REAL(r_std),DIMENSION(npts,nvm)                           :: delta_ind        !! change in number of individuals 
172    REAL(r_std),DIMENSION(0:nbdl)                             :: z_soil           !! Variable to store depth of the different
173                                                                                  !! soil layers (m)
174!_ ================================================================================================================================
175
176    IF (printlev>=3) WRITE(numout,*) 'Entering lcchange_main'
177   
178  !! 1. initialization
179   
180    prod10(:,0)         = zero
181    prod100(:,0)        = zero   
182    above               = zero
183    convflux(:)         = zero
184    cflux_prod10(:)     = zero
185    cflux_prod100(:)    = zero
186    delta_ind(:,:)      = zero
187    delta_veg(:)        = zero
188       !! 1.1.1 Copy the depth of the different soil layers (number of layers=nbdl)
189       !        previously calculated as variable diaglev in routines
190       !        sechiba.f90 and slowproc.f90 
191       z_soil(0) = zero
192       z_soil(1:nbdl) = diaglev(1:nbdl)
193   
194  !! 2. calculation of changes in carbon stocks and biomass by land cover change\n
195   
196    DO i = 1, npts ! Loop over # pixels - domain size
197       
198       !! 2.1 initialization of carbon stocks\n
199       delta_veg(:) = veget_max_new(i,:)-veget_max_old(i,:)
200       delta_veg_sum = SUM(delta_veg,MASK=delta_veg.LT.0.)
201       
202       dilu_lit_above(i,:,:) = zero
203       dilu_lit_below(i,:,:,:) = zero
204       dilu_soil_carbon(i,:,:,:) = zero
205       biomass_loss(i,:,:) = zero
206       
207       !! 2.2 if vegetation coverage decreases, compute dilution of litter, soil carbon, and biomass.\n
208
209
210       DO j=2, nvm
211          IF ( delta_veg(j) < -min_stomate ) THEN
212             dilu_lit_above(i,:,:) = dilu_lit_above(i,:,:) + delta_veg(j)*litter_above(i,:,j,:) / delta_veg_sum
213             dilu_lit_below(i,:,:,:) = dilu_lit_below(i,:,:,:) + delta_veg(j)*litter_below(i,:,j,:,:) / delta_veg_sum
214             DO k=2,nbdl
215               DO m = 1,npool
216               dilu_soil_carbon(i,:,k,icarbon) =  dilu_soil_carbon(i,:,k,icarbon) + delta_veg(j) * (carbon(i,:,j,k,icarbon) + &
217                                        (DOC(:,j,k,ifree,m,icarbon) + DOC(:,j,k,iadsorbed,m,icarbon))* moist_soil(:,k) * (z_soil(k)-z_soil(k-1))) / delta_veg_sum
218               dilu_soil_carbon(i,:,k,icarbon13) =  dilu_soil_carbon(i,:,k,icarbon13) + delta_veg(j) * (carbon(i,:,j,k,icarbon13) + &
219                                        (DOC(:,j,k,ifree,m,icarbon13) + DOC(:,j,k,iadsorbed,m,icarbon13))* moist_soil(:,k) * (z_soil(k)-z_soil(k-1))) / delta_veg_sum
220               ENDDO
221             ENDDO
222               DO m=1,npool
223               dilu_soil_carbon(i,:,1,icarbon) =  dilu_soil_carbon(i,:,1,icarbon) + delta_veg(j) * (carbon(i,:,j,1,icarbon)+ &
224                                        (DOC(:,j,1,ifree,m,icarbon) + DOC(:,j,1,iadsorbed,m,icarbon))* moist_soil(:,1) * (z_soil(1))) / delta_veg_sum
225               dilu_soil_carbon(i,:,1,icarbon13) =  dilu_soil_carbon(i,:,1,icarbon13) + delta_veg(j) * (carbon(i,:,j,1,icarbon13)+ &
226                                        (DOC(:,j,1,ifree,m,icarbon13) + DOC(:,j,1,iadsorbed,m,icarbon13))* moist_soil(:,1) * (z_soil(1))) / delta_veg_sum
227               ENDDO
228             biomass_loss(i,:,:) = biomass_loss(i,:,:) + biomass(i,j,:,:)*delta_veg(j) / delta_veg_sum
229          ENDIF
230       ENDDO
231       
232       !! 2.3
233!WRITE (numout,*) "litter_above_#1_lc_PFT6=",litter_above(:,:,6,:)
234!WRITE (numout,*) "litter_above_#1_lc_PFT10=",litter_above(:,:,10,:)
235       DO j=2, nvm ! Loop over # PFTs
236
237          !! 2.3.1 The case that vegetation coverage of PFTj increases
238          IF ( delta_veg(j) > min_stomate) THEN
239
240             !! 2.3.1.1 Initial setting of new establishment
241             IF (veget_max_old(i,j) .LT. min_stomate) THEN
242                IF (is_tree(j)) THEN
243
244                   ! cn_sapl(j)=0.5; stomate_data.f90
245                   cn_ind(i,j) = cn_sapl(j) 
246                ELSE
247                   cn_ind(i,j) = un
248                ENDIF
249                ind(i,j)= delta_veg(j) / cn_ind(i,j)
250                PFTpresent(i,j) = .TRUE.
251                everywhere(i,j) = 1.
252                senescence(i,j) = .FALSE.
253                age(i,j) = zero
254
255                ! large_value = 1.E33_r_std
256                when_growthinit(i,j) = large_value 
257                leaf_frac(i,j,1) = 1.0
258                npp_longterm(i,j) = npp_longterm_init
259                lm_lastyearmax(i,j) = bm_sapl(j,ileaf,icarbon) * ind(i,j)
260             ENDIF
261             IF ( cn_ind(i,j) > min_stomate ) THEN
262                delta_ind(i,j) = delta_veg(j) / cn_ind(i,j) 
263             ENDIF
264             
265             !! 2.3.1.2 Update of biomass in each each carbon stock component
266             !!         Update of biomass in each each carbon stock component (leaf, sapabove, sapbelow,
267             !>         heartabove, heartbelow, root, fruit, and carbres)\n
268             DO k = 1, nparts ! loop over # carbon stock components, nparts = 8; stomate_constant.f90
269                DO l = 1,nelements ! loop over # elements
270
271                   bm_new(i,l) = delta_ind(i,j) * bm_sapl(j,k,l) 
272                   IF (veget_max_old(i,j) .GT. min_stomate) THEN
273
274                      ! in the case that bm_new is overestimated compared with biomass?
275                      IF ((bm_new(i,l)/delta_veg(j)) > biomass(i,j,k,l)) THEN
276                         bm_new(i,l) = biomass(i,j,k,l)*delta_veg(j)
277                      ENDIF
278                   ENDIF
279                   biomass(i,j,k,l) = ( biomass(i,j,k,l) * veget_max_old(i,j) + bm_new(i,l) ) / veget_max_new(i,j)
280                   co2_to_bm(i,j) = co2_to_bm(i,j) + (bm_new(i,icarbon)* dt_days) / (one_year * veget_max_new(i,j))
281                END DO ! loop over # elements
282             ENDDO ! loop over # carbon stock components
283
284             !! 2.3.1.3 Calculation of dilution in litter, soil carbon, and  input of litter
285             !!        In this 'IF statement', dilu_* is zero. Formulas for litter and soil carbon
286             !!         could be shortend?? Are the following formulas correct?
287
288             ! Litter
289             litter_above(i,:,j,:)=(litter_above(i,:,j,:) * veget_max_old(i,j) + &
290                  dilu_lit_above(i,:,:) * delta_veg(j)) / veget_max_new(i,j)
291           
292             ! Soil carbon
293             carbon(i,:,j,:,:)=(carbon(i,:,j,:,:) * veget_max_old(i,j) + dilu_soil_carbon(i,:,:,:) * delta_veg(j)) / veget_max_new(i,j)
294
295             DO l = 1,nelements
296
297                ! Litter input
298                bm_to_litter(i,j,isapbelow,l) = bm_to_litter(i,j,isapbelow,l) + &
299                     &                         biomass_loss(i,isapbelow,l)*delta_veg(j) / veget_max_new(i,j)
300                bm_to_litter(i,j,iheartbelow,l) = bm_to_litter(i,j,iheartbelow,l) + biomass_loss(i,iheartbelow,l) *delta_veg(j) &
301                     &  / veget_max_new(i,j)
302                bm_to_litter(i,j,iroot,l) = bm_to_litter(i,j,iroot,l) + biomass_loss(i,iroot,l)*delta_veg(j) / veget_max_new(i,j)
303                bm_to_litter(i,j,ifruit,l) = bm_to_litter(i,j,ifruit,l) + biomass_loss(i,ifruit,l)*delta_veg(j) / veget_max_new(i,j)
304                bm_to_litter(i,j,icarbres,l) = bm_to_litter(i,j,icarbres,l) + &
305                     &                         biomass_loss(i,icarbres,l)   *delta_veg(j) / veget_max_new(i,j)
306                bm_to_litter(i,j,ileaf,l) = bm_to_litter(i,j,ileaf,l) + biomass_loss(i,ileaf,l)*delta_veg(j) / veget_max_new(i,j)
307
308             END DO
309
310             age(i,j)=age(i,j)*veget_max_old(i,j)/veget_max_new(i,j)
311             
312          !! 2.3.2 The case that vegetation coverage of PFTj is no change or decreases
313          ELSE 
314 
315             !! 2.3.2.1 Biomass export
316             ! coeff_lcchange_*:  Coeff of biomass export for the year, decade, and century
317             above = biomass(i,j,isapabove,icarbon) + biomass(i,j,iheartabove,icarbon)
318             convflux(i)  = convflux(i)  - ( coeff_lcchange_1(j) * above * delta_veg(j) ) 
319             prod10(i,0)  = prod10(i,0)  - ( coeff_lcchange_10(j) * above * delta_veg(j) )
320             prod100(i,0) = prod100(i,0) - ( coeff_lcchange_100(j) * above * delta_veg(j) )
321
322
323             !! 2.3.2.2 Total reduction
324             !! If the vegetation is to small, it has been set to 0.
325             IF ( veget_max_new(i,j) .LT. min_stomate ) THEN
326               
327                ind(i,j) = zero
328                biomass(i,j,:,:) = zero
329                PFTpresent(i,j) = .FALSE.
330                senescence(i,j) = .FALSE.
331                age(i,j) = zero
332                when_growthinit(i,j) = undef
333                everywhere(i,j) = zero
334                carbon(i,:,j,:,:) = zero
335                litter_above(i,:,j,:) = zero
336                litter_below(i,:,j,:,:) = zero
337                bm_to_litter(i,j,:,:) = zero
338                turnover_daily(i,j,:,:) = zero
339               
340             ENDIF
341 
342          ENDIF ! End if PFT's coverage reduction
343         
344       ENDDO ! Loop over # PFTs
345!WRITE (numout,*) "litter_above_#2_lc_PFT6=",litter_above(:,:,6,:)
346!WRITE (numout,*) "litter_above_#2_lc_PFT10=",litter_above(:,:,10,:)
347       
348       !! 2.4 update 10 year-turnover pool content following flux emission
349       !!     (linear decay (10%) of the initial carbon input)
350       DO  l = 0, 8
351          m = 10 - l
352          cflux_prod10(i) =  cflux_prod10(i) + flux10(i,m)
353          prod10(i,m)     =  prod10(i,m-1)   - flux10(i,m-1)
354          flux10(i,m)     =  flux10(i,m-1)
355         
356          IF (prod10(i,m) .LT. 1.0) prod10(i,m) = zero
357       ENDDO
358       
359       cflux_prod10(i) = cflux_prod10(i) + flux10(i,1) 
360       flux10(i,1)     = 0.1 * prod10(i,0)
361       prod10(i,1)     = prod10(i,0)
362       
363       !! 2.5 update 100 year-turnover pool content following flux emission\n
364       DO   l = 0, 98
365          m = 100 - l
366          cflux_prod100(i)  =  cflux_prod100(i) + flux100(i,m)
367          prod100(i,m)      =  prod100(i,m-1)   - flux100(i,m-1)
368          flux100(i,m)      =  flux100(i,m-1)
369         
370          IF (prod100(i,m).LT.1.0) prod100(i,m) = zero
371       ENDDO
372       
373       cflux_prod100(i)  = cflux_prod100(i) + flux100(i,1) 
374       flux100(i,1)      = 0.01 * prod100(i,0)
375       prod100(i,1)      = prod100(i,0)
376       prod10(i,0)        = zero
377       prod100(i,0)       = zero 
378       
379    ENDDO ! Loop over # pixels - domain size
380   
381  !! 3. history
382    convflux        = convflux/one_year*dt_days
383    cflux_prod10    = cflux_prod10/one_year*dt_days
384    cflux_prod100   = cflux_prod100/one_year*dt_days
385   
386    IF (printlev>=4) WRITE(numout,*) 'Leaving lcchange_main'
387   
388  END SUBROUTINE lcchange_main
389 
390END MODULE stomate_lcchange
Note: See TracBrowser for help on using the repository browser.