source: tags/ORCHIDEE_4_1/ORCHIDEE/src_stomate/stomate_woodharvest.f90 @ 7852

Last change on this file since 7852 was 5639, checked in by josefine.ghattas, 6 years ago

Nitrogen is in the trunk!

  • Copied branches/ORCHIDEE-CN at revision 5638 into ORCHIDEE
  • Copied branches/ORCHIDEE-CN_CONFIG at revision 5637 in ORCHIDEE_OL

See ticket #469

  1. Vuichard, J Ghattas
File size: 11.7 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,nflux_prod_total_harvest)
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    REAL(r_std), DIMENSION(npts), INTENT(out)      :: nflux_prod_total_harvest !! release of N associated to wood harvest  @tex ($gN m^{-2}$) @endtex   
100!! 0.3 Modified variables   
101   
102    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass    !! biomass @tex ($gC m^{-2}$) @endtex
103    REAL(r_std), DIMENSION(npts,0:10), INTENT(inout)  :: prod10_harvest  !! products remaining in the 10/100 year-turnover pool
104                                                                         !! after the annual release for each compartment
105                                                                         !! (10 or 100 + 1 : input from year of land cover change)
106    REAL(r_std), DIMENSION(npts,0:100), INTENT(inout) :: prod100_harvest !! see prod10_havest
107    REAL(r_std), DIMENSION(npts,10), INTENT(inout)    :: flux10_harvest  !! annual release from the 10 year-turnover pool compartments
108
109    REAL(r_std), DIMENSION(npts,100), INTENT(inout)   :: flux100_harvest !! annual release from the 100 year-turnover pool compartments
110                                               
111    !! 0.4 Local variables
112
113    INTEGER(i_std)            :: i, j, l, m    !! indices (unitless)
114    REAL(r_std)               :: harvest       !! wood harvest renomalized by forest fraction
115    REAL(r_std)               :: relharvest    !! relative harvest fraction of each pft
116    REAL(r_std)               :: relbiomass    !! relative biomass of each PFT
117    REAL(r_std)               :: vegettree     !! vegfrac of trees
118    REAL(r_std)               :: abovemean     !! mean aboveground biomass @tex ($gC m^{-2}$) @endtex
119    REAL(r_std)               :: remove_biomass_heartabove_carbon  !! Removed carbon biomass from heartabove pool
120    REAL(r_std)               :: remove_biomass_sapabove_carbon    !! Removed carbon biomass from sapabove pool
121    REAL(r_std)               :: remove_biomass_heartabove_nitrogen!! Removed nitrogen biomass from heartabove pool
122    REAL(r_std)               :: remove_biomass_sapabove_nitrogen  !! Removed nitrogen biomass from sapabove pool
123   
124!_ ================================================================================================================================
125
126    IF (printlev>=3) WRITE(numout,*) 'Entering woodharvest_main'
127   
128  !! 1. initialization
129   
130    !! 2. calculation of changes in carbon stocks and biomass by wood harvest \n
131    !! 2.1 initialization of carbon stocks\n
132    prod10_harvest(:,0)           = zero
133    prod100_harvest(:,0)          = zero   
134    convflux_harvest(:)           = zero
135    cflux_prod10_harvest(:)       = zero
136    cflux_prod100_harvest(:)      = zero
137    harvestpft(:,:)               = zero
138
139    fHarvestToProduct(:,:)        = zero
140    nflux_prod_total_harvest(:)   = zero
141
142    DO i = 1, npts 
143       abovemean=0         
144       vegettree=0
145       !! 2.3 calculation of total forest fraction and mean biomass
146       DO j=1, nvm
147          if (is_tree(j)) THEN
148             abovemean=abovemean+((biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon)))*veget_max_old(i,j)
149             vegettree=vegettree+veget_max_old(i,j)
150          ENDIF
151       ENDDO
152       IF (vegettree > min_stomate) THEN
153         
154          harvest=harvestwood(i) 
155          abovemean=abovemean/vegettree
156
157          !! 2.4 wood harvest for each PFT is estimated proportionally to the relative biomass.
158          DO j=1,nvm
159             relharvest=0
160             IF ((is_tree(j)) .AND. (abovemean > min_stomate)  .AND. &
161                  (biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon)>min_stomate))  THEN
162                relbiomass=((biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon))/abovemean) 
163                relharvest=harvest*relbiomass
164
165                !! 2.5 the harvest could not be higher than 20% of the biomasse to avoid rapid killing of vegetation
166                IF ((biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon))*0.2<relharvest) &
167                     relharvest=(biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon))*0.2; 
168               
169                remove_biomass_heartabove_carbon = relharvest*(biomass(i,j,iheartabove,icarbon) / &
170                     (biomass(i,j,iheartabove,icarbon)+biomass(i,j,isapabove,icarbon)))
171                remove_biomass_heartabove_nitrogen = remove_biomass_heartabove_carbon * biomass(i,j,iheartabove,initrogen) / &
172                     biomass(i,j,iheartabove,icarbon)
173                remove_biomass_sapabove_carbon = relharvest*(biomass(i,j,isapabove,icarbon) / &
174                     (biomass(i,j,iheartabove,icarbon) + biomass(i,j,isapabove,icarbon)))
175                remove_biomass_sapabove_nitrogen = remove_biomass_sapabove_carbon * biomass(i,j,isapabove,initrogen) / &
176                     biomass(i,j,isapabove,icarbon)
177               
178               
179                biomass(i,j,iheartabove,icarbon)= biomass(i,j,iheartabove,icarbon) - remove_biomass_heartabove_carbon
180                biomass(i,j,isapabove,icarbon)= biomass(i,j,isapabove,icarbon) - remove_biomass_sapabove_carbon
181                biomass(i,j,iheartabove,initrogen)= biomass(i,j,iheartabove,initrogen) - remove_biomass_heartabove_nitrogen
182                biomass(i,j,isapabove,initrogen)= biomass(i,j,isapabove,initrogen) - remove_biomass_sapabove_nitrogen
183
184                fHarvestToProduct(i,j)=relharvest*veget_max_old(i,j)
185                ! 2.6 calculation of different flux like for LU
186                convflux_harvest(i)  = convflux_harvest(i) +coeff_lcchange_1(j)*relharvest*veget_max_old(i,j)
187                harvestpft(i,j)=relharvest*dt_days/one_year
188                prod10_harvest(i,0)  = prod10_harvest(i,0)  +coeff_lcchange_10(j)*relharvest*veget_max_old(i,j)
189                prod100_harvest(i,0) = prod100_harvest(i,0) +coeff_lcchange_100(j)*relharvest*veget_max_old(i,j)
190
191
192                nflux_prod_total_harvest(i) = nflux_prod_total_harvest(i) + &
193                     (remove_biomass_heartabove_nitrogen+remove_biomass_sapabove_nitrogen)*veget_max_old(i,j)
194             ENDIF
195          ENDDO
196       ENDIF
197
198       cflux_prod10_harvest(i) = cflux_prod10_harvest(i) + flux10_harvest(i,1) 
199       flux10_harvest(i,1)     = 0.1 * prod10_harvest(i,0)
200       prod10_harvest(i,1)     = prod10_harvest(i,0)
201       
202       !! 3.5 update 100 year-turnover pool content following flux emission\n
203       DO   l = 0, 98
204          m = 100 - l
205          cflux_prod100_harvest(i)  =  cflux_prod100_harvest(i) + flux100_harvest(i,m)
206          prod100_harvest(i,m)      =  prod100_harvest(i,m-1)   - flux100_harvest(i,m-1)
207          flux100_harvest(i,m)      =  flux100_harvest(i,m-1)
208          IF (prod100_harvest(i,m).LT.1.0) prod100_harvest(i,m) = zero
209       ENDDO
210       
211       cflux_prod100_harvest(i)  = cflux_prod100_harvest(i) + flux100_harvest(i,1) 
212       flux100_harvest(i,1)      = 0.01 * prod100_harvest(i,0)
213       prod100_harvest(i,1)      = prod100_harvest(i,0)
214       prod10_harvest(i,0)        = 0.0
215       prod100_harvest(i,0)       = 0.0 
216       
217    ENDDO ! Loop over # pixels - domain size
218   
219  !! 4. history
220    convflux_harvest        = convflux_harvest/one_year*dt_days
221    cflux_prod10_harvest    = cflux_prod10_harvest/one_year*dt_days
222    cflux_prod100_harvest   = cflux_prod100_harvest/one_year*dt_days
223    fHarvestToProduct       = fHarvestToProduct/one_year*dt_days   
224    nflux_prod_total_harvest = nflux_prod_total_harvest/one_year*dt_days
225
226    IF (printlev>=4) WRITE(numout,*) 'Leaving woodharvest_main'
227   
228  END SUBROUTINE woodharvest_main
229 
230END MODULE stomate_woodharvest
Note: See TracBrowser for help on using the repository browser.