!*****************************************************************!
!*                                                               *!
!*              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
             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
