!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* - Calculation of annual allocation of NPP (SR PARTITION) *!
!* - Calculation of annual allocation of NPP of soil *!
!* vegetation (PARTITION_SV *!
!* - Calculation of diameter at breast height (SR CALC_DBH) *!
!* *!
!* 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 PARTITION *!
!****************************!
SUBROUTINE PARTITION( p )
!*** Declaration part ***!
USE data_out
USE data_par
USE data_stand
USE data_species
USE data_simul
IMPLICIT NONE
REAL :: lambdaf = 0., & ! partitioning functions
lambdas = 0., &
lambdar = 0., &
lambdac = 0., &
lambdaSum = 0.,& ! sum of the above three lambdas
NPP = 0., & ! annual NPP
F = 0., & ! state variables: foliage,
S = 0., & ! sapwood,
H = 0., & ! heartwood
R = 0., & ! fine roots,
B = 0., & ! bole height,
Ahb = 0., & ! cross sectional area heartwood at tree base
hs = 0., & ! sapwood height
Ht = 0., & ! total tree height
Asw = 0., & ! cross sectional area of sapwood in bole
DBH = 0., & ! tree diameter at breast height (DBH)
FNew, SNew, & ! new states
RNew, BNew, &
HtNew, &
HNew, Ahbnew, &
sigmaf = 0., & ! current leaf activity rate
sigman = 0., & ! current root activity rate
ar = 0., & ! aux vars for partitioning functions
as = 0., &
ac = 0., &
betar = 0., &
betas = 0., &
aux = 0., &
Fmax, & ! determines whether height growth or not
rsap, & ! auxiliary variable for height growth determination
growthrate ! height growthrate depends on relative light regime in the middle of the canopy
REAL :: Sf, & ! senescence rates
Ss, &
Sr, &
Gf, & ! growth rates
Gs, &
Gr
real :: DBH_help
REAL :: leaf_N_conc, & ! last years N concentration in leaves gN kgDM
tbc_root_Ndemand, & ! N demand for ghrowth of fine roots, branches and coarse roots g tree-1
Nredfak ! reduction factor for N allocation to fine roots, branches and coarse roots
TYPE(Coh_Obj) :: p ! pointer to cohort list
REAL :: term1, &
a1, a2, a3, & ! coefficients of quadratic equation
x1 = 0., &
x2 = 0. ! solutions of quadratic equation
real :: Fmax_old
! if this cohort is mistletoe infected, reduce NPP by mistletoe-specific demand
! demand is defined in PARTITION_MI. as mistletoe is always 1st cohort, the demand of mistletoe is calculated before the reduction here
if (p%coh%mistletoe.eq.1) then
p%coh%NPP = p%coh%NPP-(NPP_demand_mistletoe/p%coh%ntreea)
endif
ns = p%coh%species
F = p%coh%x_fol
Fmax = p%coh%Fmax
S = p%coh%x_sap
R = p%coh%x_frt
H = p%coh%x_hrt
B = p%coh%x_hbole
NPP = p%coh%NPP
Ht = p%coh%height
Ahb = p%coh%x_Ahb
Sf = p%coh%sfol
Ss = p%coh%ssap
Sr = p%coh%sfrt
hs = p%coh%x_hsap
Asw = p%coh%Asapw
Fmax_old = Fmax
DBH_help = p%coh%diam
if (flag_end.eq.1) then
p%coh%notViable = .TRUE.
flag_end = 0
end if
if(p%coh%notViable.neqv..TRUE.) then
select case (flag_folhei)
case (1,4)
spar(ns)%pha = spar(ns)%pha_v1 * spar(ns)%pha_v3 * &
(F)**(-1-spar(ns)%pha_v3)/(spar(ns)%pha_v2+(F)**(-spar(ns)%pha_v3))**2.
case (2)
rsap=Asw/(Asw+Ahb)
spar(ns)%pha = 2.*spar(ns)%crown_a/(pi**0.5*(rsap*spar(ns)%pnus)**1.5*F**0.5)
case (3)
! this version only for tests and pine trees
spar(ns)%pha = (3500*(10.+F**0.9)-(0.9*3500.*F**0.9))/(10.+F**0.9)**2
end select ! flag_folhei
! only allocate if enough NPP is available
IF (NPP>1.0E-9) THEN
select case (flag_folhei)
case (0)
growthrate=spar(ns)%pha*spar(ns)%pha_coeff1 + spar(ns)%pha*spar(ns)%pha_coeff2*(1./p%coh%IrelCan-1.)
case (1,3)
growthrate=spar(ns)%pha + spar(ns)%pha*(1./MAX(p%coh%IrelCan,0.25)-1.)
case (2)
growthrate=spar(ns)%pha + spar(ns)%pha*(1.-p%coh%IrelCan)*5.
case (4)
growthrate=spar(ns)%pha *0.5/MAX(p%coh%IrelCan,0.25)
end select ! flag_folhei
sigmaf = NPP/F
! calculate root activity based on drought index
! test of a relationship which modifies fine root leaf ratio with shade tolerance:
IF (flag_sign.eq.1 .or. flag_sign.eq.11) THEN
term1 = spar(ns)%sigman * 10. * (((5.-spar(ns)%stol)*1.-p%coh%crown_area) / (5.-spar(ns)%stol)*1.)
sigman = amax1(term1,spar(ns)%sigman) * p%coh%drIndAl/p%coh%nDaysGr
ELSE
sigman = spar(ns)%sigman * p%coh%drIndAl / p%coh%nDaysGr
END IF
if (flag_sign .eq. 0 .or. flag_sign .eq. 1) then
! auxiliary variables for fine roots
ar = spar(ns)%pcnr * sigmaf / sigman
betar = (Sr - R + ar*(F-Sf)) / NPP
! auxiliary variables for sapwood
as = spar(ns)%prhos / spar(ns)%pnus
aux = 2.*(B+p%coh%deltaB) + Ht
betas = ( (as/3.)*(aux - growthrate*Sf) * (F-Sf) + Ss - S ) / NPP
! solve quadratic equation for lambdaf
term1 = (1.+spar(ns)%alphac)
a1 = term1 * as/3. * growthrate * NPP
a2 = 1.0 + ar + term1 * as/3. * (aux + growthrate*(F-2.*Sf))
a3 = term1*betas + betar - 1.
x1 = (-a2 + SQRT( a2*a2 - 4.*a1*a3) ) / (2.*a1)
x2 = (-a2 - SQRT( a2*a2 - 4.*a1*a3) ) / (2.*a1)
lambdaf = x1
if (lambdaf .le. 0. .or. lambdaf .gt. 1.) then
lambdaf = 0.5
lambdar = 0.5
lambdas = 0.
lambdac = 0.
else
! calculate coefficients for sapwood and roots
lambdar = ar * lambdaf + betar;
lambdas = as/3. * (aux + growthrate*(F+lambdaf*NPP-2.*Sf)) * lambdaf + betas
lambdac = spar(ns)%alphac * lambdas
! check consistency of calculation, i.e. no negative values
IF(lambdas < 0. .or. lambdas .gt. 1.) THEN
lambdas = 0.
lambdac = 0.
lambdaf = (1.-betar)/(ar+1)
lambdar = 1.-lambdaf
if (lambdaf .le. 0. .or. lambdaf .gt. 1.) then
lambdaf = 0.5
lambdar = 0.5
else if (lambdar<0) then
lambdar=0.
lambdaf=1.
end if
ELSE
! reduced allocation schemes for lamdaf<0. or lamdar<0. still to be added
lambdaf = AMAX1( lambdaf, 0. )
lambdar = AMAX1( lambdar, 0. )
! warrant that lambdaSum = 1 if balance can not be achieved this time step
lambdaSum = lambdaf + (1.+spar(ns)%alphac)*lambdas + lambdar
lambdaf = lambdaf / lambdaSum
lambdas = lambdas / lambdaSum
lambdar = lambdar / lambdaSum
lambdac = lambdac / lambdaSum
lambdaSum = lambdaf + (1.+spar(ns)%alphac)*lambdas + lambdar ! for debugging only
END IF
end if ! lambdaf .le. 0.
else ! flag_sign = 10, 11
! auxiliary variables for fine roots
ar = spar(ns)%pcnr * sigmaf / sigman
betar = (Sr - ar*Sf) / NPP
! auxiliary variables for sapwood
as = spar(ns)%prhos / spar(ns)%pnus
betas = (Ss - 2.*as*hs*Sf ) / NPP
! auxiliary variables for coarse roots, twigs and branches
ac = spar(ns)%alphac
! linear equation system in lamda(i)
term1 = 1. + ar + 2.*as*hs*(1+ac)
lambdaf = 1. - (1.+ac)*betas - betar
lambdaf = lambdaf / term1
lambdar = ar * lambdaf + betar
lambdas = 2.*as*hs * lambdaf + betas
lambdac = ac * lambdas
if (lambdaf .le. 0. .or. lambdaf .gt. 1.) then
lambdaf = 0.5
lambdar = 0.5
lambdas = 0.
lambdac = ac * lambdas
else
! calculate coefficients for sapwood and roots
lambdar = ar * lambdaf + betar;
lambdas = 2.*as*hs * lambdaf + betas
lambdac = ac * lambdas
! check consistency of calculation, i.e. no negative values
IF(lambdas < 0. .or. lambdas .gt. 1.) THEN
lambdas = 0.
lambdac = 0.
lambdaf = (1.-betar)/(ar+1)
lambdar = 1.-lambdaf
if (lambdaf .le. 0. .or. lambdaf .gt. 1.) then
lambdaf = 0.5
lambdar = 0.5
else if (lambdar<0) then
lambdar=0.
lambdaf=1.
end if
ELSE
! reduced allocation schemes for lamdaf<0. or lamdar<0. still to be added
lambdaf = AMAX1( lambdaf, 0. )
lambdar = AMAX1( lambdar, 0. )
! warrant that lambdaSum = 1 if balance can not be achieved this time step
lambdaSum = lambdaf + (1.+spar(ns)%alphac)*lambdas + lambdar
lambdaf = lambdaf / lambdaSum
lambdas = lambdas / lambdaSum
lambdar = lambdar / lambdaSum
lambdac = lambdac / lambdaSum
lambdaSum = lambdaf + (1.+spar(ns)%alphac)*lambdas + lambdar ! for debugging only
END IF
end if ! lambdaf .le. 0.
endif ! flag_sign
ELSE
lambdaf = 0.
lambdas = 0.
lambdar = 0.
END IF ! IF NPP < 1.0E-09
! gross growth rates of compartments
Gf = lambdaf * NPP
Gr = lambdar * NPP
Gs = lambdas * NPP
p%coh%gfol = Gf
p%coh%gfrt = Gr
p%coh%gsap = Gs
p%coh%x_crt = p%coh%x_crt + Gs*spar(ns)%alphac*spar(ns)%cr_frac
p%coh%x_tb = p%coh%x_tb + Gs*spar(ns)%alphac*(1.-spar(ns)%cr_frac)
! update of state vector
FNew = F + Gf - Sf
SNew = S + Gs - Ss
RNew = R + Gr - Sr
Hnew = H + Ss
AhbNew= Ahb + Asw*spar(ns)%pss
! check whether height growth or not
IF (lambdas == 0.OR.FNew<Fmax) THEN ! treat special case where there is no height growth
HtNew = Ht
ELSE
! height growth depending on the relative light regime in the middle of the canopy
HtNew = Ht + growthrate * (FNew-Fmax)
Fmax=FNew
ENDIF
BNew = B+p%coh%deltaB
! copy back to original variables
p%coh%Fmax = Fmax
p%coh%x_fol = FNew
p%coh%x_sap = SNew
p%coh%x_frt = RNew
p%coh%x_hrt = HNew
p%coh%height = HtNew
p%coh%x_hbole= BNew
p%coh%x_Ahb = AhbNew
CALL CALC_DBH(BNew,Htnew,Snew,Hnew,Ahbnew,p%coh%Ahc,p%coh%ident,DBH,p%coh%dcrb,hs,Asw)
if (flag_end.eq.1) then
DBH = p%coh%diam
p%coh%notViable = .TRUE.
flag_end = 0
end if
! Monitoring of current values
if (time_out .gt. 0 .and. flag_cohout .eq. 2) then
CALL OUT_ALL( p%coh%ident, p%coh%ntreea, NPP, DBH, growthrate,Fnew,Fmax_old,Htnew, lambdaf,lambdas,lambdar,lambdac,x1,x2)
endif
p%coh%x_hsap = hs
p%coh%diam = DBH ! This is the new value
p%coh%Asapw = Asw
p%coh%jrb = (DBH-DBH_help)*10/2
if(((DBH-DBH_help)*10/2).lt.0.) p%coh%jrb = 0.
! variables required by mortality submodel
p%coh%fol_inc = Gf - Sf
p%coh%bio_inc = NPP - Sf - (1.+spar(ns)%alphac)*Ss - Sr
p%coh%stem_inc = Gs ! deltaH + deltaS = Ss + Gs - Ss
p%coh%frt_inc = Gr - Sr ! fine root increment
p%coh%totBio = p%coh%x_fol + (1.+spar(ns)%alphac)*(p%coh%x_sap + p%coh%x_hrt) + p%coh%x_frt
p%coh%notViable = (FNew <= 0.) .OR. (SNew <= 0.) .OR. &
(RNew <= 0.) .OR. (Htnew <= Bnew)
! Nitrogen dynamics:
leaf_N_conc = p%coh%N_fol/F
! Simple model: all (sap)wood grows with CN-ratios of branches, twigs and coarse roots.
! When sapwood senesces N is reallocated and the new heart wood is at the level of stem CN-ratios.
! Branches, twigs and coarse roots do not senesce
! first step nitrogen related processes: N in litter, N-recallocation
p%coh%N_pool = p%coh%N_pool + Sf/F*p%coh%N_fol*spar(ns)%reallo_fol &
+ Sr*cpart/spar(ns)%cnr_frt*1000.* spar(ns)%reallo_frt &
+ Ss*cpart *1000. * (1/spar(ns)%cnr_tbc - 1/spar(ns)%cnr_stem)
p%coh%N_fol = p%coh%N_fol*(1-Sf/F)
p%coh%litC_fol = p%coh%litC_fol + p%coh%ntreea * Sf * cpart
p%coh%litC_frt = p%coh%litC_frt + p%coh%ntreea * Sr * cpart
! Species specific N content and reallocation factor (see species.par)
! Caution: tbc mortallity is not a litter compartment; it is assigned as heartwood
p%coh%litN_fol = p%coh%litN_fol + p%coh%ntreea * Sf * cpart * (1.-spar(ns)%reallo_fol) / spar(ns)%cnr_fol
p%coh%litN_frt = p%coh%litN_frt + p%coh%ntreea * Sr * cpart * (1.-spar(ns)%reallo_frt) / spar(ns)%cnr_frt
! second step: allocation of N to new growth
! before bud-break allocation to leaves is 50% of the N content of last years foliage
tbc_root_Ndemand = Gs*cpart *kg_in_g / spar(ns)%cnr_tbc + Gr* cpart/spar(ns)%cnr_frt*kg_in_g
IF(tbc_root_Ndemand + Gf*p%coh%med_sla*0.5 > p%coh%N_pool) THEN
if (tbc_root_Ndemand .gt. 1E-8) then
Nredfak = AMAX1((p%coh%N_pool-Gf*p%coh%med_sla*0.5) / tbc_root_Ndemand,0.) ! Division by zero possible
else
Nredfak = 0.
endif
tbc_root_Ndemand = tbc_root_Ndemand*Nredfak
ENDIF
p%coh%N_pool = p%coh%N_pool - tbc_root_Ndemand
IF(p%coh%N_pool < Gf*0.5*leaf_N_conc) THEN
p%coh%N_fol = p%coh%N_fol + p%coh%N_pool
p%coh%N_pool = 0.
ELSE
p%coh%N_fol = p%coh%N_fol + Gf*0.5*leaf_N_conc
p%coh%N_pool = p%coh%N_pool - Gf*0.5*leaf_N_conc
ENDIF
end if
END SUBROUTINE PARTITION
!*******************************!
!* SUBROUTINE PARTITION_SV *!
!*******************************!
SUBROUTINE PARTITION_SV( p )
!*** Declaration part ***!
USE data_par
USE data_stand
USE data_species
USE data_simul
IMPLICIT NONE
REAL :: lambdaf = 0., & ! partitioning functions
lambdas = 0., &
lambdar = 0., &
NPP = 0., & ! annual NPP
F = 0., & ! state variables: foliage,
S = 0., & ! sapwood,
R = 0., & ! fine roots,
Ht = 0., & ! total tree height
FNew, SNew, & ! new states
RNew, &
sigmaf = 0., & ! current leaf activity rate
sigman = 0. ! current root activity rate
REAL :: Sf, & ! senescence rates
Ss, &
Sr, &
Gf, & ! growth rates
Gs, &
Gr
REAL :: FRsum
REAL :: tbc_root_Ndemand, & ! N demand for ghrowth of fine roots, branches and coarse roots g tree-1
Nredfak ! reduction factor for N allocation to fine roots, branches and coarse roots
REAL, EXTERNAL :: f_lf, df_lf, ddf_lf
INTEGER :: flag_SV_allo, &
rnum
TYPE(Coh_Obj) :: p ! pointer to cohort list
ns = p%coh%species
F = p%coh%x_fol
S = p%coh%x_sap
R = p%coh%x_frt
NPP = p%coh%NPP
Ht = p%coh%height
Sf = p%coh%sfol
Ss = p%coh%ssap
Sr = p%coh%sfrt
! choice of allocation model. 0 = constant allocation factors, 1 = allometric model
flag_SV_allo = 1
! only allocate if enough NPP is available
IF (NPP>1.0E-9) THEN
! calculate leaf activity based on net PS and leaf mass
sigmaf = NPP/F
! calculate root activity based on drought index
! test of a relationship which modifies fine root leaf ratio with shade tolerance
IF (flag_sign.eq.1) THEN
sigman = amax1(spar(ns)%sigman*10*(((5.-spar(ns)%stol)*1.-p%coh%crown_area)/(5.-spar(ns)%stol)*1.),spar(ns)%sigman) * p%coh%drIndAl / p%coh%nDaysGr
ELSE
sigman = spar(ns)%sigman * p%coh%drIndAl / p%coh%nDaysGr
END IF
M_avail=(NPP+F-Sf+R-Sr+S-Ss)/kpatchsize
IF(flag_SV_allo==0) THEN
! the parameters pdiam in the species.par file are used for allocation fractions
lambdaf=spar(ns)%pdiam1
lambdar=spar(ns)%pdiam2
lambdas=spar(ns)%pdiam3
ELSE
FRsum=(F+R)/kpatchsize
CALL newt (FRsum, f_lf, df_lf, ddf_lf, 1.e-6, 100, rnum)
IF(FRsum>M_avail .and. .not.flag_mult8910) CALL error_mess(time,'no solution found for allocation for groundvegetation cohort, rnum: ',real(rnum))
IF(rnum==-1) THEN
if (.not.flag_mult8910) CALL error_mess(time,'no solution found for allocation for groundvegetation cohort: ',real(p%coh%ident))
lambdaf=0.4
lambdar=0.4
lambdas=0.2
ELSE
lambdaf=(FRsum)/M_avail/2.
lambdar=(FRsum)/M_avail/2.
lambdas=1.-lambdaf-lambdar
ENDIF
ENDIF
END IF ! IF NPP < 1.0E-09
! gross growth rates of compartments
Gf = lambdaf * M_avail*kpatchsize -F +Sf
Gr = lambdar * M_avail*kpatchsize -R +Sr
Gs = lambdas * M_avail*kpatchsize -S +Ss
! preliminary solution for permanent seeding
IF(lambdaf * M_avail < 1.e-4) THEN
Gf = Gf + 1.e-4*kpatchsize
ENDIF
p%coh%gfol = Gf
p%coh%gfrt = Gr
p%coh%gsap = Gs
! update of state vector
FNew = F + Gf - Sf
SNew = S + Gs - Ss
RNew = R + Gr - Sr
p%coh%x_fol = FNew
p%coh%x_sap = SNew
p%coh%x_frt = RNew
! determine litter production from plant turnover rates
! first step nitrogen related processes: N in litter, N-recallocation
p%coh%N_pool = p%coh%N_pool + Sf/F*p%coh%N_fol*spar(ns)%reallo_fol &
+ Sr*cpart/spar(ns)%cnr_frt*1000.* spar(ns)%reallo_frt &
+ Ss*cpart *1000. * (1/spar(ns)%cnr_tbc - 1/spar(ns)%cnr_stem)
p%coh%N_fol = p%coh%N_fol*(1-Sf/F)
! Summation, due to the filling of the pool at other points as well
p%coh%litC_fol = p%coh%litC_fol + p%coh%ntreea * Sf * cpart
p%coh%litC_frt = p%coh%litC_frt + p%coh%ntreea * Sr * cpart
! New version with species specific N content and reallocation factor (see species.par)
! changed to 1-reallo
p%coh%litN_fol = p%coh%litN_fol + p%coh%ntreea * Sf * cpart * (1.-spar(ns)%reallo_fol) / spar(ns)%cnr_fol
p%coh%litN_frt = p%coh%litN_frt + p%coh%ntreea * Sr * cpart * (1.-spar(ns)%reallo_frt) / spar(ns)%cnr_frt
! second step: allocation of N to new growth
! before bud-break allocation to leaves is 50% of the N content of last years foliage
tbc_root_Ndemand = Gs*cpart *kg_in_g / spar(ns)%cnr_tbc + Gr* cpart/spar(ns)%cnr_frt*kg_in_g
IF(tbc_root_Ndemand + Gf*p%coh%med_sla*0.5 > p%coh%N_pool) THEN
if (tbc_root_Ndemand .gt. 1E-8) then
Nredfak = AMAX1((p%coh%N_pool-Gf*p%coh%med_sla*0.5) / tbc_root_Ndemand,0.) ! Div. by zero possible !
else
Nredfak = 0.
endif
tbc_root_Ndemand = tbc_root_Ndemand*Nredfak
ENDIF
p%coh%N_pool = p%coh%N_pool - tbc_root_Ndemand
END SUBROUTINE PARTITION_SV
!*******************************!
!* SUBROUTINE PARTITION_MI *!
!*******************************!
SUBROUTINE PARTITION_MI( p )
!*** Declaration part ***!
USE data_par
USE data_stand
USE data_simul
IMPLICIT NONE
TYPE(Coh_Obj) :: p ! pointer to cohort list
!no partitioning, foliage mass keeps constant
p%coh%x_fol = p%coh%x_fol ! !FNew
p%coh%x_sap = 0.!SNew
p%coh%x_frt = 0.!RNew
END SUBROUTINE PARTITION_MI
!***************************!
! FUNCTION f_lf *!
!***************************!
REAL FUNCTION f_lf(x)
USE data_stand
USE data_plant
REAL :: x
f_lf = ksi*x**kappa + x - M_avail
END ! FUNCTION f_lf
!***************************!
! FUNCTION df_lf *!
!***************************!
REAL FUNCTION df_lf(x)
USE data_stand
USE data_plant
REAL :: x
df_lf = ksi*kappa*x**(kappa-1.) + 1.
END ! FUNCTION df_lf
!***************************!
! FUNCTION ddf_lf *!
!***************************!
REAL FUNCTION ddf_lf(x)
USE data_stand
USE data_plant
REAL :: x
ddf_lf = ksi*kappa*(kappa-1.)*x**(kappa-2.)
END ! FUNCTION ddf_lf
!***************************!
! SUBROUTINE CALC_DBH *!
!***************************!
SUBROUTINE CALC_DBH(B, Ht, S, H, Ahb, Ahc, ident, dbh, dc, hs, Asw)
!*** Declaration part ***!
USE data_par
USE data_species
USE data_simul
IMPLICIT NONE
INTEGER :: ident
REAL :: Dc ! diameter at crown base
REAL :: B, & ! bole height,
Ht, & ! total tree height
S, & ! sapwood
H, & ! heartwood
hs, & ! sapwood height
D, & ! stem diameter at forest floor
DBH, & ! tree diameter at breast height
Ahb, & ! cross sectional area heartwood at tree base
Ahc, & ! cross sectional area of heartwood at crown base
Asw, & ! cross sectional area of sapwood in bole
discr, func, help, hp1, hp2,hp3, hp4
REAL :: fp, fq, & ! coefficients of quadratic equation
w1, w2, & ! solutions of quadratic equation
precision ! criterion for acceptance of solution
real :: sprhos ! sapwood density [kg/cm3]
!*** Calculation part ***!
precision = 1.e-5
sprhos = spar(ns)%prhos
! calculate Diameters
hs = (2*B +Ht)/3.
Asw = S/(spar(ns)%prhos*hs)
! if Bole height >= height trees are dead and calculations not required
IF(B .lt. Ht) THEN
select case (flag_volfunc)
case (0)
D = SQRT( (S+H)*4. / (sprhos*hs*pi) )
IF (Ht<h_breast) THEN
DBH = 0.0
ELSEIF (Ht>h_breast.and.B<h_breast) then
DBH=D-(D/(Ht-B))*(h_breast-B)
ELSE
DBH=D
ENDIF
case (1)
D = SQRT((Ahb+Asw)*4./pi)
! if Bole height = 0 then there is no need to calulate Diameter at crown base and Dc = D
IF(B.EQ.0.) THEN
Dc = D
ELSE
fp = -2. * (B/Ht) * (3.*H/(sprhos*B)-Ahb)-Ahb*(B/Ht)**2.
fp = -2. * B/Ht * (3.*H/(sprhos*B)-Ahb)-Ahb*(B/Ht)**2.
fq = ((3.*H/(sprhos) - Ahb*B) / Ht)**2.
discr = fp**2./4.-fq
if (abs(discr) .lt. zero) then
discr = zero ! avoid small values
endif
! No solution
if(discr.lt.0) then
if (.not.flag_mult8910) then
CALL error_mess(time,'discriminant < 0 in calc_dbh for cohort: ',real(ident))
CALL stop_mess(time,'discriminant < 0 in calc_dbh ')
CALL error_mess(time,'stop in calc_dbh for stand No: ',real(ip))
CALL error_mess(time,'heart wood mass H: ',H)
CALL error_mess(time,'bole height b: ',b)
CALL error_mess(time,'height Ht: ',Ht)
CALL error_mess(time,'ave. sapwood height hs: ',hs)
CALL error_mess(time,'sapwood area Asw: ',Asw)
CALL error_mess(time,'heartwood area at stem base Ahb: ',Ahb)
endif
flag_end = 1
return
end if
discr = SQRT(discr)
w1 = -fp/2. + discr
w2 = -fp/2. - discr
1313 hp1 = SQRT(w1*Ahb)
hp2 = (Ahb+SQRT(w1*Ahb))*B
hp3 = (w1*Ht + (Ahb+SQRT(w1*Ahb))*B)
help = (sprhos/3.) * (w1*Ht + (Ahb+SQRT(w1*Ahb))*B)
func = (sprhos/3.) * (w1*Ht + (Ahb+SQRT(w1*Ahb))*B) - H
hp4= H* precision
IF(abs(func) <= H * precision) THEN
Ahc = w1
if (.not.flag_mult8910) then
CALL error_mess(time,' positive root is a solution in calc_dbh for cohort: ',real(ident))
CALL error_mess(time,'stop in calc_dbh for stand No: ',real(ip))
CALL error_mess(time,'function: ',func)
endif
flag_end = 1
return
ELSE
func = (sprhos/3.) * (w2*Ht + (Ahb+SQRT(w2*Ahb))*B) - H
IF(abs(func) <= H * precision) THEN
Ahc = w2
ELSE
IF(precision.LT.1e-2) THEN
precision = precision*10.
GOTO 1313
if (.not.flag_mult8910) then
CALL error_mess(time,'no valid solution found in calc_dbh for heartwood geometry for cohort: ',real(ident))
CALL error_mess(time,': heart wood mass, H = ',H)
CALL error_mess(time,': precision requirement = ',precision)
CALL error_mess(time,'iteration in stand No: ',real(ip))
endif
ELSE
if (.not.flag_mult8910) then
CALL error_mess(time,'no valid solution found in calc_dbh for heartwood geometry for cohort: ',real(ident))
CALL stop_mess(time,'no valid solution found in calc_dbh for heartwood geometry')
CALL error_mess(time,'species No: ',real(ns))
CALL error_mess(time,'stop in calc_dbh for stand No: ',real(ip))
CALL error_mess(time,'precision requirement H*precision ',H*precision)
CALL error_mess(time,'heart wood mass H: ',H)
CALL error_mess(time,'bole height b: ',b)
CALL error_mess(time,'height Ht: ',Ht)
CALL error_mess(time,'ave. sapwood height hs: ',hs)
CALL error_mess(time,'sapwood area Asw: ',Asw)
CALL error_mess(time,'heartwood area at stem base Ahb: ',Ahb)
endif
flag_end = 1
return
ENDIF
ENDIF
ENDIF
Dc = SQRT((Ahc+Asw)*4./pi)
END IF
if (Ht<=h_breast) then
DBH = 0.0
else if (Ht>h_breast.and.B<h_breast) then
DBH=Dc*(Ht-h_breast)/(Ht-B)
else
DBH=D-(D-Dc)*h_breast/B
end if
end select
ELSE
if (.not.flag_mult8910) then
CALL error_mess(time,'no calculation of heartwood geometry for cohort (Bole height >= height trees are dead): ',real(ident))
CALL error_mess(time,'bole height: ',b)
CALL error_mess(time,'height: ',Ht)
endif
END IF ! if B > Ht
END SUBROUTINE CALC_DBH