-
Petra Lasch-Born authoredPetra Lasch-Born authored
target_thin.f 17.46 KiB
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutine *!
!* target thinning - *!
!* thinning routine with given values of biomass per *!
!* thinning year as target values *!
!* targetm i given in kg DW/ha *!
!* *!
!* Copyright (C) 1996-2018 *!
!* Potsdam Institute for Climate Impact Reserach (PIK) *!
!* Authors and contributors see AUTHOR file *!
!* This file is part of 4C and is licensed under BSD-2-Clause *!
!* See LICENSE file or under: *!
!* http://www.https://opensource.org/licenses/BSD-2-Clause *!
!* Contact: *!
!* https://gitlab.pik-potsdam.de/foresee/4C *!
!* *!
!*****************************************************************!
SUBROUTINE target_thinning(i)
use data_stand
use data_manag
use data_simul
use data_species
use data_par
implicit none
real :: targetm ! target value of stem biomass
real :: dbhmin=0, &
dbhmin_us = 0, &
wpa=0, & ! Weibull parameter
wpa_us , &
wpb=0, & ! -"-
wpb_us, &
wpc=0, & ! -"-
d63=0, &
d63_us, &
help=0, &
pequal, &
tdbh=0, &
bas_area=0, &
bas_help=0., &
rtarget_help=0, &
target_help1=0,&
dbh_h =0, &
db_l = 0., &
db_u = 0., &
d_est=0., &
w_kb=0., &
stembiom, &
stembiom_us = 0. , &
stembiom_re = 0. , &
stembiom_all = 0. , &
diff, &
mdiam, &
mdiam_us
integer :: nrmin, &
lowtree, &
undertree, &
flagth, &
taxnr, &
counth, &
min_id, & max_id, &
ih1,ih2,ncoh, &
coun1
! auxilary for thinning routine 4: selective thinning
integer :: count,i, &
idum , third, ipot, isel, ih
integer,dimension(0:anz_coh) :: cohl
integer, dimension(anz_coh) :: id_pot
real :: h1, h2
real,external :: gasdev
real:: ran0
! reacalculation of target to kg DW/patch
h1 = 0.
h2 = 0.
count = 0
cohl = -1
flagth = 0
coun1 = 0
help=0.
lowtree=0
undertree = 0
anz_tree_dbh = 0
bas_area = 0.
! stem biomass of overstorey
stembiom = 0.
! stem biomass of understorey
stembiom_us = 0.
stembiom_all = 0.
if (time.eq.73.and. ip .eq.87) then
stembiom = 0
end if
taxnr = thin_spec(i)
mdiam = 0.
mdiam_us = 0.
! calculation of mean diameter (correspondung to med_diam) and basal area of stand
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Modification for V Kint: no test for diameter
IF((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.0) THEN
! overstorey
stembiom = stembiom + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
help = help + zeig%coh%ntreea*(zeig%coh%diam**2)
bas_area = bas_area + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
if( zeig%coh%diam>0) then
anz_tree_dbh = anz_tree_dbh + zeig%coh%ntreea
mdiam = mdiam + zeig%coh%ntreea * (zeig%coh%diam**2)
end if
! Trees with DBH=0 for population and per species; Baeume mit DBH =0 fuer Bestand und pro Spezie
ELSE IF( (zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) THEN
! seedings/regeneration
stembiom_re = stembiom_re + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
lowtree = lowtree + zeig%coh%ntreea
ELSE if((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) THEN
! understorey
stembiom_us = stembiom_us + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
mdiam_us = mdiam_us + zeig%coh%ntreea * (zeig%coh%diam**2)
undertree = undertree + zeig%coh%ntreea
ENDIF
zeig => zeig%next
ENDDO
! mean diamteer for over and understorey
stembiom_all = stembiom + stembiom_usif(anz_tree_dbh.ne.0) mdiam = sqrt(mdiam/real(anz_tree_dbh))
if(undertree.ne.0) mdiam_us = sqrt(mdiam_us/undertree)
third = nint(anz_tree_dbh*0.333333)
anz_tree_ha = nint(anz_tree_dbh*10000./kpatchsize)
IF(anz_tree>0)THEN
if(lowtree<anz_tree) help = sqrt(help/(anz_tree-lowtree))
ENDIF
! setting of aux. variable target_help
rtarget_help = stembiom_all
! tending
if(thin_tysp(i).eq.4.or.(stembiom_re.ne.0. .and. stembiom_all.eq.0)) then
rtarget_help = stembiom_re
end if
! Umrechnung in Biomasse pro patch in kg
targetm = target_mass(i)*1000*kpatchsize/10000.
! target value of biomass
if(thin_tysp(i).eq.4 .or.(stembiom_re.ne.0. .and. stembiom_all.eq.0) ) then
! tending
targetm = stembiom_re - targetm*stembiom_re
else
end if
if( targetm.eq.1) targetm = 0.
! targetm (kg DW/patch)
! cuttting
if (targetm.eq.0.)then
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
!tending of regeneration
else if(thin_tysp(i).eq.4) then
min_id = 1000
max_id = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1) then
ih1 = zeig%coh%ident
if(ih1.lt.min_id) min_id = ih1
ih2 = zeig%coh%ident
if (ih2.gt.max_id) max_id = ih2
end if
zeig=> zeig%next
end do
target_help1 = 0.
do
call random_number(pequal)
ncoh = min_id +(max_id-min_id)*pequal
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1.and. zeig%coh%ident.eq.ncoh ) then
zeig%coh%ntreea = zeig%coh%ntreea - 1
zeig%coh%nta = zeig%coh%nta-1
zeig%coh%ntreem = zeig%coh%ntreem +1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit end if
zeig=>zeig%next
end do
diff = targetm - rtarget_help
if(diff.lt.0.01) exit
end do
else IF ( targetm .ne. 0.) then
if(target_mass(i).lt.1.) then
targetm = target_mass(i) * rtarget_help
end if
! different thinnings from below and above
select case(thin_tysp(i))
case(1)
! moderate lower thinning;
d_est = 1.02
w_kb = 1.8
case(2)
! intensive lower thinning;
d_est = 1.03
w_kb = 1.5
case(3)
! high thinning;
d_est = 1.04
w_kb = 1.2
end select
! calculation of Weibull-Parameter
call min_dbh_overs(nrmin,dbhmin,taxnr)
call min_dbh_unders(nrmin,dbhmin_us, taxnr)
bas_help = bas_area
wpa = dbhmin
wpa_us = dbhmin_us
d63 = mdiam*d_est
d63_us = mdiam_us * d_est
wpb = (d63 - wpa)/ w_kb
wpb_us = (d63_us-wpa_us)/w_kb
wpc = 2
if (thin_tysp(i).eq. 3) then
! starting with overstorey!, continuing with the understorey
if(targetm.lt.(stembiom_all-stembiom)) then
! total removing of overstorey
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
! defining the remaining part of stem biomass which has to be remove
rtarget_help = stembiom_us
! understorey
!selection of trees for thinning
if(mdiam_us.ne.0) then
! start understorey thinning
do
call random_number(pequal)
tdbh = wpa_us + wpb_us*(-log(1.-pequal))**(1./wpc)
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt.10000) exit
end if
zeig=> zeig%next
END DO ! list of potential thinned tees cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
count = count +1
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(count.gt.100000) exit
end do ! understorey thinning
end if ! mediam_us.ne.0
else ! targetm.lt.(stembiom_all-stembiom)
!selection of trees for thinning
do
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.2*dbh_h
db_u = dbh_h + 0.2*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
call random_number(pequal)
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
! write(1234,*) time, ipot, pequal, isel
! if(isel.eq.0) isel =1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(coun1.gt.100000) exit
end do ! total thinning
end if !targetm.lt.(stembiom_all-stembiom)
end if ! thintype 3
if(thin_tysp(i).eq.1.or.thin_tysp(i).eq.2) then
if(targetm.lt.(stembiom_all-stembiom_us)) then
! total removing of understorey
zeig => pt%first
DO IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.2) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
! definging the remaining part of stem biomass which has to remove
rtarget_help = stembiom
if(mdiam.ne.0) then
! additional thinning from the overstorey
counth = 0
do
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.3*dbh_h
db_u = dbh_h + 0.3*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(counth.gt.100000) exit
end do ! overstorey thinning
end if ! mdiam.ne.0
else ! targtem.lt.(stembiom_all-stembiom_us)
! first thinning from the understorey
!selection of trees for thinning
if(mdiam_us.ne.0) then
do
call random_number(pequal)
tdbh = wpa_us + wpb_us*(-log(1.-pequal))**(1./wpc)
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth .gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1 .or. (stembiom_all-stembiom_us).eq.rtarget_help) exit
if(coun1.gt.100000) exit
end do ! understorey thinning
end if ! mdiam_us
end if
end if !! thin_tysp.eq.1 or. thin-tysp.eq.2
END IF ! all thinnings and tending
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem>0.and.zeig%coh%species.eq.taxnr) then
! all parts without stems of trees are input for litter
h1 = zeig%coh%ntreem*zeig%coh%x_fol
h2 = h2 + h1
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreem*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreem*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
endif
zeig=>zeig%next
enddo
END SUBROUTINE target_thinning