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

Last change on this file since 8579 was 7327, checked in by josefine.ghattas, 3 years ago

Bug correction for wood harvest. See ticket #781
File size: 10.4 KB
Line 
1! =================================================================================================================================
2! MODULE       : stomate_woodharvest
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       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: svn://forge.ipsl.jussieu.fr/orchidee/trunk/ORCHIDEE/src_stomate/stomate_woodharvest.f90 $
19!! $Date: 2016-01-06 13:15:53 +0100 (mer. 06 janv. 2016) $
20!! $Revision: 3094 $
21!! \n
22!_ ================================================================================================================================
23
24
25MODULE stomate_woodharvest
26
27  USE ioipsl_para
28  USE stomate_data
29  USE pft_parameters
30  USE constantes
31 
32  IMPLICIT NONE
33 
34  PRIVATE
35  PUBLIC woodharvest_main
36 
37CONTAINS
38
39
40!! ================================================================================================================================
41!! SUBROUTINE   : woodharvest_main
42!!
43!>\BRIEF        Impact of wood harvest on carbon stocks
44!!
45!! DESCRIPTION  : This subroutine is always activate if DO_WOOD_HARVEST=y in the configuration file.
46!! woodharvest_main is called from stomateLpj every year (at the beginning)
47!! The impact of wood harvest on carbon stocks is computed in this subroutine. The land cover change is written
48!! by the difference of current and previous "maximal" coverage fraction of a PFT.
49!! On the basis of this difference, the amount of 'new establishment'/'biomass export',
50!! and increase/decrease of each component, are estimated.\n
51!!
52!! Main structure of woodharvest is:
53!! 1. Initialization
54!! 2. Calculation of changes in carbon stocks and biomass by wood harvest
55!! 3. Update 10 year- and 100 year-turnover pool contents
56!! 4. History
57!!
58!! RECENT CHANGE(S) : None
59!!
60!! MAIN OUTPUT VARIABLE(S) : ::prod10, ::prod100, ::flux10, ::flux100,
61!!   :: cflux_prod10 and :: cflux_prod100
62!!
63!! REFERENCES   : None
64!!
65!! FLOWCHART    :
66!! \n
67!_ ================================================================================================================================
68
69 
70  SUBROUTINE woodharvest_main ( npts, dt_days, veget_max_old, &
71       biomass, &       
72       flux10_harvest,flux100_harvest, prod10_harvest,prod100_harvest,&
73       convflux_harvest,cflux_prod10_harvest,cflux_prod100_harvest,&
74       harvestwood,harvestpft,fHarvestToProduct)
75   
76    IMPLICIT NONE
77   
78  !! 0. Variable and parameter declaration
79   
80    !! 0.1 Input variables
81   
82    INTEGER, INTENT(in)                            :: npts             !! Domain size - number of pixels (unitless)
83    REAL(r_std), INTENT(in)                        :: dt_days          !! Time step of vegetation dynamics for stomate (days)
84    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)   :: veget_max_old    !! Current "maximal" coverage fraction of a PFT (LAI
85                                                                       !! -> infinity) on ground
86    REAL(r_std), DIMENSION(npts), INTENT(in)       :: harvestwood      !! Harvested wood biomass (gC m-2 yr-1)
87 
88    !! 0.2 Output variables
89
90    REAL(r_std), DIMENSION(npts), INTENT(out)      :: convflux_harvest        !! flux release during first year following wood harvest
91    REAL(r_std), DIMENSION(npts), INTENT(out)      :: cflux_prod10_harvest    !! total annual release from the 10 year-turnover
92                                                                              !! pool @tex ($gC m^{-2}$) @endtex
93    REAL(r_std), DIMENSION(npts), INTENT(out)      :: cflux_prod100_harvest   !! total annual release from the 100 year-
94                                                                              !! turnover pool @tex ($gC m^{-2}$) @endtex
95    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)  :: harvestpft              !! wood harvested for each PFT
96                                                                              !! @tex ($gC m^{-2} day^{-1}$) @endtex
97    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)  :: fHarvestToProduct       !!  Harvested biomass into product pool due to anthorpogenic
98                                                                              !! land use (gC m-2 )
99!! 0.3 Modified variables   
100   
101    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass    !! biomass @tex ($gC m^{-2}$) @endtex
102    REAL(r_std), DIMENSION(npts,0:10), INTENT(inout)  :: prod10_harvest  !! products remaining in the 10/100 year-turnover pool
103                                                                         !! after the annual release for each compartment
104                                                                         !! (10 or 100 + 1 : input from year of land cover change)
105    REAL(r_std), DIMENSION(npts,0:100), INTENT(inout) :: prod100_harvest !! see prod10_havest
106    REAL(r_std), DIMENSION(npts,10), INTENT(inout)    :: flux10_harvest  !! annual release from the 10 year-turnover pool compartments
107
108    REAL(r_std), DIMENSION(npts,100), INTENT(inout)   :: flux100_harvest !! annual release from the 100 year-turnover pool compartments
109                                                   
110    !! 0.4 Local variables
111
112    INTEGER(i_std)            :: i, j, l, m    !! indices (unitless)
113    REAL(r_std)               :: harvest       !! wood harvest renomalized by forest fraction
114    REAL(r_std)               :: relharvest    !! relative harvest fraction of each pft
115    REAL(r_std)               :: relbiomass    !! relative biomass of each PFT
116    REAL(r_std)               :: vegettree     !! vegfrac of trees
117    REAL(r_std)               :: abovemean     !! mean aboveground biomass @tex ($gC m^{-2}$) @endtex
118!_ ================================================================================================================================
119
120    IF (printlev>=3) WRITE(numout,*) 'Entering woodharvest_main'
121   
122  !! 1. initialization
123   
124    !! 2. calculation of changes in carbon stocks and biomass by wood harvest \n
125    !! 2.1 initialization of carbon stocks\n
126    prod10_harvest(:,0)           = zero
127    prod100_harvest(:,0)          = zero   
128    convflux_harvest(:)           = zero
129    cflux_prod10_harvest(:)       = zero
130    cflux_prod100_harvest(:)      = zero
131    harvestpft(:,:)               = zero
132
133    fHarvestToProduct(:,:)        = zero
134
135    DO i = 1, npts 
136       abovemean=0         
137       vegettree=0
138       !! 2.3 calculation of total forest fraction and mean biomass
139       DO j=1, nvm
140          if (is_tree(j)) THEN
141             abovemean=abovemean+((biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon)))*veget_max_old(i,j)
142             vegettree=vegettree+veget_max_old(i,j)
143          ENDIF
144       ENDDO
145       IF (vegettree > min_stomate) THEN
146         
147          harvest=harvestwood(i) 
148          abovemean=abovemean/vegettree
149
150          !! 2.4 wood harvest for each PFT is estimated proportionally to the relative biomass.
151          DO j=1,nvm
152             relharvest=0
153             IF ((is_tree(j)) .AND. (abovemean > min_stomate)  .AND. &
154                  (biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon)>min_stomate))  THEN
155                relbiomass=((biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon))/abovemean) 
156                relharvest=harvest*relbiomass
157
158                !! 2.5 the harvest could not be higher than 20% of the biomasse to avoid rapid killing of vegetation
159                IF ((biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon))*0.2<relharvest) &
160                     relharvest=(biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon))*0.2; 
161                biomass(i,j,iheartabove,:)= biomass(i,j,iheartabove,:) - &
162                     relharvest*(biomass(i,j,iheartabove,:)/(biomass(i,j,iheartabove,:)+biomass(i,j,isapabove,:)))
163                biomass(i,j,isapabove,:)= biomass(i,j,isapabove,:) - &
164                     relharvest*(biomass(i,j,isapabove,:)/(biomass(i,j,iheartabove,:)+biomass(i,j,isapabove,:)))
165
166                fHarvestToProduct(i,j)=relharvest*veget_max_old(i,j)
167                ! 2.6 calculation of different flux like for LU
168                convflux_harvest(i)  = convflux_harvest(i) +coeff_lcchange_1(j)*relharvest*veget_max_old(i,j)
169                harvestpft(i,j)=relharvest*dt_days/one_year
170                prod10_harvest(i,0)  = prod10_harvest(i,0)  +coeff_lcchange_10(j)*relharvest*veget_max_old(i,j)
171                prod100_harvest(i,0) = prod100_harvest(i,0) +coeff_lcchange_100(j)*relharvest*veget_max_old(i,j)
172             ENDIF
173          ENDDO
174       ENDIF
175
176
177
178       !! 3.4 update 10 year-turnover pool content following flux emission
179       !!     (linear decay (10%) of the initial carbon input)
180       DO  l = 0, 8
181          m = 10 - l
182          cflux_prod10_harvest(i) =  cflux_prod10_harvest(i) + flux10_harvest(i,m)
183          prod10_harvest(i,m)   =  prod10_harvest(i,m-1)   - flux10_harvest(i,m-1)
184          flux10_harvest(i,m)   =  flux10_harvest(i,m-1)
185          IF (prod10_harvest(i,m) .LT. 1.0) prod10_harvest(i,m) = zero
186       ENDDO
187
188
189       cflux_prod10_harvest(i) = cflux_prod10_harvest(i) + flux10_harvest(i,1) 
190       flux10_harvest(i,1)     = 0.1 * prod10_harvest(i,0)
191       prod10_harvest(i,1)     = prod10_harvest(i,0)
192       
193       !! 3.5 update 100 year-turnover pool content following flux emission\n
194       DO   l = 0, 98
195          m = 100 - l
196          cflux_prod100_harvest(i)  =  cflux_prod100_harvest(i) + flux100_harvest(i,m)
197          prod100_harvest(i,m)      =  prod100_harvest(i,m-1)   - flux100_harvest(i,m-1)
198          flux100_harvest(i,m)      =  flux100_harvest(i,m-1)
199          IF (prod100_harvest(i,m).LT.1.0) prod100_harvest(i,m) = zero
200       ENDDO
201       
202       cflux_prod100_harvest(i)  = cflux_prod100_harvest(i) + flux100_harvest(i,1) 
203       flux100_harvest(i,1)      = 0.01 * prod100_harvest(i,0)
204       prod100_harvest(i,1)      = prod100_harvest(i,0)
205       prod10_harvest(i,0)        = 0.0
206       prod100_harvest(i,0)       = 0.0 
207       
208    ENDDO ! Loop over # pixels - domain size
209   
210  !! 4. history
211    convflux_harvest        = convflux_harvest/one_year*dt_days
212    cflux_prod10_harvest    = cflux_prod10_harvest/one_year*dt_days
213    cflux_prod100_harvest   = cflux_prod100_harvest/one_year*dt_days
214    fHarvestToProduct       = fHarvestToProduct/one_year*dt_days   
215    IF (printlev>=4) WRITE(numout,*) 'Leaving woodharvest_main'
216   
217  END SUBROUTINE woodharvest_main
218 
219END MODULE stomate_woodharvest
Note: See TracBrowser for help on using the repository browser.