!! Compute plants mortality 
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  1.0 - 5th February 2020 
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 5/Feb/2020 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Routines to compute plants mortality 
! 
MODULE PlantsMortality


! Modules used:

USE DataTypeSizes, ONLY : &
   ! Imported Type Definitions:
   short, float

USE Units, ONLY : &
    ! imported parameters:
    hectare, pi, year

USE PlantsAllometrics, ONLY : &
    !imported routines:
    CanopyCover

IMPLICIT NONE

!global routines:
PUBLIC :: KillPlants

!==============
    CONTAINS
!==============


    
!==============================================================================
!| Description:
! the number of trees is updated considering the mortality of the plants.
! for each tree that dies, a fraction mi of the mean biomass
! three types of mortality are considered: 
!
! 1. The first mortality is due to the self thinning (the one included in 3PG), 
!     which basically ensures that the mean single-tree stem biomass WS does 
!     not exceed the maximum permissible single-tree stem biomass 
!     WSx (kg·trees-1) (Sands, 2004). We assume the maximum target density
!     is reached after 1 year
!
! 2. The second mortality is age dependent mortality following the approach 
!     of `LPJ-GUESS` (Smith et al., 2014), with aging the plants become 
!     more suceptible to the wind, diseases, etc.
!
! 3. the third mortality is the so called the "crowding competition function", 
!     this mortality insures that the % of cover of pixel do not exceed the 95% 
!
!  References:
!  
!  Sands P., “Adaptation of 3-PG to novel species: guidelines for data 
!    collection and parameter assignment”, Cooperative Research 
!    Center for Sustainable Production Forestry, 2004
!
!  Smith, B., Wårlind, D., Arneth, A., Hickler, T., Leadley, P., 
!      Siltberg, J., and Zaehle, S.: Implications of incorporating N 
!      cycling and N limitations on primary production in an individual-based 
!      dynamic vegetation model, Biogeosciences, 11, 2027–2054,
!      https://doi.org/10.5194/bg-11-2027-2014, 2014.
!
SUBROUTINE  KillPlants &
!
( dt, age, agemax, ms, mf, mr, wSx1000 , dc, density,  cc, ws, wr, wf, wtot ) 
    
IMPLICIT NONE

!arguments with intent(in):
INTEGER (KIND = short), INTENT(IN) :: dt !! time step (s)
REAL (KIND = float), INTENT(IN) :: age !! current tree age (year)
REAL (KIND = float), INTENT(IN) :: agemax !! maximum age (year)
REAL (KIND = float), INTENT(IN) :: ms !! Fraction of mean stem biomass pools per tree on each dying tree 
REAL (KIND = float), INTENT(IN) :: mf !! Fraction of mean follioge biomass pools per tree on each dying tree 
REAL (KIND = float), INTENT(IN) :: mr !!  Fraction of mean roots biomass pools per tree on each dying tree
REAL (KIND = float), INTENT(IN) :: wSx1000 !! maximum permissible single-tree stem biomass  (t) 
REAL (KIND = float), INTENT(IN) :: dc !! crown diameter (m)

!arguments with intent(inout):
REAL (KIND = float), INTENT(INOUT) :: density !! number of plants per hectar
REAL (KIND = float), INTENT(INOUT) :: cc !! canopy cover
REAL (KIND = float), INTENT(INOUT) :: ws !! stem biomass (t/ha)
REAL (KIND = float), INTENT(INOUT) :: wr !! roots biomass (t/ha)
REAL (KIND = float), INTENT(INOUT) :: wf !! folliage biomass (t/ha)
REAL (KIND = float), INTENT(INOUT) :: wtot !! total biomass (t/ha)



REAL (KIND = float) :: density_after !! the updated number of the trees after considering the mortality 
REAL (KIND = float) :: CCmax = 0.95 !! maximum target canopy cover 
REAL (KIND = float) :: CCred !!canopy cover reduction within time step

! Fractions of mean foliage, root and stem biomass pools per tree on each dying tree 
! 1- the first mortality, the need for this mortality should be checked at the end of each time step 
!local declarations:

!INTEGER (KIND = short) :: i
REAL (KIND = float) :: n
!REAL (KIND = float) :: X1
!REAL (KIND = float) :: X2
!REAL (KIND = float) :: dfN
!REAL (KIND = float) :: dN
!REAL (KIND = float) :: fN
REAL (KIND = float) :: Dthinning !! the trees number reduction after the first type of mortality
REAL (KIND = float) :: Dselfthinning !! probabilistic function for age-dependent Mortality
REAL (KIND = float) :: Dage !! probabilistic function for age-dependent Mortality
! variables required for the update of the biomass 
REAL (KIND = float) :: delStems !! the number of dead trees 

!DEBUG
REAL (KIND = float) :: wSmax
REAL (KIND = float) :: thinPower = 3./2.
REAL (KIND = float) :: avStemMass

!---------------------------end of declarations--------------------------------


! 1- the first type of mortality

!max stem mass
wSmax = wSx1000 * (1000 / density) ** thinPower 

!current average stem mass
avStemMass = ws / density


IF ( wSmax < avStemMass ) THEN !compute self thinning
    
    !compte maximum target density
    n = 1000. / ( avStemMass / wSx1000 ) ** (2./3.)
    
    !compute number of trees to kill
    IF (dt <  10 * year) THEN
        Dthinning = (density - n) * dt /  ( 10 * year) 
    ELSE
         Dthinning = density - n
    END IF

ELSE
        Dthinning = 0.
END IF
Dthinning = 0. !DEBUG tolgo questa mortalità, da indagare

 
!  2- the second mortality 
Dage = - 3. * lOG10 (0.001/agemax) * (age/agemax)**2  * dt / ( 10 * year)  !DEBUG lascio solo la mortalità per età


! 3- the third mortality 
IF ( cc <= CCmax) THEN !low density
    Dselfthinning = 0.
ELSE ! high density
     CCred = (cc - CCmax) * dt / ( 30 * year )
     n = 4. * (cc - CCred) * hectare / (Pi * dc**2.) 
     Dselfthinning = (density - n )
END IF
Dselfthinning = 0.  !DEBUG tolgo questa mortalità, da indagare


! UPDATE  the number of the trees after the mortality 
Density_after = density - ( density * DAGE + Dselfthinning + Dthinning )


! update the biomass pools according to the mortality of the trees

delStems= density - Density_after
wf = wf - mf * delStems * (WF / density)
wr = wr - mr * delStems * (WR / density)
ws = ws - ms * delStems * (WS / density)
wtot = wf + wr + ws

!update density
density = density_after

!update canopy cover
cc = CanopyCover (dc, density)

RETURN
END SUBROUTINE KillPlants

    
END MODULE PlantsMortality