!!  Compute canopy interception
!|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 - 25th Mar 2019  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 25/Mar/2019 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Compute canopy interception. The method implemented in this 
! module is according to SWAT model (canopy storage)
! 
 MODULE PlantsInterception
    ! Modules used: 
! 
USE DataTypeSizes, ONLY : &
    ! Imported Type Definitions:
    short, float, double  

USE LogLib,ONLY : &
   ! imported routines:
   Catch  

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

USE IniLib, ONLY: &
    !Imported routines:
    IniOpen, SectionIsPresent, &
    IniClose, IniReadReal, KeyIsPresent, &
    IniReadInt, &
   !Imported type definitions:
   IniList

USE GridOperations, ONLY : &
   !imported routines:
   GridByIni

USE Units, ONLY: &
    !imported variables:
    millimeter

USE StringManipulation, ONLY: &
  !Imported routines:
  ToString

IMPLICIT NONE
! Global (i.e. public) Declarations:

!public variables
INTEGER (KIND = short) :: dtCanopyInterception
TYPE (grid_real) :: canopyStorage !!water canopy storage (mm)
TYPE (grid_real) :: canopyPT  !! transpiration from canopy (m/s)
TYPE (grid_real) :: laimax !! maximum leaf area index value (m2/m2)
TYPE (grid_real) :: canopymax !! maximum canopy storage capacity (m)
INTEGER (KIND = short) :: interceptionParametersByMap = 1 !!set if parameters are load from map (1) or set according to plant species properties (0)

!public routines
PUBLIC :: Throughfall
PUBLIC :: AdjustPT


! local (i.e. private) declarations:

!private variables


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

!==============================================================================
!| initialize variables for computing rainfall interception
SUBROUTINE InterceptionInit &
!
( filename, domain )  

IMPLICIT NONE

!Arguments with intent in
CHARACTER (LEN = *), INTENT(in) :: filename !!configuration file name
TYPE (grid_integer), INTENT(IN) :: domain  !!analysis domain

!local declarations:
TYPE(IniList) :: iniDB !! configuration info
REAL (KIND = float) :: ics !!initial canopy storage value
    
!-------------------------end of declarations----------------------------------

CALL IniOpen (filename, iniDB)

!read if parameters are loaded from map files
IF (KeyIsPresent ( 'parameters-by-map ', iniDB) ) THEN
  interceptionParametersByMap =  IniReadInt ('parameters-by-map', iniDB )
ELSE
   CALL Catch ('error', 'PlantsInterception', &
       'missing parameters-by-map in PlantsInterception configuration file')
END IF


IF ( interceptionParametersByMap == 1) THEN !load maps

        ! load laimax map
        IF (SectionIsPresent ( 'laimax', iniDB) ) THEN
           CALL GridByIni (iniDB, laimax, section = 'laimax')
        ELSE
           CALL Catch ('error', 'PlantsInterception', &
               'missing laimax map in PlantsInterception configuration file')
        END IF


        ! load canopymax map
        IF (SectionIsPresent ( 'canopymax', iniDB) ) THEN
           CALL GridByIni (iniDB, canopymax, section = 'canopymax')
        ELSE
           CALL Catch ('error', 'PlantsInterception', &
               'missing canopymax map in PlantsInterception configuration file')
        END IF
        
END IF

!set canopy storage initial value
IF ( KeyIsPresent ('cold', iniDB, section = 'canopy-storage') ) THEN
    
    !initial value
    ics = IniReadReal ('cold', iniDB, section = 'canopy-storage')
    CALL Catch ('info', 'PlantsInterception: ',     &
                 'initial canopy storage value (m): ', &
                 argument = ToString(ics))
    		
    CALL NewGrid (canopyStorage, domain, ics / millimeter)
ELSE !read map (hot start)
    CALL GridByIni (iniDB, canopyStorage, section = 'canopy-storage')
END IF

CALL IniClose (iniDB)


!allocate canopyET grid
CALL NewGrid (canopyPT, domain, 0.)


RETURN
END SUBROUTINE InterceptionInit
    
 
!==============================================================================
!| calculate the effective rainfall through canopy
SUBROUTINE Throughfall &
!
(rain , lai, domain, fv, raineff)
   
 IMPLICIT NONE 

!Arguments with intent in
TYPE (grid_real), INTENT(IN) :: rain !!rainfall rate (m/s)
TYPE (grid_real), INTENT(IN) :: lai   !!leaf area index (m2/m2)
TYPE (grid_integer), INTENT(IN) :: domain !!analysis domain
TYPE (grid_real), INTENT(IN) :: fv !!fraction of vegetation covering the cell (0-1)

!arguments with intent inout
TYPE (grid_real), INTENT(INOUT) :: raineff  !!effecttive rainfall rate (throughfall) (m/s)

!local declarations:

REAL (KIND = float) :: canopyinter   !! is the amount of water intercepted at current time step
REAL (KIND = float) :: canopymaxi  !! maximum canopy storage at current day's leaf area index
REAL (KIND = float) :: raindt !!rainfall amount in the given dt (mm)
INTEGER (KIND = short) :: i, j
!-------------------------end of declarations----------------------------------


DO j = 1, domain % jdim
    DO i = 1, domain % idim
        IF (domain % mat (i,j) /= domain % nodata ) THEN
    
            !compute maximum amount of water that can be held in canopy  storage as a function of lai
            canopymaxi = canopymax % mat (i,j) * lai % mat (i,j) / laimax % mat (i,j) / millimeter ! unit = mm

            !compute amount of rainfall of the given dt in mm, before canopy interception is removed
            raindt = rain % mat (i,j) * dtCanopyInterception / millimeter

            !update amount of intercepted rainfall. The canopy storage is filled before any 
            ! water is allowed to reach the ground
            IF ( raindt <= canopymaxi - canopyStorage % mat (i,j) ) THEN !rain is lower than the available storage
                canopyStorage % mat (i,j) = canopyStorage % mat (i,j) + raindt  !all rain is intercepted
                raineff % mat (i,j) = 0.
            ELSE !rain amount is greater than available storage
                raineff % mat (i,j) = raindt - ( canopymaxi - canopyStorage % mat (i,j) )
                canopyStorage % mat (i,j) = canopymaxi !canopy storage reaches maximum value
    
            END IF
 
            !convert raineff from mm to m/s
            raineff % mat (i,j) = raineff % mat (i,j) * millimeter / dtCanopyInterception
            
            !canopy affects only fraction of ground surface covered by vegetation,
            !adjust raineff for bare soil
            raineff % mat (i,j) = fv % mat (i,j) * raineff % mat (i,j) + &
                                  ( 1. - fv % mat (i,j) ) * rain % mat (i,j)
            
           
        END IF
   END DO
END DO

RETURN

END SUBROUTINE Throughfall 



!==============================================================================
!| Subroutine to adjust the actual evaporation to the intercepted rainfall.
! The amount of water intercepted by the canopy is going to contribute to 
! the evaporation rate. When calculating the evaporation within a forest,
! the model tends to remove as much as possible from the water intercepted 
! by the canopy. This step will be carried out before the adjustement of 
! the evapotranspiration to soil moisture. 
SUBROUTINE AdjustPT  &
!
(fv, domain, pt, dtpet)

   
IMPLICIT NONE 

!Arguments with intent(in)
TYPE (grid_real), INTENT(IN) :: fv !!fraction of vegetation covering the cell
TYPE (grid_integer), INTENT(IN) :: domain !!analysis domain
INTEGER (KIND = short), INTENT(IN) :: dtpet !! dt of evapotranspiration computation

!Arguments with intent inout
TYPE (grid_real),  INTENT(INOUT) ::  pt !!potential transpiration  from soil [m/s] 


!local declarations:
REAL (KIND = float) :: ptdtpet !!potential transpiration amount in the given dt of PET (mm)
REAL (KIND = float) :: ptdtcanopy !!potential transpiration amount in the given dt of canopy (mm)
INTEGER (KIND = short) :: i, j
!REAL (KIND = float) :: Evapintercept !! amount of water evaporated from canopy storage
!REAL (KIND = float)  ::  petadj !potential evapotranspiration  [m/s]
!-------------------------end of declarations----------------------------------




!      petadj = pet - CANOPIN
!      if (petadj < 0.) then
!        CANOPIN = -petadj
!        Evapintercept = pet                                                                                                                                                  
!        petadj = 0.
!      
!      else
!       Evapintercept = CANOPIN
!        CANOPIN = 0.
!      endif
!
!pet= petadj
    
!!Taking into consideration the amount of the evapotranspiration, that the amount of water is going to evaporate first from the intercepted water,
!!the remaining evaporative demand is going to evaporate from the soil 
!!(like the normal calculation that is usually implemented for the evapotranspiration).

DO j = 1, domain % jdim
    DO i = 1, domain % idim
        IF (domain % mat (i,j) /= domain % nodata ) THEN
            
            IF (canopyStorage % mat (i,j) <= 0.) THEN !no water stored on canopy:  canopyPT = 0
                canopyPT % mat (i,j) =  0.
                CYCLE !go to next cell
            END IF
             
            !compute amount of potential transpiration of the given dt in mm
            ptdtcanopy = pt % mat (i,j) * dtCanopyInterception * 1000.
            ptdtpet    = pt % mat (i,j) * dtpet * 1000.
            
            IF ( ptdtcanopy <= canopyStorage % mat (i,j) ) THEN !canopy storage satisfies the transpiration demand
                canopyStorage % mat (i,j) = canopyStorage % mat (i,j) - ptdtcanopy !remove transpiration from canopy storage
                canopyPT % mat (i,j) =  pt % mat (i,j) !evaporation from canopy is at potential rate
                pt % mat (i,j) =  ( ptdtpet - ptdtcanopy ) / dtpet * millimeter !remove canopy evaporation  from total evaporaion
                
            ELSE !canopy storage satisfies partially the transpiration demand
                canopyPT % mat (i,j) = canopyStorage % mat (i,j) / dtCanopyInterception * millimeter
                
                pt % mat (i,j) =  ( ptdtpet - canopyStorage % mat (i,j) ) / dtpet * millimeter !remove canopy evaporation  from total evaporaion
                
                
                !pt % mat (i,j) = pt % mat (i,j) - canopyStorage % mat (i,j) / dtCanopyInterception * millimeter * &
                 !                                dtCanopyInterception/dtpet!reduce potential trasnpiration from soil
                canopyStorage % mat (i,j) = 0.
            END IF
           
        END IF
   END DO
END DO        
    
RETURN


END SUBROUTINE  AdjustPT





 
END MODULE PlantsInterception