!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* - stand_mort *!
!* - int_mort intrinsic mortality rate *!
!* - stress_mort stress mortality rate *!
!* - int_mort_weib *!
!* *!
!* - Calculation of dead trees per cohort and species *!
!* deterministic approach *!
!* - relative mortality rate is determined by intr. mortality *!
!* and stress mortality *!
!* - stress mortality is calculated depending on *!
!* npp, ystress, yhealth *!
!* - intrinsic probability is optionally calculated on *!
!* age of cohort *!
!* - for each tree of the cohort mortality is decided *!
!* by means of the Mortality probability and *!
!* a uniformly distributed variable *!
!* *!
!* 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 *!
!* *!
!*****************************************************************!
! input variables:
! pro cohort NPP
! state variables:
! pro cohort nTreeA nTreeD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE stand_mort
USE data_stand
USE data_species
USE data_simul
USE data_manag
use data_out
use data_par
IMPLICIT NONE
!local variables
INTEGER :: flag_hgrowth
INTEGER :: taxnr
REAL :: intmort
REAL :: strmort
REAL :: totmort
REAL :: totmort_m
REAL :: besmort
REAL :: ntdead
REAL :: ntdead_m
REAL :: nhelp
REAL :: survpfunct
INTEGER :: hage
REAL :: besp1,besp2 ! parameters for besetting mortality
real :: help1, help2
real :: intmorth
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' stand_mort'
if (flag_standup .eq. 0) flag_standup = 1 ! call stand_balance later
ntdead=0.
nhelp=0.
ntdead_m=0
flag_hgrowth=0
sumvsdead = 0.
sumvsdead_m3 = 0.
besmort = 0
strmort = 0.
totmort = 0.
svar%sumvsdead = 0
svar%sumvsdead_m3 = 0.
zeig=>pt%first
DO
IF(.not.associated(zeig)) exit
IF (zeig%coh%height.ge.thr_height.and. zeig%coh%species.le.nspec_tree) then
taxnr=zeig%coh%species
IF(time.eq.1) then
zeig%coh%nta=zeig%coh%nTreeA
ELSE
IF (flag_mg.eq.1) then
IF (thin_year(act_thin_year-1).eq.(time-1)) zeig%coh%nta=zeig%coh%nTreeA
ELSE IF (flag_mg.eq.2) then
if(flag_adapm .eq. 1) zeig%coh%nta=zeig%coh%nTreeA
ENDIF
ENDIF
IF(zeig%coh%notViable) then
print*,time, zeig%coh%notViable
zeig%coh%nTreeD = zeig%coh%ntreeA
zeig%coh%nta = 0.
zeig%coh%ntreeA = 0
goto 1000
ENDIF
! calculation of stress and health indicator based on foliage biomass increment
hage = zeig%coh%x_age
IF (flag_hgrowth==0) THEN
IF(zeig%coh%fol_inc.le.0.0) then
zeig%coh%x_stress = zeig%coh%x_stress + 1
zeig%coh%x_health= 0
ELSE
zeig%coh%x_health = zeig%coh%x_health + 1
IF(zeig%coh%x_stress.eq.1.and. zeig%coh%x_health.gt.0) zeig%coh%nta = zeig%coh%ntreeA
IF(zeig%coh%x_stress .gt.0) zeig%coh%x_stress = zeig%coh%x_stress - 1
ENDIF
ENDIF
IF (flag_hgrowth==1) THEN
IF(zeig%coh%bio_inc.le.0.0) then
zeig%coh%x_stress = zeig%coh%x_stress + 1
zeig%coh%x_health= 0
ELSE
zeig%coh%x_health = zeig%coh%x_health + 1
IF(zeig%coh%x_stress.eq.1.and. zeig%coh%x_health.gt.0) zeig%coh%nta = zeig%coh%ntreeA
IF(zeig%coh%x_stress .gt.0) zeig%coh%x_stress = MAX(zeig%coh%x_stress - 3,0)
ENDIF
ENDIF
! calculation of relative mortality rates
! intrinsic mortality
! constant
call int_mort(taxnr,intmorth)
! age-depending using Weibull function
call int_mort_weib(taxnr,intmort,hage)
! stress mortality
IF(zeig%coh%x_stress.gt.0) then
IF(flag_hgrowth==0) strmort = weibal*spar(taxnr)%weibla*zeig%coh%x_stress**(weibal-1)
IF(flag_hgrowth==1) strmort = weibal*spar(taxnr)%weibla*(zeig%coh%x_stress/3.)**(weibal-1)
ELSE
strmort = 0.
ENDIF
!mortality depending on gross growth rate foliage
IF(strmort==0.0 .AND. flag_hgrowth==2) THEN
IF(zeig%coh%sfol/zeig%coh%gfol.GT.0.9) THEN
strmort=((zeig%coh%sfol/zeig%coh%gfol-0.9)*10.)**2
ELSE
ENDIF
ENDIF
if(strmort==0. .and. flag_hgrowth==3) then
help1 = 10**((log10(4.5)-log10((zeig%coh%x_sap + zeig%coh%x_hrt))*zeig%coh%ntreea)/1.5)
help2 = help1/zeig%coh%ntreea
end if
! mortality caused by besetting for oak
besp1= 0.018
besp2= 0.0216
if(zeig%coh%species.eq.4) then
if( zeig%coh%bes.le. 1.2) then
besmort = 0.
else
besmort = besp1*zeig%coh%bes- besp2
end if
else if (zeig%coh%species.eq.1.) then
if(zeig%coh%bes.le. 1.2) then
besmort = 0.
else
besmort = 0.04*zeig%coh%bes- 0.04
end if
end if
!total mortality rate depending on flag_mort
IF(flag_mort.eq.1) THEN
totmort = strmort
ELSE IF(zeig%coh%x_age.le.30) then
totmort=intmort+(1-intmort)*strmort
if(taxnr.eq.8) totmort = strmort
ELSE
totmort=intmort+(1-intmort)*strmort
ENDIF
! if species type oak then combination of stress mortality and besetting mortality
if(zeig%coh%species.eq.4.or.zeig%coh%species.eq.1) then
totmort = besmort + (1-besmort)* strmort
end if
survpfunct = exp(- spar(taxnr)%weibla * zeig%coh%x_stress**weibal)
ntdead = totmort*zeig%coh%nTreeA
IF(totmort.GT.1) CALL error_mess(time,"totmort greater 1: ",totmort)
! calculation of real stem number
zeig%coh%nta = zeig%coh%nta - ntdead
! calculation of integer stem number
zeig%coh%nTreeD = zeig%coh%nTreeA-NINT(zeig%coh%nta)
zeig%coh%nTreeA = NINT(zeig%coh%nta)
IF(zeig%coh%nTreeA.lt.1.) zeig%coh%nTreeA=0.
if (zeig%coh%mistletoe.eq.1) then ! in case Mist.infect. tree dies, number of mistletoes dies, too
totmort_m = zeig%coh%nTreeD/(zeig%coh%nTreeD+zeig%coh%nTreeA) ! share of trees removed of total trees. used for the share of mistletoe that dies
ntdead_m = 1. !flag
end if
! calculation of the litter pool of a died tree of a cohort
1000 IF (zeig%coh%ntreeD.ne.0) then
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreeD*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreeD*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreeD*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreeD*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreeD*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreeD*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreeD*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreeD*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
if(flag_mg.ne.0) then
! if(thin_dead.eq.0.and.thin_flag1(1).lt.0.) then
if(thin_dead.eq.0) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
end if
else if(zeig%coh%diam.le.tardiam_dstem.or. flag_mg.eq.0) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
else if(zeig%coh%diam.gt.tardiam_dstem.and.flag_mg.ne.0.or.flag_mg.eq.5) then
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
else if(thin_dead.eq.1.and.zeig%coh%diam.le.tardiam_dstem) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
end if
ENDIF
ENDIF
zeig=>zeig%next
ENDDO
! if tree cohort with mistletoe changed, change number of mistletoes too
if (ntdead_m.eq.1.) then
zeig => pt%first
do while (associated(zeig))
if (zeig%coh%species.eq.nspec_tree+2) then
zeig%coh%nta=zeig%coh%nTreeA
ntdead_m = totmort_m*zeig%coh%nTreeA
zeig%coh%nta = zeig%coh%nta - ntdead_m
zeig%coh%nTreeD = zeig%coh%nTreeA-NINT(zeig%coh%nta)
zeig%coh%nTreeA = NINT(zeig%coh%nta)
IF(zeig%coh%nTreeA.lt.1.) then
zeig%coh%nTreeA=0.
flag_mistle=0 !set flag mistletoe back to zero
ENDIF
endif
zeig=>zeig%next
enddo ! zeig (cohorts)
end if
ntdead_m=0.
! recalculation sumvsdead
sumvsdead = sumvsdead * 10000./kpatchsize ! kg/patch ---> ! kg/ha
sumvsdead_m3 = sumvsdead_m3 * 10000./kpatchsize ! kg/patch ---> ! kg/ha
cumsumvsdead = cumsumvsdead + sumvsdead
END SUBROUTINE stand_mort
SUBROUTINE int_mort(taxnr,intmort)
USE data_species
IMPLICIT NONE
REAL :: intmort
INTEGER :: taxnr
intmort=1.-exp(-spar(taxnr)%intr)
END SUBROUTINE int_mort
SUBROUTINE int_mort_weib(taxnr,intmort,hage)
USE data_species
USE data_stand
USE data_simul
IMPLICIT NONE
REAL :: intmort, weibla_int
INTEGER :: taxnr
INTEGER :: hage
! Weibull functions depending on age
weibla_int = -log(0.01)/(spar(taxnr)%max_age**weibal_int)
intmort = weibal_int*weibla_int*(hage)**(weibal_int-1.)
END SUBROUTINE int_mort_weib