source: branches/publications/ORCHIDEE_CAN_r3069/src_stomate/sapiens_kill.f90 @ 7475

Last change on this file since 7475 was 3069, checked in by sebastiaan.luyssaert, 9 years ago

PROD: tested gloabally for the zoomed European grid for 20 years. Do not use the previous commit as there were problems with svn it does not contain the described changes. This commit contains the bug fixes to species change and added functionality to change forest management following mortality or a harvest.

File size: 93.6 KB
Line 
1! =================================================================================================================================
2! MODULE       : sapiens_kill
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       Kill all vegetation scheduled for killing by human activities.
10!!
11!!\n DESCRIPTION : Actually kills biomass.  The biomass is scheduled for killing
12!!       by the variable circ_class_kill in other sapiens routines.  Here the
13!!       biomass is moved to the litter and harvest pools, as appropriate.
14!!
15!! RECENT CHANGE(S): None
16!!
17!! REFERENCE(S) :
18!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
19!!         plant geography and terrestrial carbon cycling in the LPJ dynamic
20!!         global vegetation model, Global Change Biology, 9, 161-185.\n
21!! - Waring, R. H. (1983). "Estimating forest growth and efficiency in relation
22!!         to canopy leaf area." Advances in Ecological Research 13: 327-354.\n
23!!
24!! SVN          :
25!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/lpj_gap.f90 $
26!! $Date: 2013-08-19 17:50:44 +0200 (Mon, 19 Aug 2013) $
27!! $Revision: 1425 $
28!! \n
29!_ ================================================================================================================================
30
31MODULE sapiens_kill
32
33  ! modules used:
34
35  USE stomate_data
36  USE pft_parameters
37  USE ioipsl_para 
38  USE function_library,       ONLY : distribute_mortality_biomass,&
39                                     nmax,  wood_to_dia
40  USE constantes
41
42  IMPLICIT NONE
43
44  ! private & public routines
45
46  PRIVATE
47  PUBLIC gap_prog_sapiens
48 
49  ! Variable declaration
50
51CONTAINS
52
53
54!! ================================================================================================================================
55!! SUBROUTINE   : gap_prog_sapiens
56!!
57!>\BRIEF        Transfer dead biomass to litter and update stand density as a result
58!!              of human activities.
59!!
60!! DESCRIPTION  : This routine has a simple purpose: kill individuals by human activities. 
61!!   The variable circ_class_kill indicates how many individuals need to be killed by the
62!!   various processes in the code, and is calculated in other modules.  gap_prog_sapiens
63!!   takes this variable and decides on the correct order of which to actually kill the
64!!   trees.  As these are all human activities, some of the biomass may be moved to
65!!   harvest pools, while other parts are moved to litter.
66!!
67!! RECENT CHANGE(S):
68!!
69!! MAIN OUTPUT VARIABLE(S): ::circ_class_biomass; ::biomass; ::ind
70!!
71!! REFERENCE(S)   :
72!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
73!!         plant geography and terrestrial carbon cycling in the LPJ dynamic
74!!         global vegetation model, Global Change Biology, 9, 161-185.
75!! - Waring, R. H. (1983). "Estimating forest growth and efficiency in relation
76!!         to canopy leaf area." Advances in Ecological Research 13: 327-354.
77!!
78!! FLOWCHART    : None
79!!\n
80!_ ================================================================================================================================
81 
82  SUBROUTINE gap_prog_sapiens (npts, biomass, ind,  &
83       bm_to_litter, circ_class_biomass, circ_class_kill, &
84       circ_class_n, harvest_pool_bound, harvest_pool, harvest_type, &
85       harvest_cut, harvest_area, resolution, &
86       veget_max, age_stand, last_cut, mai_count, circ_class_dist, &
87       coppice_dens)
88
89    !! 0. Variable and parameter declaration
90
91    !! 0.1 Input variables
92    ! NOTE: if harvest_pool_bound is not explicitely defined
93    ! the model uses 1:ndia_harvest+2 as its dimensions. If
94    ! you like to do so remember to change the counters in
95    ! this routine
96    INTEGER(i_std), INTENT(in)                                   :: npts               !! Domain size (-)
97    REAL(r_std), DIMENSION(0:ndia_harvest+1), INTENT(in)         :: harvest_pool_bound !! The boundaries of the diameter classes
98
99                                                                                       !! in the wood harvest pools
100                                                                                       !! @tex $(m)$ @endtex
101    REAL(r_std),DIMENSION(npts,2), INTENT(in)                    :: resolution         !! The length of each grid-box in X and Y
102                                                                                       !! @tex $(m)$ @endtex
103    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: veget_max          !! Maximum fraction of vegetation type
104    REAL(r_std), DIMENSION(ncirc), INTENT(in)                    :: circ_class_dist    !! The probability distribution of trees
105                                                                                       !! in a circ class in case of a
106                                                                                       !! redistribution (unitless).
107    REAL(r_std), DIMENSION(npts,nvm),INTENT(in)                  :: coppice_dens       !! The density of a coppice at the first
108                                                                                       !! cutting.
109                                                                                       !! @tex $( 1 / m**2 )$ @endtex
110
111    !! 0.2 Output variables
112
113    !! 0.3 Modified variables   
114    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
115                                             INTENT(inout)       :: biomass            !! Biomass @tex $(gC m^{-2}) $@endtex
116    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)              :: ind                !! Stand level number of individuals
117                                                                                       !! @tex $(m^{-2})$ @endtex
118    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &                         
119                                              INTENT(inout)      :: circ_class_biomass !! Biomass of the components of the model 
120                                                                                       !! tree within a circumference
121                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
122    REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout)        :: circ_class_n       !! Number of individuals in each circ class
123                                                                                       !! @tex $(m^{-2})$ @endtex
124    REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times),&
125                                                  INTENT(inout)  :: circ_class_kill    !! Number of trees within a circ that needs
126                                                                                       !! to be killed @tex $(ind m^{-2})$ @endtex
127    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
128                                                  INTENT(inout)  :: bm_to_litter       !! Biomass transfer to litter
129                                                                                       !! @tex $(gC m^{-2})$ @endtex
130    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)               :: harvest_pool       !! The wood and biomass that have been
131                                                                                       !! havested by humans @tex $(gC)$ @endtex
132    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_type       !! Type of management that resulted
133                                                                                       !! in the harvest (unitless)
134    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_cut        !! Type of cutting that was used for the harvest
135                                                                                       !! (unitless)
136    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_area       !! Harvested area (m^{2})
137    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: last_cut           !! Years since last thinning or clearcut (years)
138    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: age_stand          !! Age of stand (years)
139    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: mai_count          !! The number of times we've
140                                                                                       !! calculated the volume increment
141                                                                                       !! for a stand
142 
143    !! 0.4 Local variables
144    INTEGER(i_std)                                               :: ipts,ivm           !! Indices
145    INTEGER(i_std)                                               :: iele               !! Indices
146    INTEGER(i_std)                                               :: ipar, imbc         !! Indices
147   
148    REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements)           :: check_intern       !! Contains the components of the internal
149                                                                                       !! mass balance chech for this routine
150                                                                                       !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
151    REAL(r_std), DIMENSION(npts,nvm,nelements)                   :: closure_intern     !! Check closure of internal mass balance
152                                                                                       !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
153    REAL(r_std), DIMENSION(npts,nvm,nelements)                   :: pool_start         !! Start pool of this routine
154                                                                                       !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
155    REAL(r_std), DIMENSION(npts,nvm,nelements)                   :: pool_end           !! End pool of this routine
156                                                                                       !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
157!_ ================================================================================================================================
158
159    IF (bavard.GE.2) WRITE(numout,*) 'Entering gap_prog_sapiens'
160
161!! 2. Initialize variables
162
163    !! 2.1 Initialize check for mass balance closure
164    !  The mass balance is calculated at the end of this routine
165    !  in section .
166    pool_start = zero
167    DO ipar = 1,nparts
168       DO iele = 1,nelements
169          !  Initial litter and biomass
170          pool_start(:,:,iele) = pool_start(:,:,iele) + &
171               (biomass(:,:,ipar,iele) * veget_max(:,:)) + &
172               (bm_to_litter(:,:,ipar,iele) * veget_max(:,:))
173       ENDDO
174    ENDDO
175
176    ! The biomass harvest pool shouldn't be multiplied by veget_max
177    ! its units are already in gC pixel-1
178    ! Right now the harvest pool is in gC.  We will divide by the
179    ! resolution to make it per pixel.
180    DO ipts = 1, npts
181       pool_start(ipts,:,icarbon) = pool_start(ipts,:,icarbon) + &
182         SUM(harvest_pool(ipts,:,:,icarbon),2) / &
183         (resolution(ipts,1) * resolution(ipts,2))
184    ENDDO
185
186    !! 2.2 Initialize check for surface area conservation
187    !  Veget_max is a INTENT(in) variable and can therefore
188    !  not be changed during the course of this subroutine
189    !  No need to check whether the subroutine preserves the
190    !  total surface area of the pixel.
191
192
193 !! 3. Kill plants
194
195    ! All trees marked to die from forest management or salavage looging will
196    ! be killed here (trees that died under a conservation strategu are dealt
197    ! with sapiens_gap_prognostic. 
198    DO  ipts = 1,npts ! loop over land_points
199
200       DO ivm = 2,nvm  ! loop over #PFT
201
202          IF(veget_max(ipts,ivm) == zero)THEN
203              ! this vegetation type is not present, so no reason to do the
204              ! calculation
205              CYCLE
206          ENDIF
207
208          DO iele = 1,nelements
209
210             ! Clearcuts and thinning can only happen on trees.
211             IF(is_tree(ivm))THEN
212
213                ! Now see if there are any plants to be killed by a clearcut.
214                ! Only the stems are harvested.  This should be done before
215                ! a thinning so that we don't overcount.
216                CALL clearcut_harvest_aboveground_dead_roots(npts, ipts, ivm, iele, &
217                     ifm_thin, icut_clear, biomass, ind, bm_to_litter, &
218                     circ_class_biomass, circ_class_kill, circ_class_n, harvest_pool, &
219                     harvest_type, harvest_cut, harvest_area, &
220                     harvest_pool_bound, resolution, veget_max, last_cut, age_stand, &
221                     mai_count, (un - branch_harvest(ivm)*branch_ratio(ivm)))
222
223                ! Not only do we harvest stems during thinning, we have to adjust
224                ! the mortality.  We'll do all this in a subroutine.
225                CALL thinning_harvest_aboveground_dead_roots(npts, ipts, ivm, iele, &
226                     ifm_thin, icut_thin, biomass, ind, bm_to_litter, &
227                     circ_class_biomass, circ_class_kill, circ_class_n, harvest_pool, &
228                     harvest_type, harvest_cut, harvest_area, &
229                     harvest_pool_bound, resolution, veget_max, last_cut, &
230                     (un - branch_harvest(ivm)*branch_ratio(ivm)))
231
232                ! And the final cut of an SRC, where the underground wood dies.  This is
233                ! exactly the same as the clearcut above, except that we take all
234                ! of the aboveground biomass.
235                CALL clearcut_harvest_aboveground_dead_roots(npts, ipts, ivm, iele, &
236                     ifm_src, icut_cop3, biomass, ind, bm_to_litter, &
237                     circ_class_biomass, circ_class_kill, circ_class_n, harvest_pool, &
238                     harvest_type, harvest_cut, harvest_area, &
239                     harvest_pool_bound, resolution, veget_max, last_cut, age_stand, &
240                     mai_count, un)
241
242               ! Have we cut a coppice system for the first time?
243                CALL first_coppice_cut(npts, ipts, ivm, iele, ifm_cop, &
244                     icut_cop1, harvest_pool_bound, resolution, veget_max, biomass, &
245                     ind, circ_class_biomass, circ_class_n, circ_class_kill, &
246                     bm_to_litter, harvest_pool, harvest_type, harvest_cut, &
247                     harvest_area, last_cut)
248
249                CALL first_coppice_cut(npts, ipts, ivm, iele, ifm_src, &
250                     icut_cop1, harvest_pool_bound, resolution, veget_max, biomass, &
251                     ind, circ_class_biomass, circ_class_n, circ_class_kill, &
252                     bm_to_litter, harvest_pool, harvest_type, harvest_cut, &
253                     harvest_area, last_cut)
254
255                ! What if we coppice for the second (and subsequent) cuts?
256                CALL second_coppice_cut(npts, ipts, ivm, iele, ifm_cop, &
257                     icut_cop2, harvest_pool_bound, resolution, veget_max, biomass, &
258                     ind, circ_class_biomass, circ_class_n, circ_class_kill, &
259                     bm_to_litter, harvest_pool, harvest_type, harvest_cut, &
260                     harvest_area, circ_class_dist, last_cut, coppice_dens)
261
262                ! For short rotation coppices
263                CALL second_coppice_cut(npts, ipts, ivm, iele, ifm_src, &
264                     icut_cop2, harvest_pool_bound, resolution, veget_max, biomass, &
265                     ind, circ_class_biomass, circ_class_n, circ_class_kill, &
266                     bm_to_litter, harvest_pool, harvest_type, harvest_cut, &
267                     harvest_area, circ_class_dist, last_cut, coppice_dens)
268
269             ENDIF ! is_tree(ivm)
270
271             ! Here some circ classes may be empty.  We do not want to do
272             ! anything about that yet, since we have not killed the natural
273             ! mortality biomass.  If we redistribute the biomass now we will
274             ! change the circ_class_biomass, but the old values for this
275             ! were used to determine circ_class_kill(inatural) in
276             ! stomate_mark_kill.  If we change those values we will kill
277             ! the wrong amount of biomass in lpj_gap.
278
279          ENDDO ! loop over elements
280
281       ENDDO  ! loop over pfts
282
283    END DO ! loop over land points
284
285 !! 4. Check mass balance closure
286   
287    !  Consider this whole routine as a black box with incoming and outgoing
288    !  fluxes and a change in the mass of the box. Express in absolute
289    !  units gC (or gN), hence, multiply with dt and veget_max. In most routines
290    !  veget_max does not change and could be omitted but a general approach
291    !  was prefered.
292
293    !! 4.1 Calculate pools at the end of the routine
294    pool_end = zero
295    DO ipar = 1,nparts
296       DO iele = 1,nelements
297          pool_end(:,:,iele) = pool_end(:,:,iele) + &
298               (biomass(:,:,ipar,iele) * veget_max(:,:)) + &
299               (bm_to_litter(:,:,ipar,iele) * veget_max(:,:))
300       ENDDO
301    ENDDO
302   
303    ! The biomass harvest pool shouldn't be multiplied by veget_max
304    ! its units are already in gC pixel-1
305    ! Right now the harvest pool is in gC.  We will divide by the
306    ! resolution to make it per pixel.
307    !
308    DO ipts = 1, npts
309       pool_end(ipts,:,icarbon) = pool_end(ipts,:,icarbon) + &
310            SUM(harvest_pool(ipts,:,:,icarbon),2) / &
311            (resolution(ipts,1) * resolution(ipts,2))
312    ENDDO
313
314    !! 4.2 Calculate components of the mass balance
315    check_intern(:,:,iatm2land,icarbon) = zero
316    check_intern(:,:,iland2atm,icarbon) = -un * zero
317    check_intern(:,:,ilat2out,icarbon) = zero
318    check_intern(:,:,ilat2in,icarbon) = -un * zero
319    check_intern(:,:,ipoolchange,icarbon) = -un * (pool_end(:,:,icarbon) - &
320         pool_start(:,:,icarbon))
321    closure_intern = zero
322    DO imbc = 1,nmbcomp
323       closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + &
324            check_intern(:,:,imbc,icarbon)
325    ENDDO
326 
327    ! Write outcome
328    DO ipts=1,npts
329       DO ivm=1,nvm
330          IF(ABS(closure_intern(ipts,ivm,icarbon)) .LE. min_stomate)THEN
331             IF (ld_massbal) WRITE(numout,*) 'Mass balance closure in sapiens_kill.f90'
332          ELSE
333             WRITE(numout,*) 'Error: mass balance is not closed in sapiens_kill.f90'
334             WRITE(numout,*) '   ipts,ivm; ', ipts,ivm
335             WRITE(numout,*) '   Difference is, ', closure_intern(ipts,ivm,icarbon)
336             WRITE(numout,*) '   pool_end,pool_start: ', pool_end(ipts,ivm,icarbon), &
337                  pool_start(ipts,ivm,icarbon)
338             IF(ld_stop)THEN
339                CALL ipslerr_p (3,'gap_prog_sapiens', 'Mass balance error.','','')
340             ENDIF
341          ENDIF
342       ENDDO
343    ENDDO
344
345    IF (bavard.GE.2) WRITE(numout,*) 'Leaving gap_prog_sapiens'
346
347  END SUBROUTINE gap_prog_sapiens
348
349
350!! ================================================================================================================================
351!! SUBROUTINE   : thinning_harvest_aboveground_dead_roots
352!!
353!>\BRIEF        Transfer dead biomass to litter and update stand density due
354!!              to thinning operations
355!!
356!! DESCRIPTION  : Moves dead biomass to the litter pool and the harvest pool
357!!                where appropriate.  Since not all the trees are killed,
358!!                fewer counters have to be reset here than in the clearcut
359!!                case.  In addition, thinning will change the fraction of
360!!                trees killed by mortality, so we account for that as well.
361!!
362!!                NOTE: If harvest_fraction is 1, all the aboveground biomass is
363!!                      harvested.  If it's 1-branch_ratio, only stems are
364!!                      harvested.  Belowground wood, roots, fruits, and
365!!                      leaves all become litter.  A portion of aboveground
366!!                      wood is put in the harvest pool.  Reserves and labile
367!!                      carbon are assumed to be distributed evenly across
368!!                      all non-heartwood components.
369!!
370!! RECENT CHANGE(S):
371!!
372!! MAIN OUTPUT VARIABLE(S): ::circ_class_biomass; biomass, ::ind density of individuals,
373!! ::circ_class_kill
374!!
375!! REFERENCE(S)   :
376!!
377!! FLOWCHART    : None
378!!\n
379!_ ================================================================================================================================
380 
381  SUBROUTINE thinning_harvest_aboveground_dead_roots (npts, ipts, ivm, iele, &
382       ifm, icut, biomass, ind, bm_to_litter, &
383       circ_class_biomass, circ_class_kill, circ_class_n, harvest_pool, &
384       harvest_type, harvest_cut, harvest_area, &
385       harvest_pool_bound, resolution, veget_max, last_cut, harvest_fraction)
386
387
388    !! 0. Variable and parameter declaration
389
390    !! 0.1 Input variables
391    ! NOTE: if harvest_pool_bound is not explicitely defined
392    ! the model uses 1:ndia_harvest+2 as its dimensions. If
393    ! you like to do so remember to change the counters in
394    ! this routine
395    INTEGER(i_std), INTENT(in)                                   :: npts               !! Domain size (-)
396    INTEGER(i_std), INTENT(in)                                   :: ipts               !! current grid square
397    INTEGER(i_std), INTENT(in)                                   :: ivm                !! current PFT
398    INTEGER(i_std), INTENT(in)                                   :: iele               !! current element type
399    INTEGER(i_std), INTENT(in)                                   :: ifm                !! current fm strategy
400    INTEGER(i_std), INTENT(in)                                   :: icut               !! current cut time
401    REAL(r_std), DIMENSION(0:ndia_harvest+1), INTENT(in)         :: harvest_pool_bound !! The boundaries of the diameter classes
402                                                                                       !! in the wood harvest pools
403                                                                                       !! @tex $(m)$ @endtex
404    REAL(r_std),DIMENSION(npts,2), INTENT(in)                    :: resolution         !! The length of each grid-box in X and Y
405                                                                                       !! @tex $(m)$ @endtex
406    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: veget_max          !! Maximum fraction of vegetation type
407    REAL(r_std),INTENT(in)                                       :: harvest_fraction   !! The fraction of the aboveground
408                                                                                       !! biomass that is harvested.  This is
409                                                                                       !! generally 1-branch_ratio if you only
410                                                                                       !! want to harvest the stems, and 1 if
411                                                                                       !! you want to harvest everything.
412
413
414    !! 0.2 Output variables
415
416
417    !! 0.3 Modified variables   
418
419    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
420                                             INTENT(inout)       :: biomass            !! Biomass @tex $(gC m^{-2}) $@endtex
421    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)              :: ind                !! Stand level number of individuals
422                                                                                       !! @tex $(m^{-2})$ @endtex
423    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &                         
424                                              INTENT(inout)      :: circ_class_biomass !! Biomass of the components of the model 
425                                                                                       !! tree within a circumference
426                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
427    REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout)        :: circ_class_n       !! Number of individuals in each circ class
428                                                                                       !! @tex $(m^{-2})$ @endtex
429    REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times), &
430                                                  INTENT(inout)  :: circ_class_kill    !! Number of trees within a circ that needs
431                                                                                       !! to be killed @tex $(ind m^{-2})$ @endtex
432    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
433                                                  INTENT(inout)  :: bm_to_litter       !! Biomass transfer to litter
434                                                                                       !! @tex $(gC m^{-2})$ @endtex
435    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)               :: harvest_pool       !! The wood and biomass that have been
436                                                                                       !! havested by humans @tex $(gC)$ @endtex
437    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_type       !! Type of management that resulted
438                                                                                       !! in the harvest (unitless)
439    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_cut        !! Type of cutting that was used for the harvest
440                                                                                       !! (unitless)
441    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_area       !! Harvested area (m^{2})
442    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: last_cut           !! Years since last thinning or clearcut (years)
443
444    !! 0.4 Local variables
445
446    INTEGER(i_std)                                               :: icir               !! Indices
447    REAL(r_std)                                                  :: total_bm
448    REAL(r_std)                                                  :: thinned_bm
449    REAL(r_std)                                                  :: dead_bm
450    REAL(r_std)                                                  :: bm_difference
451
452    !_ ================================================================================================================================
453
454    IF (bavard.GE.2) WRITE(numout,*) &
455         'Entering thinning_harvest_aboveground_dead_roots sapiens_kill.f90'
456
457    !  Initialize check for surface area conservation
458    !  Veget_max is a INTENT(in) variable and can therefore
459    !  not be changed during the course of this subroutine
460    !  No need to check whether the subroutine preserves the
461    !  total surface area of the pixel.
462
463    IF( SUM(circ_class_kill(ipts,ivm,:,ifm,icut)) .GT. min_stomate ) THEN
464
465       CALL harvest_aboveground_dead_roots(npts, ipts, ivm, iele, &
466            ifm, icut, biomass, ind, bm_to_litter, &
467            circ_class_biomass, circ_class_kill, circ_class_n, harvest_pool, &
468            harvest_type, harvest_cut, harvest_area, &
469            harvest_pool_bound, resolution, veget_max, harvest_fraction)
470
471       ! thinning trees generally avoids natural mortality and self-thinning, so
472       ! let's remove all thinned trees from the pools.  This is tricky because
473       ! mortality depends on a fraction of the total biomass, which may be spread
474       ! out over different circ classes than what were killed by thinning.
475       thinned_bm = zero
476       dead_bm = zero
477       DO icir=1,ncirc
478          total_bm=SUM(circ_class_biomass(ipts,ivm,icir,:,iele))
479          thinned_bm=thinned_bm+circ_class_kill(ipts,ivm,icir,ifm,icut)*total_bm
480          dead_bm=dead_bm+circ_class_kill(ipts,ivm,icir,ifm_none,icut_thin)*total_bm
481       ENDDO
482       bm_difference = dead_bm - thinned_bm
483
484       ! if the thinned biomass is more than what would have died from natural causes,
485       ! we can just assume that thinning killed all these trees.
486       IF(bm_difference .LE. min_stomate)THEN
487
488          circ_class_kill(ipts,ivm,:,ifm_none,icut_thin)=zero
489
490       ELSE
491
492          ! if not, we have to reduce the amount of biomass that is killed naturally
493          ! by what was removed by thinning, redistributing the remainder properly
494          ! across all circ classes as we did in stomate_mark_kill.
495          circ_class_kill(ipts,ivm,:,ifm_none,icut_thin) = zero
496          CALL distribute_mortality_biomass( bm_difference, &
497               death_distribution_factor(ivm), &
498               circ_class_n(ipts,ivm,:), &
499               circ_class_biomass(ipts,ivm,:,:,icarbon), &
500               circ_class_kill(ipts,ivm,:,ifm_none,icut_thin) )
501       ENDIF ! bm_difference .LE. min_stomate
502
503
504       ! make sure nothing has gone negative
505       WHERE ( circ_class_kill(ipts,ivm,:,ifm_none,icut_thin) .LT. zero )
506          circ_class_kill(ipts,ivm,:,ifm_none,icut_thin)=zero
507       END WHERE
508
509       ! now zero out this circ_class_kill, since we should never try
510       ! to kill these same trees again.
511       circ_class_kill(ipts,ivm,:,ifm,icut)=zero
512
513       ! Reset the values for the time since the last operation
514       last_cut(ipts,ivm)=0
515       
516    ENDIF ! checking to see if we have biomass to thin
517
518
519    IF (bavard.GE.2) WRITE(numout,*) &
520         'Leaving thinning_harvest_aboveground_dead_roots sapiens_kill.f90'
521
522  END SUBROUTINE thinning_harvest_aboveground_dead_roots
523
524!! ================================================================================================================================
525!! SUBROUTINE   : harvest_aboveground_dead_roots
526!!
527!>\BRIEF        Transfer dead biomass to litter and wood to harvest pool.
528!!
529!! DESCRIPTION  : Moves dead biomass to the litter pool and the harvest pool
530!!                where appropriate. 
531!!
532!!                NOTE: If harvest_fraction is 1, all the aboveground biomass is
533!!                      harvested.  If it's 1-branch_ratio, only stems are
534!!                      harvested.  Belowground wood, roots, fruits, and
535!!                      leaves all become litter.  A portion of aboveground
536!!                      wood is put in the harvest pool.  Reserves and labile
537!!                      carbon are assumed to be distributed evenly across
538!!                      all non-heartwood components.
539!!
540!! RECENT CHANGE(S):
541!!
542!! MAIN OUTPUT VARIABLE(S): ::bm_to_litter, ::harvest_pool
543!!
544!! REFERENCE(S)   :
545!!
546!! FLOWCHART    : None
547!!\n
548!_ ================================================================================================================================
549 
550  SUBROUTINE harvest_aboveground_dead_roots (npts, ipts, ivm, iele, &
551       ifm, icut, biomass, ind, bm_to_litter, &
552       circ_class_biomass, circ_class_kill, circ_class_n, harvest_pool, &
553       harvest_type, harvest_cut, harvest_area, &
554       harvest_pool_bound, resolution, veget_max, harvest_fraction)
555
556
557    !! 0. Variable and parameter declaration
558
559    !! 0.1 Input variables
560    ! NOTE: if harvest_pool_bound is not explicitely defined
561    ! the model uses 1:ndia_harvest+2 as its dimensions. If
562    ! you like to do so remember to change the counters in
563    ! this routine
564    INTEGER(i_std), INTENT(in)                                   :: npts               !! Domain size (-)
565    INTEGER(i_std), INTENT(in)                                   :: ipts               !! current grid square
566    INTEGER(i_std), INTENT(in)                                   :: ivm                !! current PFT
567    INTEGER(i_std), INTENT(in)                                   :: iele               !! current element type
568    INTEGER(i_std), INTENT(in)                                   :: ifm                !! current fm strategy
569    INTEGER(i_std), INTENT(in)                                   :: icut               !! current cut time
570    REAL(r_std), DIMENSION(0:ndia_harvest+1), INTENT(in)         :: harvest_pool_bound !! The boundaries of the diameter classes
571                                                                                       !! in the wood harvest pools
572                                                                                       !! @tex $(m)$ @endtex
573    REAL(r_std),DIMENSION(npts,2), INTENT(in)                    :: resolution         !! The length of each grid-box in X and Y
574                                                                                       !! @tex $(m)$ @endtex
575    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: veget_max          !! Maximum fraction of vegetation type
576    REAL(r_std),INTENT(in)                                       :: harvest_fraction   !! The fraction of the aboveground
577                                                                                       !! biomass that is harvested.  This is
578                                                                                       !! generally 1-branch_ratio if you only
579                                                                                       !! want to harvest the stems, and 1 if
580                                                                                       !! you want to harvest everything.
581
582    !! 0.2 Output variables
583
584
585    !! 0.3 Modified variables   
586
587    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
588                                             INTENT(inout)       :: biomass            !! Biomass @tex $(gC m^{-2}) $@endtex
589    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)              :: ind                !! Stand level number of individuals
590                                                                                       !! @tex $(m^{-2})$ @endtex
591    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &                         
592                                              INTENT(inout)      :: circ_class_biomass !! Biomass of the components of the model 
593                                                                                       !! tree within a circumference
594                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
595    REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout)        :: circ_class_n       !! Number of individuals in each circ class
596                                                                                       !! @tex $(m^{-2})$ @endtex
597    REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times), &
598                                                  INTENT(inout)  :: circ_class_kill    !! Number of trees within a circ that needs
599                                                                                       !! to be killed @tex $(ind m^{-2})$ @endtex
600    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
601                                                  INTENT(inout)  :: bm_to_litter       !! Biomass transfer to litter
602                                                                                       !! @tex $(gC m^{-2})$ @endtex
603    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)               :: harvest_pool       !! The wood and biomass that have been
604                                                                                       !! havested by humans @tex $(gC)$ @endtex
605    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_type       !! Type of management that resulted
606                                                                                       !! in the harvest (unitless)
607    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_cut        !! Type of cutting that was used for the harvest
608                                                                                       !! (unitless)
609    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_area       !! Harvested area (m^{2})
610
611
612    !! 0.4 Local variables
613
614    INTEGER(i_std)                                               :: idia, icir         !! Indices
615    REAL(r_std)                                                  :: total_nonheart_biomass_harvest
616    REAL(r_std)                                                  :: total_nonheart_biomass_kill
617    REAL(r_std)                                                  :: harvest_ratio
618    REAL(r_std), DIMENSION(ncirc)                                :: diameters_temp
619
620    !_ ================================================================================================================================
621
622    IF (bavard.GE.2) WRITE(numout,*) &
623         'Entering harvest_aboveground_dead_roots sapiens_kill.f90'
624
625    !  Initialize check for surface area conservation
626    !  Veget_max is a INTENT(in) variable and can therefore
627    !  not be changed during the course of this subroutine
628    !  No need to check whether the subroutine preserves the
629    !  total surface area of the pixel.
630   
631    ! we will need the trunk diameters for all the circ classes in order to figure
632    ! out which harvest pool to put them in.
633    diameters_temp(:)=wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm) 
634
635    DO icir=1,ncirc
636
637       ! It's possible that we only kill individuals in certain circ classes.
638       ! In this case, we need to check for each circ class to make sure
639       ! we don't divide by zero.
640       IF(circ_class_kill(ipts,ivm,icir,ifm,icut) .LT. min_stomate) CYCLE
641
642       IF(ld_kill .AND. ipts == test_grid .AND. ivm == test_pft)THEN
643          WRITE(numout,*) 'Killing trees in harvest_aboveground_dead_roots'
644          WRITE(numout,*) 'ipts,ivm,icir,ifm,icut: ',ipts,ivm,icir,ifm,icut
645          WRITE(numout,*) 'circ_class_kill,circ_class_n: ',&
646               circ_class_kill(ipts,ivm,icir,ifm,icut),circ_class_n(ipts,ivm,icir)
647       ENDIF
648
649       ! The decision of what to do with the carboyhydrate reserves and the
650       ! labile carbon pool are tricky.  In turnover, they are not touched
651       ! when leaves/roots/fruits die, which implies that they are stored
652       ! primarily in the above and belowground sapwood (heartwood is dead).
653       ! However, here a tree does not know thinning is coming, so it cannot
654       ! move this carbon around.  We will therefore determine where is goes
655       ! based on the total amount of non-heartwood pools and where they go.
656       total_nonheart_biomass_kill=circ_class_kill(ipts,ivm,icir,ifm,icut)*&
657            (circ_class_biomass(ipts,ivm,icir,isapbelow,iele)+&
658            circ_class_biomass(ipts,ivm,icir,iroot,iele)+&
659            circ_class_biomass(ipts,ivm,icir,ileaf,iele)+&
660            circ_class_biomass(ipts,ivm,icir,ifruit,iele)+&
661            circ_class_biomass(ipts,ivm,icir,isapabove,iele))
662
663       total_nonheart_biomass_harvest=harvest_fraction*&
664            circ_class_kill(ipts,ivm,icir,ifm,icut)*&
665            circ_class_biomass(ipts,ivm,icir,isapabove,iele)
666
667       ! This is now the fraction of the nonheartwood biomass that will go
668       ! into the harvest pool.  We will use this to partition the carbohydrate
669       ! reserve and labile pools into litter and harvest.
670       harvest_ratio=total_nonheart_biomass_harvest/total_nonheart_biomass_kill
671
672       ! all of the following biomass becomes litter
673       bm_to_litter(ipts,ivm,isapbelow,iele) = &
674            bm_to_litter(ipts,ivm,isapbelow,iele) + &
675            circ_class_biomass(ipts,ivm,icir,isapbelow,iele)*&
676            circ_class_kill(ipts,ivm,icir,ifm,icut)
677       bm_to_litter(ipts,ivm,iheartbelow,iele) = &
678            bm_to_litter(ipts,ivm,iheartbelow,iele) + &
679            circ_class_biomass(ipts,ivm,icir,iheartbelow,iele)*&
680            circ_class_kill(ipts,ivm,icir,ifm,icut)
681       bm_to_litter(ipts,ivm,iroot,iele) = bm_to_litter(ipts,ivm,iroot,iele) + &
682            circ_class_biomass(ipts,ivm,icir,iroot,iele)*&
683            circ_class_kill(ipts,ivm,icir,ifm,icut)
684       bm_to_litter(ipts,ivm,ileaf,iele) = bm_to_litter(ipts,ivm,ileaf,iele) + &
685            circ_class_biomass(ipts,ivm,icir,ileaf,iele)*&
686            circ_class_kill(ipts,ivm,icir,ifm,icut)
687       bm_to_litter(ipts,ivm,ifruit,iele) = bm_to_litter(ipts,ivm,ifruit,iele) + &
688            circ_class_biomass(ipts,ivm,icir,ifruit,iele)*&
689            circ_class_kill(ipts,ivm,icir,ifm,icut)
690
691       ! Include the branches as litter, if necessary.
692       bm_to_litter(ipts,ivm,isapabove,iele) = &
693            bm_to_litter(ipts,ivm,isapabove,iele) + &
694            (un - harvest_fraction)*circ_class_biomass(ipts,ivm,icir,isapabove,iele)*&
695            circ_class_kill(ipts,ivm,icir,ifm,icut)
696       bm_to_litter(ipts,ivm,iheartabove,iele) = &
697            bm_to_litter(ipts,ivm,iheartabove,iele) + &
698            (un - harvest_fraction)*circ_class_biomass(ipts,ivm,icir,iheartabove,iele)*&
699            circ_class_kill(ipts,ivm,icir,ifm,icut)
700
701       ! Don't forget to include the proper proportion of the labile and
702       ! carbres pools.
703       bm_to_litter(ipts,ivm,icarbres,iele) = &
704            bm_to_litter(ipts,ivm,icarbres,iele) + &
705            (un - harvest_ratio )*circ_class_biomass(ipts,ivm,icir,icarbres,iele)*&
706            circ_class_kill(ipts,ivm,icir,ifm,icut)
707       bm_to_litter(ipts,ivm,ilabile,iele) = &
708            bm_to_litter(ipts,ivm,ilabile,iele) + &
709            (un - harvest_ratio )*circ_class_biomass(ipts,ivm,icir,ilabile,iele)*&
710            circ_class_kill(ipts,ivm,icir,ifm,icut)
711
712       ! now we have the stems.  These are what we want to save, and we are
713       ! going to classify them by PFT type, location, and diameter so that
714       ! we can compare to FAO statistics and do some more interesting
715       ! post-treatment.
716
717       ! UNITS: The units of the biomass pools are per meter squared.
718       ! It seems to make more sense to store the harvest pool in absolute
719       ! units, which means we have to multiple the result by the fraction of
720       ! this grid cell covered by this PFT (veget_max) and the size of the
721       ! grid cell.  It seems that the fraction of nobio land does not need
722       ! to be used here, since frac_nobio and veget_max are summed together
723       ! in slowproc and they must equal one.  This gives
724       ! us the grams of carbon harvest, which is a big number but that's
725       ! why people invented scientific notation.
726
727       dia_loop:DO idia=1,ndia_harvest
728
729          ! Move the harvest to the correct diameter class and
730          ! account for the largest diameter class because
731          ! harvest_pool_bound has only ndia_harvest+1 dimensions
732          IF ( (diameters_temp(icir) .LE. harvest_pool_bound(idia)) .OR. &
733               ((diameters_temp(icir) .GT. harvest_pool_bound(ndia_harvest)) .AND. &
734               (diameters_temp(icir) .LT. harvest_pool_bound(ndia_harvest+1))) ) THEN
735
736             ! this is the correct diameter class
737             harvest_pool(ipts,ivm,idia,icarbon) = &
738                  harvest_pool(ipts,ivm,idia,icarbon) + &
739                  harvest_fraction*circ_class_kill(ipts,ivm,icir,ifm,icut)*&
740                  (circ_class_biomass(ipts,ivm,icir,isapabove,iele)+&
741                  circ_class_biomass(ipts,ivm,icir,iheartabove,iele)) * &
742                  veget_max(ipts,ivm) * &
743                  resolution(ipts,1) * resolution(ipts,2)
744
745             ! also add the proper fractions of the labile and carbres pools
746             harvest_pool(ipts,ivm,idia,icarbon) = &
747                  harvest_pool(ipts,ivm,idia,icarbon)+ &
748                  harvest_ratio * circ_class_kill(ipts,ivm,icir,ifm,icut)*&
749                  ( circ_class_biomass(ipts,ivm,icir,ilabile,iele) + &
750                  circ_class_biomass(ipts,ivm,icir,icarbres,iele) )* &
751                  veget_max(ipts,ivm) * &
752                  resolution(ipts,1) * resolution(ipts,2)
753             
754             ! Associated harvest variables
755             harvest_type(ipts,ivm) = ifm
756             harvest_cut(ipts,ivm) = icut
757             harvest_area(ipts,ivm) = veget_max(ipts,ivm) * &
758                  resolution(ipts,1) * resolution(ipts,2)
759               
760             IF((ld_kill .OR. ld_forestry) .AND. &
761                  ipts == test_grid .AND. ivm == test_pft)THEN
762                WRITE(numout,*) 'Harvesting harvest_aboveground_dead_roots 1 ',ipts,ivm,icir,idia
763                WRITE(numout,*) 'harvest_area,harvest_type,harvest_cut: ',&
764                     harvest_area(ipts,ivm),harvest_type(ipts,ivm),harvest_cut(ipts,ivm)
765                WRITE(numout,*) 'Woody biomass: ',&
766                     harvest_fraction*circ_class_kill(ipts,ivm,icir,ifm,icut)*&
767                  (circ_class_biomass(ipts,ivm,icir,isapabove,iele)+&
768                  circ_class_biomass(ipts,ivm,icir,iheartabove,iele))
769                WRITE(numout,*) 'Reserves: ',&
770                     harvest_ratio * circ_class_kill(ipts,ivm,icir,ifm,icut)*&
771                  ( circ_class_biomass(ipts,ivm,icir,ilabile,iele) + &
772                  circ_class_biomass(ipts,ivm,icir,icarbres,iele) )
773                WRITE(numout,*) 'cckill,harvest_pool ',&
774                     circ_class_kill(ipts,ivm,icir,ifm,icut),&
775                     harvest_pool(ipts,ivm,idia,icarbon)
776             ENDIF
777
778             EXIT dia_loop
779
780          ENDIF ! diameters_temp(icir) .LE. harvest_pool_bound(idia)
781
782          ! If we are here, it means the diameter of this trunk is larger
783          ! than the upper boundary of our harvest bins.  Considering the
784          ! upper boundary is set to be 99999 meters in constants.f90,
785          ! something has gone wrong!
786          IF (diameters_temp(icir) .GE. harvest_pool_bound(ndia_harvest+1)) THEN
787
788             WRITE(numout,*) 'ERROR: Our tree is too big to fit into the harvest pools.'
789             WRITE(numout,*) 'ipts, ivm, icir: ',ipts, ivm, icir
790             WRITE(numout,*) 'diameters_temp(icir), harvest_pool_bound(idia)',&
791                  diameters_temp(icir),harvest_pool_bound(ndia_harvest+1)
792             CALL ipslerr_p (3,'harvest_aboveground_dead_roots', &
793                  'Tree too big to fit into the harvest pool.',&
794                  'Look into the output file for details.', &
795                  '')
796          ENDIF
797
798       ENDDO dia_loop
799
800       ! remove the number of individuals that died
801       circ_class_n(ipts,ivm,icir) = circ_class_n(ipts,ivm,icir) - &
802            circ_class_kill(ipts,ivm,icir,ifm,icut)
803
804       ! remove the dead biomass
805       biomass(ipts,ivm,:,iele) = biomass(ipts,ivm,:,iele) - &
806            circ_class_biomass(ipts,ivm,icir,:,iele)*&
807            circ_class_kill(ipts,ivm,icir,ifm,icut)
808
809       ! adjust the number of individuals
810       ind(ipts,ivm)=ind(ipts,ivm)-circ_class_kill(ipts,ivm,icir,ifm,icut)
811
812    ENDDO ! loop over circ classes
813
814    ! We should not have negative numbers.  I don't want to zero out all
815    ! of them, because it could be sign of a real problem, but if they
816    ! are a very small negative number it's just a numerical issue.
817    WHERE( ABS(biomass(ipts,ivm,:,iele)) .LE. min_stomate)
818       biomass(ipts,ivm,:,iele) = zero
819    END WHERE
820    WHERE( ABS(circ_class_biomass(ipts,ivm,:,:,iele)) .LE. min_stomate)
821      circ_class_biomass(ipts,ivm,:,:,iele) = zero
822    END WHERE
823    WHERE( ABS(circ_class_n(ipts,ivm,:)) .LE. min_stomate)
824       circ_class_n(ipts,ivm,:) = zero
825    END WHERE
826
827    IF (bavard.GE.2) WRITE(numout,*) &
828         'Leaving harvest_aboveground_dead_roots sapiens_kill.f90'
829
830  END SUBROUTINE harvest_aboveground_dead_roots
831
832!! ================================================================================================================================
833!! SUBROUTINE   : harvest_aboveground_wood_living_roots
834!!
835!>\BRIEF        Transfer dead biomass to litter and wood to harvest pool.
836!!
837!! DESCRIPTION  : Moves dead biomass to the litter pool and the harvest pool
838!!                where appropriate. 
839!!
840!!                NOTE: All branches and stems are harvested.  Basically, this means
841!!                      all aboveground woody biomass.  The fine roots
842!!                      are killed and left to litter.  The sapwood and heartwood
843!!                      belowground is left alive.  The leaves are killed, and the
844!!                      labile and carbres pools are left untouched in the living
845!!                      biomass.
846!!
847!! RECENT CHANGE(S):
848!!
849!! MAIN OUTPUT VARIABLE(S): ::bm_to_litter, ::harvest_pool
850!!
851!! REFERENCE(S)   :
852!! - Bellassen, V., Le Maire, G., Dhote, J.F., Viovy, N., Ciais, P., 2010.
853!! Modeling forest management within a global vegetation model – Part 1:
854!! model structure and general behaviour. Ecological Modelling 221, 2458–2474.
855!!
856!! FLOWCHART    : None
857!!\n
858!_ ================================================================================================================================
859 
860  SUBROUTINE harvest_aboveground_wood_living_roots (npts, ipts, ivm, iele, &
861       ifm, icut, biomass, bm_to_litter, &
862       circ_class_biomass, circ_class_kill, harvest_pool, &
863       harvest_type, harvest_cut, harvest_area, &
864       harvest_pool_bound, resolution, veget_max)
865
866
867    !! 0. Variable and parameter declaration
868
869    !! 0.1 Input variables
870    ! NOTE: if harvest_pool_bound is not explicitely defined
871    ! the model uses 1:ndia_harvest+2 as its dimensions. If
872    ! you like to do so remember to change the counters in
873    ! this routine
874    INTEGER(i_std), INTENT(in)                                   :: npts               !! Domain size (-)
875    INTEGER(i_std), INTENT(in)                                   :: ipts               !! current grid square
876    INTEGER(i_std), INTENT(in)                                   :: ivm                !! current PFT
877    INTEGER(i_std), INTENT(in)                                   :: iele               !! current element type
878    INTEGER(i_std), INTENT(in)                                   :: ifm                !! current fm strategy
879    INTEGER(i_std), INTENT(in)                                   :: icut               !! current cut time
880    REAL(r_std), DIMENSION(0:ndia_harvest+1), INTENT(in)         :: harvest_pool_bound !! The boundaries of the diameter classes
881                                                                                       !! in the wood harvest pools
882                                                                                       !! @tex $(m)$ @endtex
883    REAL(r_std),DIMENSION(npts,2), INTENT(in)                    :: resolution         !! The length of each grid-box in X and Y
884                                                                                       !! @tex $(m)$ @endtex
885    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: veget_max          !! Maximum fraction of vegetation type
886
887
888    !! 0.2 Output variables
889
890
891    !! 0.3 Modified variables   
892
893    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
894                                             INTENT(inout)       :: biomass            !! Biomass @tex $(gC m^{-2}) $@endtex
895    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &                         
896                                              INTENT(inout)      :: circ_class_biomass !! Biomass of the components of the model 
897                                                                                       !! tree within a circumference
898                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
899    REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times), &
900                                                  INTENT(inout)  :: circ_class_kill    !! Number of trees within a circ that needs
901                                                                                       !! to be killed @tex $(ind m^{-2})$ @endtex
902    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
903                                                  INTENT(inout)  :: bm_to_litter       !! Biomass transfer to litter
904                                                                                       !! @tex $(gC m^{-2})$ @endtex
905    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)               :: harvest_pool       !! The wood and biomass that have been
906                                                                                       !! havested by humans @tex $(gC)$ @endtex
907    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_type       !! Type of management that resulted
908                                                                                       !! in the harvest (unitless)
909    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_cut        !! Type of cutting that was used for the harvest
910                                                                                       !! (unitless)
911    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_area       !! Harvested area (m^{2})
912
913    !! 0.4 Local variables
914
915    INTEGER(i_std)                                               :: idia, icir         !! Indices
916    REAL(r_std)                                                  :: total_nonheart_biomass_harvest
917    REAL(r_std)                                                  :: total_nonheart_biomass_living
918    REAL(r_std)                                                  :: total_nonheart_biomass_kill
919    REAL(r_std)                                                  :: harvest_ratio
920    REAL(r_std)                                                  :: living_ratio
921    REAL(r_std)                                                  :: litter_ratio
922    REAL(r_std), DIMENSION(ncirc)                                :: diameters_temp
923    REAL(r_std)                                                  :: dead_heartwood_be
924    REAL(r_std)                                                  :: dead_sapwood_be
925
926    !_ ================================================================================================================================
927
928    IF (bavard.GE.2) WRITE(numout,*) &
929         'Entering harvest_aboveground_wood_living_roots sapiens_kill.f90'
930
931    !  Initialize check for surface area conservation
932    !  Veget_max is a INTENT(in) variable and can therefore
933    !  not be changed during the course of this subroutine
934    !  No need to check whether the subroutine preserves the
935    !  total surface area of the pixel.
936
937    ! we will need the trunk diameters for all the circ classes in order to figure
938    ! out which harvest pool to put them in.
939    diameters_temp(:)=wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm) 
940
941    DO icir=1,ncirc
942
943       ! It's possible that we only kill individuals in certain circ classes.
944       ! In this case, we need to check for each circ class to make sure
945       ! we don't divide by zero.
946       IF(circ_class_kill(ipts,ivm,icir,ifm,icut) .LT. min_stomate) CYCLE
947
948       IF(ld_kill .AND. ipts == test_grid .AND. ivm == test_pft)THEN
949          WRITE(numout,*) 'Killing trees in harvest_aboveground_wood_living_roots'
950          WRITE(numout,*) 'ipts,ivm,icir,ifm,icut: ',ipts,ivm,icir,ifm,icut
951          WRITE(numout,*) 'circ_class_kill: ',&
952               circ_class_kill(ipts,ivm,icir,ifm,icut)
953       ENDIF
954
955       ! The decision of what to do with the carboyhydrate reserves and the
956       ! labile carbon pool are tricky.  In turnover, they are not touched
957       ! when leaves/roots/fruits die, which implies that they are stored
958       ! primarily in the above and belowground sapwood (heartwood is dead).
959       ! However, here a tree does not know thinning is coming, so it cannot
960       ! move this carbon around.  In theory, we should remove part of the
961       ! labile and reserve pools to the different product pools, but
962       ! these pools are very sensitive for the allocation.  Therefore
963       ! we don't touch them, which is equivalent to saying the reserve
964       ! and labile carbon are all found in the roots.
965
966       total_nonheart_biomass_kill=circ_class_kill(ipts,ivm,icir,ifm,icut)*&
967            (circ_class_biomass(ipts,ivm,icir,isapbelow,iele)+&
968            circ_class_biomass(ipts,ivm,icir,iroot,iele)+&
969            circ_class_biomass(ipts,ivm,icir,ileaf,iele)+&
970            circ_class_biomass(ipts,ivm,icir,ifruit,iele)+&
971            circ_class_biomass(ipts,ivm,icir,isapabove,iele))
972
973       total_nonheart_biomass_harvest=&
974            circ_class_kill(ipts,ivm,icir,ifm,icut)*&
975            (circ_class_biomass(ipts,ivm,icir,isapabove,iele))
976
977       total_nonheart_biomass_living=&
978            circ_class_kill(ipts,ivm,icir,ifm,icut)*&
979            (circ_class_biomass(ipts,ivm,icir,isapbelow,iele))
980
981       ! This is now the fraction of the nonheartwood biomass that will go
982       ! into the harvest pool.  We will use this to partition the carbohydrate
983       ! reserve and labile pools into litter and harvest.
984       harvest_ratio=total_nonheart_biomass_harvest/total_nonheart_biomass_kill
985       living_ratio=total_nonheart_biomass_living/total_nonheart_biomass_kill
986       litter_ratio=un-harvest_ratio-living_ratio
987       IF(litter_ratio .LT. -min_stomate)THEN
988          WRITE(numout,*) 'We have a problem with coppicing!'
989          WRITE(numout,*) 'ipts,ivm,icir,ifm,icut ',ipts,ivm,icir,ifm,icut
990          WRITE(numout,*) 'harvest_ratio, living_ratio, litter_ratio ',&
991               harvest_ratio, living_ratio, litter_ratio
992          CALL ipslerr_p (3,'harvest_aboveground_wood_living_roots', &
993               'litter_ratio is negative.',&
994               'Look into the output file for details.', &
995               '')
996       ELSEIF(litter_ratio .LT. -min_stomate)THEN
997          litter_ratio=zero
998       ENDIF
999
1000       ! all of the following biomass becomes litter
1001       bm_to_litter(ipts,ivm,iroot,iele) = bm_to_litter(ipts,ivm,iroot,iele) + &
1002            circ_class_biomass(ipts,ivm,icir,iroot,iele)*&
1003            circ_class_kill(ipts,ivm,icir,ifm,icut)
1004       bm_to_litter(ipts,ivm,ileaf,iele) = bm_to_litter(ipts,ivm,ileaf,iele) + &
1005            circ_class_biomass(ipts,ivm,icir,ileaf,iele)*&
1006            circ_class_kill(ipts,ivm,icir,ifm,icut)
1007       bm_to_litter(ipts,ivm,ifruit,iele) = bm_to_litter(ipts,ivm,ifruit,iele) + &
1008            circ_class_biomass(ipts,ivm,icir,ifruit,iele)*&
1009            circ_class_kill(ipts,ivm,icir,ifm,icut)
1010
1011       ! We also turn part of the belowground wood into litter.  The reason for this
1012       ! is that some coppiced forests were not being able to regenerate due to
1013       ! the amount of belowground wood, which kept the trees out of balance for
1014       ! the whole next year.  The leaves grew, but never enough to put wood into
1015       ! the stems, and the leaves died at the end of the year.  Then everything
1016       ! started over exactly the same the following years.  The biomass stayed
1017       ! constant (the belowground mass was two orders of magnitude larger than any
1018       ! leaf mass which grew) and it never died.  Killing part of the belowground
1019       ! mass should allow the tree to allometrically allocate faster.  At the
1020       ! same time, it's not entirely unrealistic.  Parts of the roots stay alive
1021       ! during coppicing, but do they all?  Unknown.  This is probably only an issue
1022       ! for full coppicing, not for short rotation coppices, since those stands
1023       ! are much younger.
1024       dead_sapwood_be=circ_class_biomass(ipts,ivm,icir,isapbelow,iele)&
1025            *coppice_kill_be_wood(ivm)
1026       bm_to_litter(ipts,ivm,isapbelow,iele) = bm_to_litter(ipts,ivm,isapbelow,iele) + &
1027            dead_sapwood_be*&
1028            circ_class_kill(ipts,ivm,icir,ifm,icut)
1029       dead_heartwood_be=circ_class_biomass(ipts,ivm,icir,iheartbelow,iele)&
1030            *coppice_kill_be_wood(ivm)
1031       bm_to_litter(ipts,ivm,iheartbelow,iele) = bm_to_litter(ipts,ivm,iheartbelow,iele) + &
1032            dead_heartwood_be*&
1033            circ_class_kill(ipts,ivm,icir,ifm,icut)
1034
1035       ! now we have the stems.  These are what we want to save, and we are
1036       ! going to classify them by PFT type, location, and diameter so that
1037       ! we can compare to FAO statistics and do some more interesting
1038       ! post-treatment.
1039
1040       ! UNITS: The units of the biomass pools are per meter squared.
1041       ! It seems to make more sense to store the harvest pool in absolute
1042       ! units, which means we have to multiple the result by the fraction of
1043       ! this grid cell covered by this PFT (veget_max) and the size of the
1044       ! grid cell.  It seems that the fraction of nobio land does not need
1045       ! to be used here, since frac_nobio and veget_max are summed together
1046       ! in slowproc and they must equal one.  This gives
1047       ! us the grams of carbon harvest, which is a big number but that's
1048       ! why people invented scientific notation.
1049
1050       dia_loop:DO idia=1,ndia_harvest
1051
1052          ! Move the harvest to the correct diameter class and
1053          ! account for the largest diameter class because
1054          ! harvest_pool_bound has only ndia_harvest+1 dimensions
1055          IF ( (diameters_temp(icir) .LE. harvest_pool_bound(idia)) .OR. &
1056               ((diameters_temp(icir) .GT. harvest_pool_bound(ndia_harvest)) .AND. &
1057               (diameters_temp(icir) .LT. harvest_pool_bound(ndia_harvest+1))) ) THEN
1058
1059             ! this is the correct diameter class
1060             harvest_pool(ipts,ivm,idia,icarbon)=&
1061                  harvest_pool(ipts,ivm,idia,icarbon)+ &
1062                  circ_class_kill(ipts,ivm,icir,ifm,icut)*&
1063                  (circ_class_biomass(ipts,ivm,icir,isapabove,iele)+&
1064                  circ_class_biomass(ipts,ivm,icir,iheartabove,iele)) * &
1065                  veget_max(ipts,ivm) * &
1066                  resolution(ipts,1) * resolution(ipts,2)
1067             
1068             ! Associated harvest variables
1069             harvest_type(ipts,ivm) = ifm
1070             harvest_cut(ipts,ivm) = icut
1071             harvest_area(ipts,ivm) = veget_max(ipts,ivm) * &
1072                  resolution(ipts,1) * resolution(ipts,2)
1073
1074             IF((ld_kill .OR. ld_forestry) .AND. &
1075                  ipts == test_grid .AND. ivm == test_pft)THEN
1076                WRITE(numout,*) 'Harvesting harvest_aboveground_living_roots 1 ',ipts,ivm,icir
1077                WRITE(numout,*) 'harvest_area,harvest_type,harvest_cut: ',&
1078                     harvest_area(ipts,ivm),harvest_type(ipts,ivm),harvest_cut(ipts,ivm)
1079                WRITE(numout,*) 'Woody biomass: ',&
1080                     circ_class_kill(ipts,ivm,icir,ifm,icut)*&
1081                  (circ_class_biomass(ipts,ivm,icir,isapabove,iele)+&
1082                  circ_class_biomass(ipts,ivm,icir,iheartabove,iele))
1083             ENDIF
1084
1085             EXIT dia_loop
1086
1087          ENDIF ! diameters_temp(icir) .LE. harvest_pool_bound(idia)
1088
1089          ! If we are here, it means the diameter of this trunk is larger
1090          ! than the upper boundary of our harvest bins.  Considering the
1091          ! upper boundary is set to be 99999 meters in constants.f90,
1092          ! something has gone wrong!
1093          IF (diameters_temp(icir) .GE. harvest_pool_bound(ndia_harvest+1)) THEN
1094
1095             WRITE(numout,*) 'Our tree is too big to fit into the harvest pools.'
1096             WRITE(numout,*) 'ipts, ivm, icir: ',ipts, ivm, icir
1097             WRITE(numout,*) 'diameters_temp(icir), harvest_pool_bound(idia)',&
1098                  diameters_temp(icir),harvest_pool_bound(ndia_harvest+1)
1099             CALL ipslerr_p (3,'harvest_aboveground_wood_living_roots', &
1100                  'Tree is too big for the harvest pool.',&
1101                  'Look into the output file for details.', &
1102                  '')
1103
1104          ENDIF
1105
1106       ENDDO dia_loop
1107
1108
1109       ! No individuals die here since the roots are still alive.
1110
1111       ! Remove the dead biomass.  This includes everything except
1112       ! the belowground biomass. 
1113       biomass(ipts,ivm,ileaf,iele) = biomass(ipts,ivm,ileaf,iele) - &
1114            circ_class_biomass(ipts,ivm,icir,ileaf,iele)*&
1115            circ_class_kill(ipts,ivm,icir,ifm,icut)
1116       biomass(ipts,ivm,isapabove,iele) = biomass(ipts,ivm,isapabove,iele) - &
1117            circ_class_biomass(ipts,ivm,icir,isapabove,iele)*&
1118            circ_class_kill(ipts,ivm,icir,ifm,icut)
1119       biomass(ipts,ivm,iheartabove,iele) = biomass(ipts,ivm,iheartabove,iele) - &
1120            circ_class_biomass(ipts,ivm,icir,iheartabove,iele)*&
1121            circ_class_kill(ipts,ivm,icir,ifm,icut)
1122       biomass(ipts,ivm,iroot,iele) = biomass(ipts,ivm,iroot,iele) - &
1123            circ_class_biomass(ipts,ivm,icir,iroot,iele)*&
1124            circ_class_kill(ipts,ivm,icir,ifm,icut)
1125       biomass(ipts,ivm,ifruit,iele) = biomass(ipts,ivm,ifruit,iele) - &
1126            circ_class_biomass(ipts,ivm,icir,ifruit,iele)*&
1127            circ_class_kill(ipts,ivm,icir,ifm,icut)
1128
1129       ! Remove the dead biomass.  This includes everything except
1130       ! the belowground biomass. 
1131       circ_class_biomass(ipts,ivm,icir,ileaf,iele)=zero
1132       circ_class_biomass(ipts,ivm,icir,isapabove,iele)=zero
1133       circ_class_biomass(ipts,ivm,icir,iheartabove,iele)=zero
1134       circ_class_biomass(ipts,ivm,icir,iroot,iele)=zero
1135       circ_class_biomass(ipts,ivm,icir,ifruit,iele)=zero
1136
1137       ! Reduce the root biomass.
1138       biomass(ipts,ivm,isapbelow,iele) = biomass(ipts,ivm,isapbelow,iele) - &
1139            dead_sapwood_be*&
1140            circ_class_kill(ipts,ivm,icir,ifm,icut)
1141       circ_class_biomass(ipts,ivm,icir,isapbelow,iele)=circ_class_biomass(ipts,ivm,icir,isapbelow,iele)&
1142            *(un-coppice_kill_be_wood(ivm))
1143       biomass(ipts,ivm,iheartbelow,iele) = biomass(ipts,ivm,iheartbelow,iele) - &
1144            dead_heartwood_be*&
1145            circ_class_kill(ipts,ivm,icir,ifm,icut)
1146       circ_class_biomass(ipts,ivm,icir,iheartbelow,iele)=circ_class_biomass(ipts,ivm,icir,iheartbelow,iele)&
1147            *(un-coppice_kill_be_wood(ivm))
1148
1149    ENDDO ! loop over circ classes
1150
1151
1152    IF (bavard.GE.2) WRITE(numout,*) &
1153         'Leaving harvest_aboveground_wood_living_roots sapiens_kill.f90'
1154
1155  END SUBROUTINE harvest_aboveground_wood_living_roots
1156
1157!! ================================================================================================================================
1158!! SUBROUTINE   : second_coppice_cut
1159!!
1160!>\BRIEF        Simulates cutting a coppice between the first cut and the
1161!!              final cut.
1162!!
1163!! DESCRIPTION  : Moves dead biomass to the litter pool and the harvest pool
1164!!                where appropriate.  A few counters are reset.  The leftover
1165!!                biomass (the woody root system) is patitioned among enough
1166!!                individuals to give the same density as is found in coppice_dens,
1167!!                which is the density of stumps left after the first cut.
1168!!
1169!!                NOTE: All branches and stems are harvested.  Basically, this means
1170!!                      all aboveground woody biomass.  The fine roots
1171!!                      are killed and left to litter.  The sapwood and heartwood
1172!!                      belowground is left alive.  The leaves are killed, and the
1173!!                      labile and carbres pools are left untouched in the living
1174!!                      biomass.
1175!!
1176!! RECENT CHANGE(S):
1177!!
1178!! MAIN OUTPUT VARIABLE(S): ::bm_to_litter, ::harvest_pool
1179!!
1180!! REFERENCE(S)   :
1181!! - Bellassen, V., Le Maire, G., Dhote, J.F., Viovy, N., Ciais, P., 2010.
1182!! Modeling forest management within a global vegetation model – Part 1:
1183!! model structure and general behaviour. Ecological Modelling 221, 2458–2474.
1184!!
1185!! FLOWCHART    : None
1186!!\n
1187!_ ================================================================================================================================
1188 
1189  SUBROUTINE second_coppice_cut (npts, ipts, ivm, iele, ifm, &
1190       icut, harvest_pool_bound, resolution, veget_max, biomass, &
1191       ind, circ_class_biomass, circ_class_n, circ_class_kill, &
1192       bm_to_litter, harvest_pool, harvest_type, harvest_cut, &
1193       harvest_area, circ_class_dist, last_cut, coppice_dens)
1194
1195
1196    !! 0. Variable and parameter declaration
1197
1198    !! 0.1 Input variables
1199    ! NOTE: if harvest_pool_bound is not explicitely defined
1200    ! the model uses 1:ndia_harvest+2 as its dimensions. If
1201    ! you like to do so remember to change the counters in
1202    ! this routine
1203    INTEGER(i_std), INTENT(in)                                   :: npts               !! Domain size (-)
1204    INTEGER(i_std), INTENT(in)                                   :: ipts               !! current grid square
1205    INTEGER(i_std), INTENT(in)                                   :: ivm                !! current PFT
1206    INTEGER(i_std), INTENT(in)                                   :: iele               !! current element type
1207    INTEGER(i_std), INTENT(in)                                   :: ifm                !! current fm strategy
1208    INTEGER(i_std), INTENT(in)                                   :: icut               !! current cut time
1209    REAL(r_std), DIMENSION(0:ndia_harvest+1), INTENT(in)         :: harvest_pool_bound !! The boundaries of the diameter classes
1210                                                                                       !! in the wood harvest pools
1211                                                                                       !! @tex $(m)$ @endtex
1212    REAL(r_std),DIMENSION(npts,2), INTENT(in)                    :: resolution         !! The length of each grid-box in X and Y
1213                                                                                       !! @tex $(m)$ @endtex
1214    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: veget_max          !! Maximum fraction of vegetation type
1215    REAL(r_std), DIMENSION(ncirc), INTENT(in)                    :: circ_class_dist    !! The probability distribution of trees
1216                                                                                       !! in a circ class in case of a
1217                                                                                       !! redistribution (unitless).
1218    REAL(r_std), DIMENSION(npts,nvm),INTENT(in)                  :: coppice_dens       !! The density of a coppice at the first
1219                                                                                       !! cutting.
1220                                                                                       !! @tex $( 1 / m**2 )$ @endtex
1221
1222    !! 0.2 Output variables
1223
1224
1225    !! 0.3 Modified variables   
1226   
1227    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
1228                                             INTENT(inout)       :: biomass            !! Biomass @tex $(gC m^{-2}) $@endtex
1229    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)              :: ind                !! Stand level number of individuals
1230                                                                                       !! @tex $(m^{-2})$ @endtex
1231    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &                         
1232                                              INTENT(inout)      :: circ_class_biomass !! Biomass of the components of the model 
1233                                                                                       !! tree within a circumference
1234                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
1235    REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout)        :: circ_class_n       !! Number of individuals in each circ class
1236                                                                                       !! @tex $(m^{-2})$ @endtex
1237    REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times), &
1238                                                  INTENT(inout)  :: circ_class_kill    !! Number of trees within a circ that needs
1239                                                                                       !! to be killed @tex $(ind m^{-2})$ @endtex
1240    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
1241                                                  INTENT(inout)  :: bm_to_litter       !! Biomass transfer to litter
1242                                                                                       !! @tex $(gC m^{-2})$ @endtex
1243    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)               :: harvest_pool       !! The wood and biomass that have been
1244                                                                                       !! havested by humans @tex $(gC)$ @endtex
1245    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_type       !! Type of management that resulted
1246                                                                                       !! in the harvest (unitless)
1247    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_cut        !! Type of cutting that was used for the harvest
1248                                                                                       !! (unitless)
1249    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_area       !! Harvested area (m^{2})
1250    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: last_cut           !! Years since last thinning or clearcut (years)
1251
1252
1253    !! 0.4 Local variables
1254
1255    REAL(r_std), DIMENSION(ncirc,nparts)                         :: circ_class_biomass_new !! Biomass of the components of the model 
1256                                                                                       !! tree within a circumference
1257                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
1258    REAL(r_std), DIMENSION(ncirc)                                :: circ_class_n_new   !! Number of individuals in each circ class
1259                                                                                       !! @tex $(m^{-2})$ @endtex
1260    REAL(r_std), DIMENSION(ncirc)                                :: old_trees_left
1261    REAL(r_std)                                                  :: trees_needed
1262    REAL(r_std)                                                  :: scale_factor
1263    INTEGER                                                      :: icir, jcir
1264    !_ ================================================================================================================================
1265
1266    IF (bavard.GE.2) WRITE(numout,*) &
1267         'Entering second_coppice_cut sapiens_kill.f90'
1268
1269    !  Initialize check for surface area conservation
1270    !  Veget_max is a INTENT(in) variable and can therefore
1271    !  not be changed during the course of this subroutine
1272    !  No need to check whether the subroutine preserves the
1273    !  total surface area of the pixel.
1274
1275    IF( SUM(circ_class_kill(ipts,ivm,:,ifm, icut))&
1276         .GT. min_stomate ) THEN
1277
1278       IF(ld_forestry)THEN
1279          WRITE(numout,*) 'Doing a second coppice cut'
1280          WRITE(numout,*) 'ipts,ivm,ifm,icut ',ipts,ivm,ifm,icut
1281       ENDIF
1282
1283       ! Again, we harvest all aboveground biomass and leave the roots
1284       CALL harvest_aboveground_wood_living_roots(npts, ipts, ivm, iele, &
1285            ifm, icut, biomass, bm_to_litter, &
1286            circ_class_biomass, circ_class_kill, harvest_pool, &
1287            harvest_type, harvest_cut, harvest_area, &
1288            harvest_pool_bound, resolution, veget_max)
1289
1290       ! reset the time since the last cut
1291       last_cut(ipts,ivm)=0
1292
1293       ! We also reset circ_class_kill for this grid square, since
1294       ! we clearly don't need to do mortality if we have no more
1295       ! aboveground biomass left.
1296       circ_class_kill(ipts,ivm,:,:,:)=zero
1297
1298       ! Now we have a certain number density and a certain amount of
1299       ! biomass.  We want to have the same number of trees that we had
1300       ! during the first cut, since the stumps from the first part aren't
1301       ! going to change (maybe a couple will die, but not a large difference).
1302
1303       ! This code is copied from lpj_gap.f90.
1304       circ_class_n_new(:)=circ_class_dist(:)*ind(ipts,ivm)
1305       circ_class_biomass_new(:,:)=zero
1306
1307       ! now we populate the new classes
1308       old_trees_left(:)=circ_class_n(ipts,ivm,:)
1309       DO icir=1,ncirc
1310
1311          trees_needed=circ_class_n_new(icir)
1312          DO jcir=1,ncirc
1313
1314             IF(trees_needed .LE. zero)EXIT
1315
1316             IF(old_trees_left(jcir) .GE. trees_needed )THEN
1317                ! we can get all the trees we need from this class
1318                ! and don't have to continue searching
1319                circ_class_biomass_new(icir,:)=circ_class_biomass_new(icir,:)+&
1320                     circ_class_biomass(ipts,ivm,jcir,:,icarbon)*trees_needed
1321                old_trees_left(jcir)=old_trees_left(jcir)-trees_needed
1322                trees_needed = zero
1323                EXIT
1324             ELSE
1325
1326                ! the trees in this class are not sufficient, so we
1327                ! need to move all of them to the new class.
1328                circ_class_biomass_new(icir,:)=circ_class_biomass_new(icir,:)+&
1329                     circ_class_biomass(ipts,ivm,jcir,:,icarbon)*&
1330                     old_trees_left(jcir)
1331                trees_needed = trees_needed -old_trees_left(jcir)
1332                old_trees_left(jcir)=zero
1333
1334             ENDIF ! checking if we have enough trees
1335
1336          ENDDO ! jcir=1,ncirc
1337
1338       ENDDO ! icir=1,ncirc
1339
1340       ! right now, circ_class_biomass_new gives the total biomass
1341       ! in each class, but it should only be for a model tree.
1342       ! So let's normalize it.
1343       DO icir=1,ncirc
1344          circ_class_biomass_new(icir,:)=circ_class_biomass_new(icir,:)&
1345               /circ_class_n_new(icir)
1346       ENDDO
1347
1348       ! now update the variables that pass this information around
1349       DO icir=1,ncirc
1350          circ_class_biomass(ipts,ivm,icir,:,icarbon)=circ_class_biomass_new(icir,:)
1351          circ_class_n(ipts,ivm,icir)=circ_class_n_new(icir)
1352       ENDDO
1353
1354       !----- end of copied code..make a subroutine?
1355
1356
1357       ! reset the time since the last cut
1358       last_cut(ipts,ivm)=0
1359
1360       ! We also reset circ_class_kill for this grid square, since
1361       ! we clearly don't need to do mortality if we have no more
1362       ! aboveground biomass left.
1363       circ_class_kill(ipts,ivm,:,:,:)=zero
1364
1365       ! After a first coppicing, we will regrow a certain number
1366       ! of shoots per stool.  For a second coppice cut, our starting
1367       ! density is exactly the same as the first cut.
1368       
1369       scale_factor=coppice_dens(ipts,ivm)/SUM(circ_class_n(ipts,ivm,:))*&
1370            REAL(shoots_per_stool(ivm),r_std)
1371
1372       circ_class_n(ipts,ivm,:)=circ_class_n(ipts,ivm,:)*&
1373            scale_factor
1374       ind(ipts,ivm)=ind(ipts,ivm)*&
1375            scale_factor
1376
1377       ! in order not to change the total biomass, we need to
1378       ! make all our model trees smaller.
1379       circ_class_biomass(ipts,ivm,:,:,:)=circ_class_biomass(ipts,ivm,:,:,:)/&
1380            scale_factor
1381
1382    ENDIF ! any individuals to second cut coppice?
1383
1384    IF (bavard.GE.2) WRITE(numout,*) 'Leaving second_coppice_cut sapiens_kill.f90'
1385
1386  END SUBROUTINE second_coppice_cut
1387
1388!! ================================================================================================================================
1389!! SUBROUTINE   : first_coppice_cut
1390!!
1391!>\BRIEF        Transfer dead biomass to litter and wood to harvest pool.
1392!!
1393!! DESCRIPTION  : This is a total harvest of the aboveground biomass of a
1394!!                stand in such a way that the belowground wood continues to
1395!!                live.  This results in a sprouting in the following year of
1396!!                several shoots per stump.  Biomass is turned to litter, harvested,
1397!!                and left alive, depending on the pool.  The number of individuals
1398!!                goes up and the size of the model tree goes down.
1399!!
1400!!                NOTE: All branches and stems are harvested.  Basically,
1401!!                      this means all aboveground woody biomass.  The fine roots
1402!!                      are killed and left to litter.  The sapwood and heartwood
1403!!                      belowground is left alive.  The leaves are killed, and the
1404!!                      labile and carbres pools are left untouched in the living
1405!!                      biomass.
1406!!
1407!! RECENT CHANGE(S):
1408!!
1409!! MAIN OUTPUT VARIABLE(S): ::bm_to_litter, ::harvest_pool
1410!!
1411!! REFERENCE(S)   :
1412!!
1413!! FLOWCHART    : None
1414!!\n
1415!_ ================================================================================================================================
1416 
1417  SUBROUTINE first_coppice_cut (npts, ipts, ivm, iele, ifm, &
1418       icut, harvest_pool_bound, resolution, veget_max, biomass, &
1419       ind, circ_class_biomass, circ_class_n, circ_class_kill, &
1420       bm_to_litter, harvest_pool, harvest_type, harvest_cut, &
1421       harvest_area, last_cut)
1422
1423
1424    !! 0. Variable and parameter declaration
1425
1426    !! 0.1 Input variables
1427    ! NOTE: if harvest_pool_bound is not explicitely defined
1428    ! the model uses 1:ndia_harvest+2 as its dimensions. If
1429    ! you like to do so remember to change the counters in
1430    ! this routine
1431    INTEGER(i_std), INTENT(in)                                   :: npts               !! Domain size (-)
1432    INTEGER(i_std), INTENT(in)                                   :: ipts               !! current grid square
1433    INTEGER(i_std), INTENT(in)                                   :: ivm                !! current PFT
1434    INTEGER(i_std), INTENT(in)                                   :: iele               !! current element type
1435    INTEGER(i_std), INTENT(in)                                   :: ifm                !! current fm strategy
1436    INTEGER(i_std), INTENT(in)                                   :: icut               !! current cut time
1437    REAL(r_std), DIMENSION(0:ndia_harvest+1), INTENT(in)         :: harvest_pool_bound !! The boundaries of the diameter classes
1438                                                                                       !! in the wood harvest pools
1439                                                                                       !! @tex $(m)$ @endtex
1440    REAL(r_std),DIMENSION(npts,2), INTENT(in)                    :: resolution         !! The length of each grid-box in X and Y
1441                                                                                       !! @tex $(m)$ @endtex
1442    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: veget_max          !! Maximum fraction of vegetation type
1443
1444
1445    !! 0.2 Output variables
1446
1447
1448    !! 0.3 Modified variables   
1449
1450    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
1451                                             INTENT(inout)       :: biomass            !! Biomass @tex $(gC m^{-2}) $@endtex
1452    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)              :: ind                !! Stand level number of individuals
1453                                                                                       !! @tex $(m^{-2})$ @endtex
1454    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &                         
1455                                              INTENT(inout)      :: circ_class_biomass !! Biomass of the components of the model 
1456                                                                                       !! tree within a circumference
1457                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
1458    REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout)        :: circ_class_n       !! Number of individuals in each circ class
1459                                                                                       !! @tex $(m^{-2})$ @endtex
1460    REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times), &
1461                                                  INTENT(inout)  :: circ_class_kill    !! Number of trees within a circ that needs
1462                                                                                       !! to be killed @tex $(ind m^{-2})$ @endtex
1463    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
1464                                                  INTENT(inout)  :: bm_to_litter       !! Biomass transfer to litter
1465                                                                                       !! @tex $(gC m^{-2})$ @endtex
1466    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)               :: harvest_pool       !! The wood and biomass that have been
1467                                                                                       !! havested by humans @tex $(gC)$ @endtex
1468    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_type       !! Type of management that resulted
1469                                                                                       !! in the harvest (unitless)
1470    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_cut        !! Type of cutting that was used for the harvest
1471                                                                                       !! (unitless)
1472    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_area       !! Harvested area (m^{2})
1473    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: last_cut           !! Years since last thinning or clearcut (years)
1474
1475    !! 0.4 Local variables
1476
1477    !_ ================================================================================================================================
1478
1479    IF (bavard.GE.2) WRITE(numout,*) 'Entering first_coppice_cut sapiens_kill.f90'
1480
1481    ! Initialize check for surface area conservation
1482    ! Veget_max is a INTENT(in) variable and can therefore
1483    ! not be changed during the course of this subroutine
1484    ! No need to check whether the subroutine preserves the
1485    ! total surface area of the pixel.
1486
1487    ! What if we coppice for the first cut?
1488    IF( SUM(circ_class_kill(ipts,ivm,:,ifm, icut))&
1489         .GT. min_stomate ) THEN
1490
1491       IF(ld_forestry)THEN
1492          WRITE(numout,*) 'Doing a first coppice cut'
1493          WRITE(numout,*) 'ipts,ivm,ifm,icut ',ipts,ivm,ifm,icut
1494       ENDIF
1495
1496       CALL harvest_aboveground_wood_living_roots(npts, ipts, ivm, iele, &
1497            ifm, icut, biomass, bm_to_litter, &
1498            circ_class_biomass, circ_class_kill, harvest_pool, &
1499            harvest_type, harvest_cut, harvest_area, &
1500            harvest_pool_bound, resolution, veget_max)
1501
1502       ! reset the time since the last cut
1503       last_cut(ipts,ivm)=0
1504
1505       ! We also reset circ_class_kill for this grid square, since
1506       ! we clearly don't need to do mortality if we have no more
1507       ! aboveground biomass left.
1508       circ_class_kill(ipts,ivm,:,:,:)=zero
1509
1510       ! After a first coppicing, we will regrow a certain number
1511       ! of shoots per stool.  Therefore, our density starting
1512       ! next year will be much higher.
1513       circ_class_n(ipts,ivm,:)=circ_class_n(ipts,ivm,:)*&
1514            REAL(shoots_per_stool(ivm),r_std)
1515       ind(ipts,ivm)=ind(ipts,ivm)*&
1516            REAL(shoots_per_stool(ivm),r_std)
1517
1518       ! in order not to change the total biomass, we need to
1519       ! make all our model trees smaller.
1520       circ_class_biomass(ipts,ivm,:,:,:)=circ_class_biomass(ipts,ivm,:,:,:)/&
1521            REAL(shoots_per_stool(ivm),r_std)
1522
1523    ENDIF ! any trees to first cut coppice?
1524
1525
1526    IF (bavard.GE.2) WRITE(numout,*) 'Leaving first_coppice_cut sapiens_kill.f90'
1527
1528  END SUBROUTINE first_coppice_cut
1529
1530!! ================================================================================================================================
1531!! SUBROUTINE   : clearcut_harvest_aboveground_dead_roots
1532!!
1533!>\BRIEF        Transfer dead biomass to litter and wood to harvest pool
1534!!              for a whole stand.
1535!!
1536!! DESCRIPTION  : Moves dead biomass to the litter pool and the harvest pool
1537!!                where appropriate, assuming the whole stand is harvested.  Various
1538!!                counters are reset.
1539!!
1540!!                NOTE: If harvest_fraction is 1, all the aboveground biomass is
1541!!                      harvested.  If it's 1-branch_ratio, only stems are
1542!!                      harvested.  Belowground wood, roots, fruits, and
1543!!                      leaves all become litter.  A portion of aboveground
1544!!                      wood is put in the harvest pool.  Reserves and labile
1545!!                      carbon are assumed to be distributed evenly across
1546!!                      all non-heartwood components.
1547!!
1548!! RECENT CHANGE(S):
1549!!
1550!! MAIN OUTPUT VARIABLE(S): ::bm_to_litter, ::harvest_pool
1551!!
1552!! REFERENCE(S)   :
1553!!
1554!! FLOWCHART    : None
1555!!\n
1556!_ ================================================================================================================================
1557 
1558  SUBROUTINE clearcut_harvest_aboveground_dead_roots (npts, ipts, ivm, iele, &
1559       ifm, icut, biomass, ind, bm_to_litter, &
1560       circ_class_biomass, circ_class_kill, circ_class_n, harvest_pool, &
1561       harvest_type, harvest_cut, harvest_area, &
1562       harvest_pool_bound, resolution, veget_max, last_cut, age_stand, &
1563       mai_count, harvest_fraction)
1564
1565
1566    !! 0. Variable and parameter declaration
1567
1568    !! 0.1 Input variables
1569    ! NOTE: if harvest_pool_bound is not explicitely defined
1570    ! the model uses 1:ndia_harvest+2 as its dimensions. If
1571    ! you like to do so remember to change the counters in
1572    ! this routine
1573    INTEGER(i_std), INTENT(in)                                   :: npts               !! Domain size (-)
1574    INTEGER(i_std), INTENT(in)                                   :: ipts               !! current grid square
1575    INTEGER(i_std), INTENT(in)                                   :: ivm                !! current PFT
1576    INTEGER(i_std), INTENT(in)                                   :: iele               !! current element type
1577    INTEGER(i_std), INTENT(in)                                   :: ifm                !! current fm strategy
1578    INTEGER(i_std), INTENT(in)                                   :: icut               !! current cut time
1579    REAL(r_std), DIMENSION(0:ndia_harvest+1), INTENT(in)         :: harvest_pool_bound !! The boundaries of the diameter classes
1580                                                                                       !! in the wood harvest pools
1581                                                                                       !! @tex $(m)$ @endtex
1582    REAL(r_std),DIMENSION(npts,2), INTENT(in)                    :: resolution         !! The length of each grid-box in X and Y
1583                                                                                       !! @tex $(m)$ @endtex
1584    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                   :: veget_max          !! Maximum fraction of vegetation type
1585    REAL(r_std),INTENT(in)                                       :: harvest_fraction   !! The fraction of the aboveground
1586                                                                                       !! biomass that is harvested.  This is
1587                                                                                       !! generally 1-branch_ratio if you only
1588                                                                                       !! want to harvest the stems, and 1 if
1589                                                                                       !! you want to harvest everything.
1590
1591
1592    !! 0.2 Output variables
1593
1594
1595    !! 0.3 Modified variables   
1596    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
1597                                             INTENT(inout)       :: biomass            !! Biomass @tex $(gC m^{-2}) $@endtex
1598    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)              :: ind                !! Stand level number of individuals
1599                                                                                       !! @tex $(m^{-2})$ @endtex
1600    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &                         
1601                                              INTENT(inout)      :: circ_class_biomass !! Biomass of the components of the model 
1602                                                                                       !! tree within a circumference
1603                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
1604    REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout)        :: circ_class_n       !! Number of individuals in each circ class
1605                                                                                       !! @tex $(m^{-2})$ @endtex
1606    REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times), &
1607                                                  INTENT(inout)  :: circ_class_kill    !! Number of trees within a circ that needs
1608                                                                                       !! to be killed @tex $(ind m^{-2})$ @endtex
1609    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
1610                                                  INTENT(inout)  :: bm_to_litter       !! Biomass transfer to litter
1611                                                                                       !! @tex $(gC m^{-2})$ @endtex
1612    REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout)               :: harvest_pool       !! The wood and biomass that have been
1613                                                                                       !! havested by humans @tex $(gC)$ @endtex
1614    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_type       !! Type of management that reulted
1615                                                                                       !! in the harvest (unitless)
1616    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_cut        !! Type of cutting that was used for the harvest
1617                                                                                       !! (unitless)
1618    REAL(r_std), DIMENSION(:,:), INTENT(inout)                   :: harvest_area       !! Harvested area (m^{2})
1619    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: last_cut           !! Years since last thinning or clearcut (years)
1620    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: age_stand          !! Age of stand (years)
1621    INTEGER(i_std), DIMENSION(npts,nvm),INTENT(inout)            :: mai_count          !! The number of times we've
1622                                                                                       !! calculated the volume increment
1623                                                                                       !! for a stand
1624    !! 0.4 Local variables
1625
1626    !_ ================================================================================================================================
1627
1628    IF (bavard.GE.2) WRITE(numout,*) &
1629         'Entering clearcut_harvest_aboveground_dead_roots sapiens_kill.f90'
1630
1631    ! Initialize check for surface area conservation
1632    ! Veget_max is a INTENT(in) variable and can therefore
1633    ! not be changed during the course of this subroutine
1634    ! No need to check whether the subroutine preserves the
1635    ! total surface area of the pixel.
1636
1637    IF( SUM(circ_class_kill(ipts,ivm,:,ifm, icut))&
1638         .GT. min_stomate ) THEN
1639
1640       CALL harvest_aboveground_dead_roots(npts, ipts, ivm, iele, &
1641            ifm, icut, biomass, ind, bm_to_litter, &
1642            circ_class_biomass, circ_class_kill, circ_class_n, harvest_pool, &
1643            harvest_type, harvest_cut, harvest_area, &
1644            harvest_pool_bound, resolution, veget_max, &
1645            harvest_fraction)
1646
1647       ! We have to reset the age of the stand and the time since the
1648       ! last human intervention.
1649       age_stand(ipts,ivm)=0
1650       last_cut(ipts,ivm)=0
1651       mai_count(ipts,ivm)=0
1652
1653       ! We also reset circ_class_kill for this grid square, since
1654       ! we clearly don't need to do mortality if we have no more
1655       ! trees left.
1656       circ_class_kill(ipts,ivm,:,:,:)=zero
1657
1658    ENDIF
1659
1660    IF (bavard.GE.2) WRITE(numout,*) &
1661         'Leaving clearcut_harvest_aboveground_dead_roots sapiens_kill.f90'
1662
1663  END SUBROUTINE clearcut_harvest_aboveground_dead_roots
1664
1665END MODULE sapiens_kill
Note: See TracBrowser for help on using the repository browser.