!! Simulate vegetation dynamic
!|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 - 31st October 2018  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 31/Oct/2018 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Module to model vegetation dynamic
! 
MODULE Plants

! 
! Modules used:

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

USE LogLib, ONLY: &
   !Imported routines:
   Catch

USE GridLib, ONLY: &
    !imported definitions:
    grid_real, grid_integer, &
    !imported routines:
    NewGrid, GridDestroy, &
    !imported parameters:
    NET_CDF

USE GridOperations, ONLY : &
   !imported routines:
   GridByIni, CellArea

USE Chronos, ONLY: &
    !imported definitions:
    DateTime, &
    !Imported routines:
    DayOfYear, GetYear, &
    !Imported operations:
    OPERATOR (==), &
    OPERATOR ( > )

USE IniLib, ONLY : &
    !Imported derived types:
    IniList, &
    !Imported routines:
    IniOpen, IniClose , &
    IniReadInt , IniReadReal, &
    IniReadString, SectionIsPresent, &
    KeyIsPresent

USE StringManipulation, ONLY : &
    !imported routines
    ToString

USE PlantsAllometrics, ONLY : &
    !imported routines:
    CrownDiameter, CanopyCover

USE PlantsModifiers, ONLY: &
    !imported routines:
    AGEmod, TEMPmod, SWCmod, VPDmod, CO2mod

USE PlantsInterception, ONLY: &
    !imported variables
    interceptionParametersByMap, &
    laimax, canopymax

USE PlantsManagement, ONLY: &
    !imported definition:
    Practice, &
    !imported variables:
    plants_management, &
    management_map, &
    !imported routines:
    SetPlantsManagement, &
    ApplyPlantsManagement, &
    SetPractice

USE PlantsMortality, ONLY : &
    !imported routines:
    KillPlants

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

IMPLICIT NONE

! Global (i.e. public) Declarations:
PUBLIC :: PlantsConfig
PUBLIC :: PlantsGrow
PUBLIC :: PlantsParameterUpdate

INTEGER (KIND = short) :: dtPlants

TYPE (grid_real) :: lai !!leaf area index (m2/m2)
TYPE (grid_real) :: fvcover !!fraction of cell covered by vegetation
TYPE (grid_real) :: rsMin !!minimum stomatal resistance
TYPE (grid_real) :: gpp !! carbon gross primary production (t)
TYPE (grid_real) :: npp !! carbon net primary production (t)
TYPE (grid_real) :: carbonroot !! carbon mass of root (t)
TYPE (grid_real) :: carbonstem !! carbon mass of stem (t)
TYPE (grid_real) :: carbonleaf !! carbon mass of leaf (t)
TYPE (grid_real) :: dbh !!diameter at brest height (cm)
TYPE (grid_real) :: plantsHeight !! tree height (m)
TYPE (grid_real) :: density !! Trees/hectare
TYPE (grid_real) :: stemyield !! stem yield (t)

LOGICAL          :: rsMinLoaded = .FALSE.
LOGICAL          :: fvcoverLoaded = .FALSE.
LOGICAL          :: laiLoaded = .FALSE.
LOGICAL          :: plantsHeightLoaded = .FALSE.
LOGICAL          :: updatePlantsParameters = .FALSE.


TYPE PlantsSpecies
    !Descriptor
    CHARACTER (LEN = 50)  :: name  !!species name, Fagus, Abies, ...
    CHARACTER (LEN = 1)   :: phenology !! evergreen (E), deciduous (D)
    REAL (KIND = float)   :: alpha !! canopy quantum efficiency (molC/molPAR)
    REAL (KIND = float)   :: GPPtoNPP !!GPP/NPP ratio
    REAL (KIND = float)   :: k !! extinction coefficient for absorption of PAR by canopy
    REAL (KIND = float)   :: cra !! Chapman-Richards asymptotic maximum height
    REAL (KIND = float)   :: crb !! Chapman-Richards exponential decay parameter
    REAL (KIND = float)   :: crc !! Chapman-Richards shape parameter
    REAL (KIND = float)   :: as !! scaling coefficient in stem mass v diameter relationship
    REAL (KIND = float)   :: ns !! scaling exponent in stem mass v diameter relationship
    !REAL (KIND = float)   :: saf !! stem allocation factor
    !REAL (KIND = float)   :: raf !! root allocation factor
    !REAL (KIND = float)   :: laf !! leaf (foliage) allocation factor
    REAL (KIND = float)   :: dbhdcmin !! minimum ratio between stem and crown diameters (m/cm)
    REAL (KIND = float)   :: dbhdcmax !! maximum ratio between stem and crown diameters (m/cm)
    REAL (KIND = float)   :: denmin !! minimum tree density (trees/ha)
    REAL (KIND = float)   :: denmax !! minimum tree density (trees/ha)
    REAL (KIND = float)   :: agemax !! maximum age
    REAL (KIND = float)   :: tmax !! maximum temperature for vegetation growing (°C)
    REAL (KIND = float)   :: tmin !! minimum temperature for vegetation growing (°C)
    REAL (KIND = float)   :: topt !! optimum temperature for vegetation growing (°C)
    REAL (KIND = float)   :: theta_fswc !! parameter to compute soil water content modifier 
    REAL (KIND = float)   :: theta_fvpd !! parameter to compute vapor pressure deficit modifier
    REAL (KIND = float)   :: fpra !! parameter 1 to compute allocation factors
    REAL (KIND = float)   :: fprn !! parameter 2 to compute allocation factors
    REAL (KIND = float)   :: spra !! parameter 3 to compute allocation factors
    REAL (KIND = float)   :: sprn !! parameter 4 to compute allocation factors
    REAL (KIND = float)   :: ltr !!leaf turnover rate [s-1]
    REAL (KIND = float)   :: rtr !!root turnover rate [s-1]
    REAL (KIND = float)   :: tcold_leaf !! temperature threshold that accelerates leaf turnover (°C)
    REAL (KIND = float)   :: sla !! specific leaf area  [m2/Kg]
    REAL (KIND = float)   :: hdmin !! H/D ratio in carbon partitioning for low density
    REAL (KIND = float)   :: hdmax !! H/D ratio in carbon partitioning for high density
    REAL (KIND = float)   :: albedo !!plant albedo (0-1)
    REAL (KIND = float)   :: laimax !!maximum leaf area index used for precipitation interception (m2/m2)
    REAL (KIND = float)   :: canopymax !! maximum canopy storage capacity (m)
    REAL (KIND = float)   :: wood_density !! wood density (kg/m3)
    REAL (KIND = float)   :: ms !! Fractions of mean stem biomass pools per tree on each dying tree
    REAL (KIND = float)   :: mr !!  Fractions of mean root biomass pools per tree on each dying tree
    REAL (KIND = float)   :: mf !!  Fractions of mean leaf biomass pools per tree on each dying tree
    REAL (KIND = float)   :: wSx1000 !!  Fractions of mean leaf biomass pools per tree on each dying tree
END TYPE PlantsSpecies

TYPE PlantsCohort !group of plants of the same species that are born during the same year
    TYPE (PlantsSpecies) :: species
    !state variables
    REAL (KIND = float) :: age !!age [years]
    REAL (KIND = float) :: gpp !!  gross primary production (t/ha)
    REAL (KIND = float) :: npp !!  net primary production (t/ha)
    REAL (KIND = float) :: lai !! leaf area index [m2/m2]
    REAL (KIND = float) :: mass_root !!  mass of root (t/ha)
    REAL (KIND = float) :: mass_stem !!  mass of stem (t/ha)
    REAL (KIND = float) :: mass_stem_year_previous !! mass of stem of the previous year (t/ha)
    REAL (KIND = float) :: mass_leaf !!  mass of leaf (t/ha)
    REAL (KIND = float) :: mass_total !!  total biomass root + stem + foliage (t/ha)
    REAL (KIND = float) :: stem_yield !!  stem mass produced by plants cut (t/ha)
    REAL (KIND = float) :: density !! number of plants per hectar
    REAL (KIND = float) :: height !! plants heigth [m]
    REAL (KIND = float) :: dbh !!diameter at brest heigth [cm]
    REAL (KIND = float) :: canopy_cover !! percentage of surface stand covered by canopy [0-1]
    REAL (KIND = float) :: crown_diameter !! diameter of crown [m]
    REAL (KIND = float) :: apar !! absorbed photosynthetically active radiation (molPAR m-2)
    TYPE(PlantsCohort), POINTER :: next =>null() ! pointer to the next cohort in the same stand
END TYPE PlantsCohort


TYPE PlantsStand ! list of cohorts
    TYPE (PlantsCohort), POINTER, PRIVATE :: first => null()
    !TYPE (PlantsCohort), POINTER, PRIVATE :: last => null()
    !TYPE (PlantsCohort), POINTER, PRIVATE :: iter => null()
    INTEGER (KIND = short)      , PRIVATE :: lenght = 0
    INTEGER (KIND = short)      , PRIVATE :: i !row of the stand in grid plants_mask
    INTEGER (KIND = short)      , PRIVATE :: j !column of the stand in grid plants_mask
     TYPE (Practice)             :: thinning
END TYPE PlantsStand

!TYPE PlantsAge
!    REAL (KIND = float) :: value  !!age [yesrs]
!    INTEGER (KIND = short) :: period !! 0 for adult tree, 1 for young tree
!    INTEGER (KIND = short) :: species_count !! number of species in the cell
!    TYPE (PlantsSpecies), ALLOCATABLE :: species (:)
!END TYPE PlantsAge
!
!TYPE PlantsHeight
!    REAL (KIND = float) :: value  !!value of height class [m]
!    REAL (KIND = float) :: layer_coverage !! cell coverage of that layer
!    INTEGER (KIND = short) :: dominance !! -1 no trees in veg period, 1 trees in veg period
!    INTEGER (KIND = short) :: ages_count !! counter of age classes in the cell
!    TYPE (PlantsAge), ALLOCATABLE :: ages (:)
!END TYPE PlantsHeight




!local (i.e. private) declarations:

!variables and parameters
LOGICAL, PRIVATE :: simulatePlants  != 1 when plants dynamic is simulated, = 0 read vegetation parameters from files
LOGICAL, PRIVATE :: useCO2modifier ! define wheter use the CO2 modifier
LOGICAL, PRIVATE :: mortality ! wheter to apply mortality or not
TYPE (grid_integer), PRIVATE :: plants_mask !!define cells active for plants dynamic simulation
TYPE (PlantsStand), PRIVATE, ALLOCATABLE :: forest (:)  !vector of stands
INTEGER (KIND = short), PRIVATE :: count_stands  !number of stands in the forest (equal to the number of cells in the domain)
TYPE (PlantsSpecies), ALLOCATABLE, PRIVATE :: species (:) !
!REAL (KIND = float), PRIVATE :: rs_par !! fraction of shortwave radiation to compute PAR [0-1]
INTEGER (KIND = short), PRIVATE :: year_prev, year_new
REAL (KIND = float), PRIVATE, PARAMETER :: C_molar_mass = 0.0120107 ! kg/mol

!routines
PRIVATE :: AparCalc
PRIVATE :: ReadSpecies
PRIVATE :: SetIC
PRIVATE :: GrowLeaf
PRIVATE :: GrowDBH
PRIVATE :: ForestToGrid


!=======
CONTAINS
!=======
    
!==============================================================================
!| Description:
!  Initialize Plants module
SUBROUTINE  PlantsConfig &
!
(iniFile, mask, begin, end)


IMPLICIT NONE

! Arguments with intent(in):
CHARACTER (LEN = 300), INTENT(IN) :: iniFile !!configuration file 
TYPE (grid_integer), INTENT(IN) :: mask  !!mask of simulation domain
TYPE (DateTime), INTENT (IN) :: begin !!simulation starting date
TYPE (DateTime), INTENT(IN) :: end !!simulation ending date


!local declarations:
TYPE(IniList) :: iniDB !!store configuration info
INTEGER (KIND = short) :: i, j, k
CHARACTER (LEN = 300) :: species_file
CHARACTER (LEN = 300) :: ICfile
CHARACTER (LEN = 300) :: PMfile
REAL (KIND = float)   :: scalar


!------------end of declaration------------------------------------------------    

CALL IniOpen (iniFile, iniDB)

!need plants dynamic simulation?
simulatePlants = IniReadInt ('plants-simulation', iniDB)

IF ( simulatePlants == 1) THEN !plants dynamic is simulated: configure forest model
    !plants mask 
    IF ( SectionIsPresent ( 'plants-mask', iniDB)  ) THEN
        CALL GridByIni (iniDB, plants_mask, section = 'plants-mask')
    
    ELSE
        CALL Catch ('error', 'Plants', 'Mask missing in configuration file')
    END IF
    
    !CO2 modifier
    IF ( SectionIsPresent ( 'use-co2-modifier', iniDB)  ) THEN
        IF ( IniReadInt ('use-co2-modifier', iniDB) == 1 ) THEN
            useCO2modifier = .TRUE.
        ELSE 
            useCO2modifier = .FALSE.
        END IF
    ELSE
       !if key is not found, set it to default false
       useCO2modifier = .FALSE.
    END IF
    
    
    !mortality
    IF ( KeyIsPresent ( 'mortality', iniDB)  ) THEN
        IF ( IniReadInt ('mortality', iniDB) == 1 ) THEN
            mortality = .TRUE.
        ELSE 
            mortality = .FALSE.
        END IF
    ELSE
       !if mortality is not found, assume it is not considered
       mortality = .FALSE.
    END IF
    
    
    !load species
    species_file = IniReadString ('plants-species-file', iniDB)
    CALL ReadSpecies (species_file)
    
    !compute number of stands and allocate 
    count_stands = 0
    DO i = 1, plants_mask % idim
        DO j = 1, plants_mask % jdim
            IF ( plants_mask % mat (i,j) /= plants_mask % nodata ) THEN
               count_stands = count_stands + 1
            END IF
        END DO
    END DO

    ALLOCATE ( forest ( count_stands) )

    !initialize list
    k = 1
    DO i = 1, plants_mask % idim
        DO j = 1, plants_mask % jdim
            IF ( plants_mask % mat (i,j) /= plants_mask % nodata ) THEN
                !initialize the first cohort in the stand
                ALLOCATE ( forest (k) % first)
                NULLIFY ( forest (k) % first % next)
                forest (k) % lenght = 1
                !store position of stand
               forest (k) % i = i
               forest (k) % j = j
               k = k + 1
            END IF
        END DO
    END DO
    
    
    !set initial condition
    ICfile = IniReadString ('plants-ic-file', iniDB)
    CALL SetIC (ICfile)
    
    !set plants management options
    IF (KeyIsPresent ('plants-management-file',iniDB) ) THEN 
        plants_management = .TRUE.
        PMfile = IniReadString ('plants-management-file', iniDB)
        CALL SetPlantsManagement (PMfile, begin, end)
        !set management pratice for each stand
        DO k = 1, count_stands
            CALL SetPractice (  management_map % mat (forest (k) % i, forest (k) % j) , forest (k) % thinning )
        END DO

    ELSE !management is turned off
        plants_management = .FALSE.
    END IF
    
    
    !set interception maps if required
    IF (interceptionParametersByMap == 0 ) THEN
         CALL Catch ('warning', 'Plants', 'interception parameters are set only where plants map is not nodata')
         
         !allocate maps
         CALL NewGrid (laimax, plants_mask)	
         CALL NewGrid (canopymax, plants_mask)	
         
         DO k = 1, count_stands
             i = forest (k) % i
             j = forest (k) % j 
             laimax % mat (i,j) = forest (k) % first % species % laimax
             canopymax % mat (i,j) = forest (k) % first % species % canopymax
         END DO
        
    END IF
    
    
    !allocate grids for state variables
    CALL NewGrid (lai, plants_mask)
    CALL NewGrid (fvcover, plants_mask)
    CALL NewGrid (gpp, plants_mask)
    CALL NewGrid (npp, plants_mask)
    CALL NewGrid (carbonroot, plants_mask)
    CALL NewGrid (carbonstem, plants_mask)
    CALL NewGrid (carbonleaf, plants_mask)
    CALL NewGrid (dbh, plants_mask)
    CALL NewGrid (plantsHeight, plants_mask)
    CALL NewGrid (density, plants_mask)
    CALL NewGrid (stemyield, plants_mask)
    
    !fill in grids
    CALL ForestToGrid (lai, 'lai') 
    CALL ForestToGrid (fvcover, 'fv')
    CALL ForestToGrid (gpp, 'gpp')
    CALL ForestToGrid (npp, 'npp')
    CALL ForestToGrid (carbonroot, 'root')
    CALL ForestToGrid (carbonstem, 'stem')
    CALL ForestToGrid (carbonleaf, 'leaf')
    CALL ForestToGrid (dbh, 'dbh')
    CALL ForestToGrid (plantsHeight, 'height')
    CALL ForestToGrid (density, 'density')
    CALL ForestToGrid (stemyield, 'stemyield')
    
    fvcoverLoaded = .TRUE.
    laiLoaded = .TRUE.
    plantsHeightLoaded = .TRUE.
    
ELSE ! no plants dynamic simulation required: read parameters from files
    
    !leaf area index
    IF ( SectionIsPresent ( 'lai', iniDB)  ) THEN
        IF (KeyIsPresent ('scalar', iniDB, 'lai') ) THEN
             scalar = IniReadReal ('scalar', iniDB, 'lai')
             CALL NewGrid (lai, mask, scalar)
        ELSE
            !check if parameter may change in time
             IF (KeyIsPresent ('format', iniDB, 'lai') ) THEN
                 IF ( IniReadString ('format', iniDB, 'lai') == 'net-cdf' ) THEN
                       updatePlantsParameters = .TRUE.
                 END IF
             END IF
             !read map
             CALL GridByIni (iniDB, lai, section = 'lai')
            
        END IF
        laiLoaded = .TRUE.
    ELSE
        CALL Catch ('warning', 'Plants', 'LAI missing in configuration file')
    END IF

    !fraction of vegetation coverage
    IF ( SectionIsPresent ( 'vegetation-fraction', iniDB)  ) THEN
        IF (KeyIsPresent ('scalar', iniDB, 'vegetation-fraction') ) THEN
             scalar = IniReadReal ('scalar', iniDB, 'vegetation-fraction')
             CALL NewGrid (fvcover, mask, scalar)
        ELSE
             !check if parameter may change in time
             IF (KeyIsPresent ('format', iniDB, 'vegetation-fraction') ) THEN
                 IF ( IniReadString ('format', iniDB, 'vegetation-fraction') == 'net-cdf' ) THEN
                       updatePlantsParameters = .TRUE.
                 END IF
             END IF
             !read map
             CALL GridByIni (iniDB, fvcover, section = 'vegetation-fraction')
        END IF
        fvcoverLoaded = .TRUE.
    ELSE
        CALL Catch ('warning', 'Plants', 'vegetation-fraction is missing')
    END IF
    
    
    !plants height
    IF ( SectionIsPresent ( 'vegetation-height', iniDB)  ) THEN
        IF (KeyIsPresent ('scalar', iniDB, 'vegetation-height') ) THEN
             scalar = IniReadReal ('scalar', iniDB, 'vegetation-height')
             CALL NewGrid (plantsHeight, mask, scalar)
        ELSE
             !check if parameter may change in time
             IF (KeyIsPresent ('format', iniDB, 'vegetation-height') ) THEN
                 IF ( IniReadString ('format', iniDB, 'vegetation-height') == 'net-cdf' ) THEN
                       updatePlantsParameters = .TRUE.
                 END IF
             END IF
             !read map
             CALL GridByIni (iniDB, plantsHeight, section = 'vegetation-height')
        END IF
        plantsHeightLoaded = .TRUE.
    ELSE
        CALL Catch ('warning', 'Plants', 'vegetation-height is missing')
    END IF
    
END IF

!minimum stomatal resistance. this is not computed by forest model and it may be used for computing PET
IF ( SectionIsPresent ( 'min-stomatal-resistance', iniDB)  ) THEN
        IF (KeyIsPresent ('scalar', iniDB, 'min-stomatal-resistance') ) THEN
             scalar = IniReadReal ('scalar', iniDB, 'min-stomatal-resistance')
             CALL NewGrid (rsMin, mask, scalar)
        ELSE
             !check if parameter may change in time
             IF (KeyIsPresent ('format', iniDB, 'min-stomatal-resistance') ) THEN
                 IF ( IniReadString ('format', iniDB, 'min-stomatal-resistance') == 'net-cdf' ) THEN
                       updatePlantsParameters = .TRUE.
                 END IF
             END IF
             !read map
             CALL GridByIni (iniDB, rsMin, section = 'min-stomatal-resistance')
        END IF
        rsMinLoaded = .TRUE.
    ELSE
        CALL Catch ('warning', 'Plants', 'min-stomatal-resistance is missing')
    END IF


CALL IniClose (iniDB)

RETURN
END SUBROUTINE PlantsConfig







!==============================================================================
!| Description: 
!  update parameter map that change in time
SUBROUTINE PlantsParameterUpdate   & 
  !
  (time)       

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time

!local declarations:
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname

!------------------------------end of declarations-----------------------------
  
  
  
!leaf area index
IF (  time == lai % next_time ) THEN
   !destroy current grid
   filename = lai % file_name
   varname = lai % var_name
   CALL GridDestroy (lai)
   !read new grid in netcdf file
   CALL NewGrid (lai, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF

!fraction of vegetation coverage
IF (  time == fvcover % next_time ) THEN
   !destroy current grid
   filename = fvcover % file_name
   varname = fvcover % var_name
   CALL GridDestroy (fvcover)
   !read new grid in netcdf file
   CALL NewGrid (fvcover, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF


!plants height
IF (  time == plantsHeight % next_time ) THEN
   !destroy current grid
   filename = plantsHeight % file_name
   varname = plantsHeight % var_name
   CALL GridDestroy (plantsHeight)
   !read new grid in netcdf file
   CALL NewGrid (plantsHeight, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF


!minimum stomatal resistance
IF (  time == rsMin % next_time ) THEN
   !destroy current grid
   filename = rsMin % file_name
   varname = rsMin % var_name
   CALL GridDestroy (rsMin)
   !read new grid in netcdf file
   CALL NewGrid (rsMin, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF


  
RETURN
END SUBROUTINE PlantsParameterUpdate





!==============================================================================
!| Description:
!  Update plants state variables
SUBROUTINE  PlantsGrow &
!
(time, radiation, temperature, swc, sfc, swp, rh, co2)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(IN) :: time !!current time 
TYPE (grid_real), INTENT(IN) :: radiation !!shortwave radiation [w/m2]
TYPE (grid_real), INTENT(IN) :: temperature !!air temperature [°C]
TYPE (grid_real), INTENT(IN) :: swc !! soil water content [m3/m3]
TYPE (grid_real), INTENT(IN) :: sfc !! soil field capacity [m3/m3]
TYPE (grid_real), INTENT(IN) :: swp !! soil wilting point [m3/m3]
TYPE (grid_real), INTENT(IN) :: rh !! air relative humidity  [0-1]
REAL (KIND = float), OPTIONAL, INTENT(IN) :: co2 !! CO2 [ppm]


!local declarations:
CHARACTER (LEN = 300) :: filename, varname
INTEGER (KIND = short) :: k, l
TYPE (PlantsCohort), POINTER :: cohort
REAL (KIND = float) :: fSWC !! GPP soil water content modifier 
REAL (KIND = float) :: fAGE !! GPP age modifier
REAL (KIND = float) :: fTEMP !! GPP temperature modifier
REAL (KIND = float) :: fVPD !! GPP vapor pressure deficit modififer
REAL (KIND = float) :: fCO2 !! GPP CO2 modifier
!REAL (KIND = float) :: deltaMroot !! root mass increment
!REAL (KIND = float) :: deltaMstem !! stem mass increment
REAL (KIND = float) :: rootaf !!root allocation factor
REAL (KIND = float) :: stemaf !! stem allocation factor
REAL (KIND = float) :: leafaf !! leaf allocation factor
REAL (KIND = float) :: pfs !!used to compute stem allocation factor
LOGICAL             :: is_new_year !! flag to mark new year has begun
TYPE (Practice)     :: practice


!------------end of declaration------------------------------------------------  

! need to update vegetation parameters from file?
IF ( simulatePlants == 0) THEN
    !check leaf area index
    IF (  time == lai % next_time ) THEN
        !destroy current grid
        filename = lai % file_name
        varname = lai % var_name
        CALL GridDestroy (lai)
        !read grid in netcdf file
        CALL NewGrid (lai, TRIM(filename), NET_CDF, &
                       variable = TRIM(varname), time = time)	           
    END IF
    
    !check vegetation fraction
    IF (  time == fvcover % next_time ) THEN
        !destroy current grid
        filename = fvcover % file_name
        varname = fvcover % var_name
        CALL GridDestroy (fvcover)
        !read grid in netcdf file
        CALL NewGrid (fvcover, TRIM(filename), NET_CDF, &
                       variable = TRIM(varname), time = time)           
    END IF
    
    !check plants height
    IF (  time == plantsHeight % next_time ) THEN
        !destroy current grid
        filename = plantsHeight % file_name
        varname = plantsHeight % var_name
        CALL GridDestroy (plantsHeight)
        !read grid in netcdf file
        CALL NewGrid (plantsHeight, TRIM(filename), NET_CDF, &
                       variable = TRIM(varname), time = time)           
    END IF
    
    !check minimum stomatal resistance
    IF (  rsMinLoaded .AND. time == rsMin % next_time ) THEN
        !destroy current grid
        filename = rsMin % file_name
        varname = rsMin % var_name
        CALL GridDestroy (rsMin)
        !read grid in netcdf file
        CALL NewGrid (rsMin, TRIM(filename), NET_CDF, &
                       variable = TRIM(varname), time = time)           
    END IF

ELSE !simulate forest dynamics
    year_new = GetYear (time)
    DO k = 1, count_stands 
        
        cohort => forest (k) % first
        DO l = 1, forest (k) % lenght !loop through all cohorts
            
            cohort % stem_yield = 0. 
            !check management practices
            IF ( plants_management ) THEN
               IF ( ALLOCATED (forest (k) % thinning % cuts) ) THEN
                  practice = forest (k) % thinning
                  
                  IF (time > practice % next .AND. practice % current < SIZE ( practice % cuts )) THEN 
    
                     practice % current = practice % current + 1
                     
                     IF ( practice % cuts ( practice % current ) % reforestation ) THEN
                         !set species parameters
                         cohort % species = species ( practice % cuts ( practice % current ) % species)
                     END IF
                     
                  
                     CALL ApplyPlantsManagement (time, practice, cohort % density, &
                                              cohort % mass_root, cohort % mass_stem, &
                                              cohort % mass_leaf, cohort % mass_total, &
                                              cohort % lai, cohort % canopy_cover, &
                                              cohort % age, cohort % height, cohort % dbh, &
                                              cohort % stem_yield )
                     
                     
                      !update crown diameter
                      cohort % crown_diameter = CrownDiameter ( dbh = cohort % dbh, &
                                                 den = cohort % density, &
                                                 denmin = cohort % species % denmin, &
                                                 denmax = cohort % species % denmax, &
                                                 dbhdcmin = cohort % species % dbhdcmin, &
                                                 dbhdcmax = cohort % species % dbhdcmax ) 
                     !update canopy cover
                     cohort % canopy_cover = CanopyCover (  cohort % crown_diameter, cohort % density ) 
                     
                     !update time of next cut
                     IF ( practice % current < SIZE ( practice % cuts ) ) THEN
                         practice % next = practice % cuts (  practice % current + 1 ) % time
                     END IF
                     
                     !copy practice back 
                     forest (k) % thinning = practice 
                     
                  END IF
               END IF
            
            END IF
            
            
            !update age
            IF ( DayOfYear (time) == 1 .AND. year_new /= year_prev ) THEN
                is_new_year = .TRUE.
                cohort % age = cohort % age + 1
                cohort % mass_stem_year_previous = cohort % mass_stem
            ELSE
                is_new_year = .FALSE.
            END IF
            
             !compute APAR
             cohort % apar = AparCalc (rad = radiation % mat (forest (k) % i, forest (k) % j) , &
                                       lai = cohort % lai, &
                                       k = cohort % species % k, &
                                       alb = cohort % species % albedo, &
                                       dt = dtPlants)
             
             
             !compute CO2 modifier
             IF ( useCO2modifier ) THEN
                 fCO2 = CO2mod (co2,  cohort % age )
             ELSE
                 fCO2 = 1.
             END IF
             
             
             
             !compute age modifier
             fAGE = AGEmod ( cohort % age, cohort % species % agemax)
             
             
             !compute temperature modifier 
             fTEMP = TEMPmod (Ta = temperature % mat (forest (k) % i, forest (k) % j), &
                              Tmin =  cohort % species % Tmin, &
                              Tmax =  cohort % species % Tmax, &
                              Topt =  cohort % species % Topt )
             
        
             !compute soil water content modifier
             fSWC = SWCmod (swc = swc % mat (forest (k) % i, forest (k) % j), &
                            wp = swp % mat (forest (k) % i, forest (k) % j), &
                            fc = sfc % mat (forest (k) % i, forest (k) % j), &
                            theta = cohort % species % theta_fswc )
             
             !compute vapor pressure deficit modifier
             fVPD = VPDmod (Ta = temperature % mat (forest (k) % i, forest (k) % j), &
                            RH = rh % mat (forest (k) % i, forest (k) % j), &
                            kd = cohort % species % theta_fvpd)
             
             
             !compute gross primary production
             cohort % gpp = cohort % apar * & ! absorbed photosynthetically active radiation molPAR m-2
                            cohort % species % alpha * &! canopy quantum efficiency (molC/molPAR)
                            C_molar_mass * & !conversion to Kg
                            hectare / & ! m2 in one hectare
                            1000. * & !conversion to t
                            fAGE * & !age modifier
                            fTEMP * & !temperature modifier
                            fSWC * & !soil water content modifier
                            fVPD * &!vapor pressure deficit modifier
                            fCO2
            
             
             !compute net primary production
             cohort % npp = cohort % gpp * cohort % species % GPPtoNPP
             
             !compute root allocation factor 
             rootaf = 0.5 / ( 1. + 2.5 * fAGE * fTEMP * fSWC )
             
             !compute stem allocation factor 
             pfs = (cohort % species % fprn * cohort % species % fpra) / &
                   (cohort % species % sprn * cohort % species % spra) * &
                   (cohort % density * cohort % dbh / 100.) ** &
                   (cohort % species % fprn - cohort % species % sprn)
             
             stemaf = ( 1. - rootaf ) / ( 1. + pfs )
             
             !compute leaf allocation factor
             leafaf = 1. - rootaf - stemaf
             
             ! update stem biomass
             cohort % mass_stem = cohort % mass_stem + stemaf * cohort % npp
             
             !update root biomass
             cohort % mass_root =  cohort % mass_root + & !current mass
                                   rootaf * cohort % npp - & !mass increment
                                   cohort % species % rtr * cohort % mass_root * dtPlants !turnover rate
             
             ! update leaf biomass and leaf area index
             CALL GrowLeaf (npp = cohort % npp , &
                            af = leafaf, &
                            Ta = temperature % mat (forest (k) % i, forest (k) % j), &
                            Tcold = cohort % species % tcold_leaf, &
                            swc = swc % mat (forest (k) % i, forest (k) % j), &
                            swp = swp % mat (forest (k) % i, forest (k) % j), &
                            sfc = sfc % mat (forest (k) % i, forest (k) % j), &
                            tr = cohort % species % ltr, &
                            sla = cohort % species % sla, &
                            mleaf = cohort % mass_leaf, &
                            lai = cohort % lai     )
             
             !update total biomass
             cohort % mass_total = cohort % mass_root + cohort % mass_stem + cohort % mass_leaf
             
             !update dbh and height tree every new year
             !IF (is_new_year) THEN
                 !CALL GrowDBH (cc = cohort % canopy_cover, &
                 !              hdmin = cohort % species % hdmin, & 
                 !              hdmax = cohort % species % hdmax, &
                 !              ws = cohort % mass_stem, &
                 !              dws = cohort % mass_stem - cohort % mass_stem_year_previous, &
                 !              !dws = stemaf * cohort % npp, &
                 !              DBH = cohort % dbh, &
                 !              height = cohort % height)
             
                  CALL GrowDBHech2o (cc = cohort % canopy_cover, &
                               hdmin = cohort % species % hdmin, & 
                               hdmax = cohort % species % hdmax, &
                               ws = cohort % mass_stem, &
                               !dws = cohort % mass_stem - cohort % mass_stem_year_previous, &
                               dws = stemaf * cohort % npp, &
                               DBH = cohort % dbh, &
                               height = cohort % height, &
                               tree_density = cohort % density, &
                               wood_density = cohort % species % wood_density, &
                               age = cohort % age, &
                               maxage = cohort % species % agemax )
             !END IF
                  
             !update crown diameter
             cohort % crown_diameter = CrownDiameter ( dbh = cohort % dbh, den = cohort % density, &
                                                 denmin = cohort % species % denmin, &
                                                 denmax = cohort % species % denmax, &
                                                 dbhdcmin = cohort % species % dbhdcmin, &
                                                 dbhdcmax = cohort % species % dbhdcmax ) 
             
            !update canopy cover
            cohort % canopy_cover = CanopyCover (  cohort % crown_diameter, cohort % density ) 
            
            
            !apply mortality
            IF (mortality) THEN
                CALL KillPlants (dtPlants, cohort % age, cohort % species % agemax, &
                                 cohort % species % ms, cohort % species % mf, &
                                 cohort % species % mr, cohort % species % wSx1000, &
                                 cohort % crown_diameter, cohort % density, &
                                 cohort % canopy_cover, cohort % mass_stem, &
                                 cohort % mass_root, cohort % mass_leaf, &
                                 cohort % mass_total)
            END IF
            
            !jump to next cohort
            cohort => cohort % next
        END DO
         
    END DO
    year_prev = year_new
    
    !update state variable grids
    CALL ForestToGrid (lai, 'lai') 
    CALL ForestToGrid (fvcover, 'fv')
    CALL ForestToGrid (gpp, 'gpp')
    CALL ForestToGrid (npp, 'npp')
    CALL ForestToGrid (carbonroot, 'root')
    CALL ForestToGrid (carbonstem, 'stem')
    CALL ForestToGrid (carbonleaf, 'leaf')
    CALL ForestToGrid (dbh, 'dbh')
    CALL ForestToGrid (plantsHeight, 'height')
    CALL ForestToGrid (density, 'density')
    CALL ForestToGrid (stemyield, 'stemyield')
END IF


RETURN
END SUBROUTINE PlantsGrow


!==============================================================================
!| Description:
! read species parameters
SUBROUTINE  ReadSpecies &
!
(inifile)

IMPLICIT NONE

CHARACTER (LEN = *), INTENT(in) :: inifile

!local declarations:
TYPE (IniList) :: speciesDB
INTEGER (KIND = short) :: nspecies
INTEGER (KIND = short) :: i
!------------------------------end of declarations----------------------------

! open and read configuration file
CALL IniOpen (inifile, speciesDB)

nspecies = IniReadInt ('species', speciesDB)

ALLOCATE (species (nspecies) )

DO i = 1, nspecies
    species (i) % name = IniReadString ('name', speciesDB, section = ToString(i))
    species (i) % phenology = IniReadString ('phenology', speciesDB, section = ToString(i))
    species (i) % alpha = IniReadReal ('alpha', speciesDB, section = ToString(i))
    species (i) % GPPtoNPP = IniReadReal ('GPP-NPP', speciesDB, section = ToString(i))
    species (i) % k = IniReadReal ('k', speciesDB, section = ToString(i))
    species (i) % cra = IniReadReal ('cra', speciesDB, section = ToString(i))
    species (i) % crb = IniReadReal ('crb', speciesDB, section = ToString(i))
    species (i) % crc = IniReadReal ('crc', speciesDB, section = ToString(i))
    species (i) % as = IniReadReal ('as', speciesDB, section = ToString(i))
    species (i) % ns = IniReadReal ('ns', speciesDB, section = ToString(i))
    !species (i) % saf = IniReadReal ('stem-allocation-factor', speciesDB, section = ToString(i))
    !species (i) % raf = IniReadReal ('root-allocation-factor', speciesDB, section = ToString(i))
    !species (i) % laf = IniReadReal ('leaf-allocation-factor', speciesDB, section = ToString(i))
    species (i) % dbhdcmin = IniReadReal ('dbhdcmin', speciesDB, section = ToString(i))
    species (i) % dbhdcmax = IniReadReal ('dbhdcmax', speciesDB, section = ToString(i))
    species (i) % denmin = IniReadReal ('denmin', speciesDB, section = ToString(i))
    species (i) % denmax = IniReadReal ('denmax', speciesDB, section = ToString(i))
    species (i) % agemax = IniReadReal ('agemax', speciesDB, section = ToString(i))
    species (i) % tmin = IniReadReal ('tmin', speciesDB, section = ToString(i))
    species (i) % tmax = IniReadReal ('tmax', speciesDB, section = ToString(i))
    species (i) % topt = IniReadReal ('topt', speciesDB, section = ToString(i))
    species (i) % theta_fswc = IniReadReal ('phi-theta', speciesDB, section = ToString(i))
    species (i) % theta_fvpd = IniReadReal ('phi-ea', speciesDB, section = ToString(i))
    species (i) % fpra = IniReadReal ('fpra', speciesDB, section = ToString(i))
    species (i) % fprn = IniReadReal ('fprn', speciesDB, section = ToString(i))
    species (i) % spra = IniReadReal ('spra', speciesDB, section = ToString(i))
    species (i) % sprn = IniReadReal ('sprn', speciesDB, section = ToString(i))
    species (i) % ltr = IniReadReal ('leaf-turnover', speciesDB, section = ToString(i)) 
    species (i) % rtr = IniReadReal ('root-turnover', speciesDB, section = ToString(i)) 
    species (i) % tcold_leaf = IniReadReal ('tcold-leaf', speciesDB, section = ToString(i))  
    species (i) % sla = IniReadReal ('sla', speciesDB, section = ToString(i))  
    species (i) % hdmin = IniReadReal ('hdmin', speciesDB, section = ToString(i)) 
    species (i) % hdmax = IniReadReal ('hdmax', speciesDB, section = ToString(i)) 
    species (i) % albedo = IniReadReal ('albedo', speciesDB, section = ToString(i)) 
    species (i) % laimax = IniReadReal ('laimax', speciesDB, section = ToString(i))
    species (i) % canopymax = IniReadReal ('canopymax', speciesDB, section = ToString(i))
    species (i) % wood_density = IniReadReal ('wood-density', speciesDB, section = ToString(i))
    species (i) % ms = IniReadReal ('ms', speciesDB, section = ToString(i))
    species (i) % mr = IniReadReal ('mr', speciesDB, section = ToString(i))
    species (i) % mf = IniReadReal ('mf', speciesDB, section = ToString(i))
    species (i) % wSx1000 = IniReadReal ('wSx1000', speciesDB, section = ToString(i))
END DO

!  species DB is complete. Deallocate ini database
CALL IniClose (speciesDB)  
    
RETURN
END SUBROUTINE ReadSpecies


!==============================================================================
!| Description:
! set initial condition
SUBROUTINE  SetIC &
!
(inifile)

IMPLICIT NONE

CHARACTER (LEN = *), INTENT(in) :: inifile

!local declarations:
TYPE (IniList) :: icDB
INTEGER (KIND = short) :: ncohorts
TYPE (grid_integer) :: cohort_map
INTEGER (KIND = short) :: i, j, k, l
INTEGER (KIND = short) :: ispecies
REAL (KIND = float) :: dbh 
TYPE (PlantsCohort), POINTER :: cohort

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

! open and read configuration file
CALL IniOpen (inifile, icDB)

!number of cohorts
ncohorts = IniReadInt ('cohorts', icDB)

!starting year
year_new = IniReadInt ('year', icDB)
year_prev = IniReadInt ('year', icDB)

!cohort map
CALL GridByIni (icDB, cohort_map, section = 'cohorts-map') 


DO k = 1, count_stands
    !aggiungere configurazione ricorsiva quando ci sono più coorti per cella
    !al momento solo una coorte per stand (cella)
    cohort => forest (k) % first
    DO l = 1, forest (k) % lenght
       ispecies = IniReadInt ('species', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % species = species (ispecies)
       cohort % age = IniReadInt ('age', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % density = IniReadReal ('density', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % height = IniReadReal ('height', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % dbh = IniReadReal ('dbh', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % lai = IniReadReal ('lai', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % mass_stem = IniReadReal ('stem-biomass', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % mass_root = IniReadReal ('root-biomass', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % mass_leaf = IniReadReal ('leaf-biomass', icDB, section = ToString   (cohort_map % mat (forest (k) % i, forest (k) % j )  ) )
       cohort % mass_total = cohort % mass_stem + cohort % mass_root + cohort % mass_leaf
       cohort % crown_diameter = CrownDiameter ( dbh = cohort % dbh, den = cohort % density, &
                                                 denmin = species (ispecies) % denmin, &
                                                 denmax = species (ispecies) % denmax, &
                                                 dbhdcmin = species (ispecies) % dbhdcmin, &
                                                 dbhdcmax = species (ispecies) % dbhdcmax ) 
       cohort % canopy_cover = CanopyCover (  cohort % crown_diameter, cohort % density )
       cohort => cohort % next
    END DO
END DO

! Deallocate memory
CALL IniClose (icDB)  
    
RETURN
END SUBROUTINE SetIC



!==============================================================================
!| Description:
! Computes the absorbed photosynthetically active radiation.
! The amount of light available for the plant that is going to define 
! its growth rate. The light that the plant could absorb is determined according 
! to the canopy cover. The amount of absorbed photosynthetically active radiation 
! is usually computed following  Lambert-Beer law . 
!
! References:
!
!      McCree, K.J.: Test of current definitions of photosynthetically active 
!      radiation against leaf photosynthesis data. Agr. Forest. Meteorol, 
!      10 , 443-453, 1972.
FUNCTION  AparCalc &
!
(rad, lai, k, alb, dt) &
!
RESULT (apar)

IMPLICIT NONE

!arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: rad !! shortwave direct radiation (w/m2)
REAL (KIND = float), INTENT(IN) :: lai !! leaf area index (m2/m2)
REAL (KIND = float), INTENT(IN) :: k !! extinction coefficient for absorption of PAR by canopy
REAL (KIND = float), INTENT(IN) :: alb !! plant albedo (0-1)
INTEGER (KIND = short), INTENT(IN) :: dt  !!time step (s)

!local declarations:
REAL (KIND = float) :: apar
REAL (KIND = float) :: par0  !!photosynthetically active radiation (molPAR m-2)
REAL (KIND = float), PARAMETER :: radToPAR = 4.57 !conversion factor 1 W m-2 = 4.57 µmol m-2 s-1 (McCree, 1972)
!------------------------------end of declarations----------------------------



!devo considerare l'albedo cioè è una short wave net radiation ? o quella  incidente ?
par0 = rad * ( 1. - alb ) * radToPAR * dt ! (µmol m-2)

apar = par0 * ( 1. - EXP ( -k * lai) ) / 1000000. ! (molPAR m-2)
    
RETURN
END FUNCTION AparCalc


!==============================================================================
!| Description:
! update leaf biomass and leaf area index
!
! References:
!
!Maneta, M. P., and N. L. Silverman, 2013: A spatially distributed model to 
! simulate water, energy, and vegetation dynamics using information from 
! regional climate models. Earth Interact., 17
! doi:10.1175/2012EI000472.1.
! 
! Arora, V. K., and G. J. Boer, 2005: A parameterization of leaf phenology 
! for the terrestrial ecosystem component of climate models. 
! Global Change Biol., 11, 39–59.
!
SUBROUTINE  GrowLeaf &
!
(npp, af, Ta, Tcold, swc, swp, sfc, tr, sla, mleaf, lai) 

IMPLICIT NONE

!arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: npp !!net primary production [t/ha]
REAL (KIND = float), INTENT(IN) :: af !!allocation factor [0-1]
REAL (KIND = float), INTENT(IN) :: Ta !! air temperature [°C]
REAL (KIND = float), INTENT(IN) :: Tcold !! temperature threshold that accelerates leaf turnover (°C)
REAL (KIND = float), INTENT(IN) :: swc !! soil water content [m3/m3]
REAL (KIND = float), INTENT(IN) :: swp !! soil wilting point [m3/m3]
REAL (KIND = float), INTENT(IN) :: sfc !! soil field capacity [m3/m3]
REAL (KIND = float), INTENT(IN) :: tr !! leaf turnover rate [s-1]
REAL (KIND = float), INTENT(IN) :: sla !! specific leaf area  [m2/Kg]

!arguments with intent(inout):
REAL (KIND = float), INTENT(INOUT) :: mleaf !!mass of leaf [t/ha]
REAL (KIND = float), INTENT(INOUT) :: lai !! leaf area index [m2/m2]

!local declarations:
REAL (KIND = float) :: deltaMleaf !! leaf mass increment
REAL (KIND = float) :: swc_stress !! soil water content stress that affects leaf increment [0-1]
REAL (KIND = float) :: temp_stress !! temperature stress that affects leaf increment [0-1]
REAL (KIND = float) :: actual_tr !! actual leaf tirnover rate with temperature and soil water content stresses

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

!compute temperature stress
IF (Ta >= Tcold) THEN
    temp_stress = 1.
ELSE IF ((Tcold - 5.) < Ta < Tcold) THEN
    temp_stress = (Ta - Tcold - 5.) / 5.
ELSE IF (Ta <= (Tcold - 5.)) THEN
    temp_stress = 0.
END IF
!temp_stress = MAX (0., MIN (1., ((Ta - (Tcold - 5.) / 5. ) ) ))

!compute hydrologic stress
swc_stress = MAX (0., MIN (1., (swc - swp) / (sfc - swp) ) )  

!compute actual leaf turnover rate: base turnover + 0.001 hydrologic stress + temperature stress
actual_tr = tr * ( 1. + 0.001 * (1. - swc_stress)**3. + (1. - temp_stress)**3. )

!leaf mass increment
deltaMleaf = af * npp 

!update leaf biomass
mleaf = mleaf + deltaMleaf
             
!update leaf area index
lai = lai + (deltaMleaf * sla/ 10.) -  (actual_tr * dtPlants) *LAI !10 = conversion factor 1000 [kg/t] / 10000 [m2/hectare]

!check lower boundary
IF (mleaf < 0.) THEN
    mleaf = 0.
END IF

IF (lai < 0.) THEN
    lai = 0.
END IF

RETURN
END SUBROUTINE GrowLeaf


!==============================================================================
!| Description:
! update DBH and height tree according to ech2o model by Maneta
!
SUBROUTINE  GrowDBHech2o &
!
(cc, hdmin, hdmax, ws, dws, DBH, height, tree_density, wood_density, age, maxage) 
 
    
IMPLICIT NONE

!arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: cc !! canopy cover
REAL (KIND = float), INTENT(IN) :: hdmin !! H/D ratio in carbon partitioning for low density
REAL (KIND = float), INTENT(IN) :: hdmax !! H/D ratio in carbon partitioning for high density
REAL (KIND = float), INTENT(IN) :: ws !! stem biomass (t/ha)
REAL (KIND = float), INTENT(IN) :: dws !! stem biomass increment (t/ha)
REAL (KIND = float), INTENT(IN) :: tree_density !! tree density (trees/ha)
REAL (KIND = float), INTENT(IN) :: wood_density !! wood density (kg/m3)
REAL (KIND = float), INTENT(IN) :: age !! tree age (year)
REAL (KIND = float), INTENT(IN) :: maxage !!tree max age (year)

!arguments with intent inout
REAL (KIND = float), INTENT(INOUT) :: DBH  !diameter at brest height (cm)
REAL (KIND = float), INTENT(INOUT) :: height  !tree height (m)

!local declarations
REAL (KIND = float) :: delta_tree_mass !!tree mass increment (kg/tree)
REAL (KIND = float) :: fhd !! grow factor that depends on the height-to-diameter ratio
REAL (KIND = float) :: dDBH !DBH increment (cm)
REAL (KIND = float) :: dheight !height increment (m)

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

!compute mass increment per tree
delta_tree_mass = dws / tree_density * 1000. ! mass increment per tree (kg)

!set grow factor
IF (height / DBH >= hdmin .AND. cc < 0.95) THEN
    fhd = hdmin
ELSE IF (height / DBH <= hdmax .AND. cc >= 0.95 .AND. age < maxage / 2) THEN
    fhd = hdmin
ELSE IF (height / DBH <= hdmax .AND. cc >= 0.95 .AND. age > maxage / 2 ) THEN
    fhd = hdmax
ELSE IF (height / DBH < hdmin) THEN
    fhd = hdmax
ELSE IF (height / DBH > hdmax) THEN
    fhd = 0.5 * hdmin
ELSE IF ( age > 0.7 * maxage) THEN
    fhd = 0.
END IF

!compute DBH increment
dDBH = 4. * delta_tree_mass / ( wood_density * pi * (DBH / 100.)**2. * ( 2. * height / (DBH / 100.) + fhd * 100. ) ) * 100.


!compute tree height increment
dheight = dDBH * fhd

!update dbh and height tree
 DBH = DBH + dDBH
 height = height + dheight

!write(*,*) dDBH, dheight
!pause


RETURN
END SUBROUTINE GrowDBHech2o


!==============================================================================
!| Description:
! update DBH and height tree every new year
!
SUBROUTINE  GrowDBH &
!
(cc, hdmin, hdmax, ws, dws, DBH, height) 
    
IMPLICIT NONE

!arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: cc !! canopy cover
REAL (KIND = float), INTENT(IN) :: hdmin !! H/D ratio in carbon partitioning for low density
REAL (KIND = float), INTENT(IN) :: hdmax !! H/D ratio in carbon partitioning for high density
REAL (KIND = float), INTENT(IN) :: ws !! stem biomass (t/ha)
REAL (KIND = float), INTENT(IN) :: dws !! stem biomass increment (t/ha)

!arguments with intent inout
REAL (KIND = float), INTENT(INOUT) :: DBH  !diameter at brest height (cm)
REAL (KIND = float), INTENT(INOUT) :: height  !tree height (m)

!local declarations
REAL (KIND = float) :: hdeff
REAL (KIND = float) :: dDBH !DBH increment (cm)
REAL (KIND = float) :: dheight !height increment (m)

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

IF ( cc <= 0.95) THEN !low density
    hdeff = hdmin + (hdmax - hdmin) * cc
    
ELSE ! high density
    hdeff = hdmax  
END IF

!compute height increment
!dheight = hdeff / (  hdeff / height + 200. / dbh ) * dws / ws

dheight = hdeff / (  hdeff / height) + (200. / dbh ) * dws / ws

!compute DBH increment
dDBH = ( height / hdeff ) + ( 200. / dbh ) * ( dws / ws )


!update dbh and height tree
 DBH = DBH + dDBH
 height = height + dheight
 
 !write(*,*) dheight, hdeff, dws, ws
 !pause

RETURN
END SUBROUTINE GrowDBH


!==============================================================================
!| Description:
! fill in a grid with state variable values in forest stands
!
SUBROUTINE  ForestToGrid &
!
(grid, statevar) 
    
IMPLICIT NONE

!arguments with intent(in):
CHARACTER (LEN = *), INTENT(IN) :: statevar

!arguments with intent inout
TYPE (grid_real), INTENT (INOUT) :: grid

!local declarations
INTEGER (KIND = short) :: i, j, k


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

DO k = 1, count_stands
    i = forest (k) % i
    j = forest (k) % j 
    SELECT CASE (statevar)
        
    CASE ('lai') !leaf area index
        grid % mat (i,j) = forest (k) % first % lai
        
    CASE ('fv') !canopy cover
        grid % mat (i,j) = forest (k) % first % canopy_cover
        
    CASE ('gpp') !GPP
        grid % mat (i,j) = forest (k) % first % gpp * &
                           CellArea (grid,i,j) / hectare
        
    CASE ('npp') !NPP 
        grid % mat (i,j) = forest (k) % first % npp * &
                           CellArea (grid,i,j) / hectare
    CASE ('root') !root mass
        grid % mat (i,j) = forest (k) % first % mass_root * &
                           CellArea (grid,i,j) / hectare
        
    CASE ('stem') !stem mass
        grid % mat (i,j) = forest (k) % first % mass_stem * &
                           CellArea (grid,i,j) / hectare
        
    CASE ('leaf') !foliage mass
        grid % mat (i,j) = forest (k) % first % mass_leaf * &
                           CellArea (grid,i,j) / hectare
        
     CASE ('stemyield') !stem mass
        grid % mat (i,j) = forest (k) % first % stem_yield * &
                           CellArea (grid,i,j) / hectare
             
    CASE ('dbh') !diameter at brest heigth (cm)
        grid % mat (i,j) = forest (k) % first % dbh 
        
    CASE ('height') !tree height (m)
        grid % mat (i,j) = forest (k) % first % height
    
    CASE ('density') !tree density (trees/hectare)
        grid % mat (i,j) = forest (k) % first % density
        
   
        
    END SELECT
    
   
END DO


RETURN
END SUBROUTINE ForestToGrid

    
END MODULE Plants