!! Compute evapotranspiration
!|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.2 - 8th March 2023  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 08/Jun/2021 | Original code |
! | 1.1      | 27/Jan/2023 | subroutine HargreavesSamaniModified added |
! | 1.2      | 08/Mar/2023 | etEnergyBalanceGround |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
!   Module to compute  evapotranspiration
!
MODULE Evapotranspiration


! Modules used: 

USE DataTypeSizes, ONLY : &  
  ! Imported Parameters:
  float, short


USE Loglib, ONLY : &
  !Imported routines:
  Catch

USE StringManipulation, ONLY : &
    !imported routines
    ToString

USE GridLib, ONLY : &
  !Imported types:
  grid_real, grid_integer, &
  !Imported routines:
  NewGrid , GetDtGrid !, &
  !GridDestroy, ExportGrid, &
  !SetLongName, SetStandardName, &
  !SetUnits, SetCurrentTime, &
  !Imported parameters:
  !NET_CDF, ESRI_ASCII, &
  !ESRI_BINARY
    
USE GridOperations, ONLY : &
    !imported routines:
    CRSisEqual, GridByIni, &
    !imported operators:
    ASSIGNMENT (=)

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

USE Chronos, ONLY : &
    !imported definitions:
    DateTime, &
    !Imported routines:
    DayOfYear, &
    !Imported operations:
    OPERATOR (+), &
    OPERATOR (==), &
    OPERATOR (<), &
    ASSIGNMENT (=)

USE Units, ONLY : &
    !imported parameters:
    pi, millimeter, day

USE DomainProperties, ONLY : &
    !imported variables:
    mask, latCentroid

USE MorphologicalProperties, ONLY : &
    !imported variables:
    dem, dem_loaded

USE AirTemperature, ONLY : &
    !Imported variables:
    dtTemperature, temperatureAir

USE AirTemperatureDailyMax , ONLY : &
    !imported variables:
    dtTemperatureDailyMax, temperatureAirDailyMax

USE AirTemperatureDailyMin , ONLY : &
    !imported variables:
    dtTemperatureDailyMin, temperatureAirDailyMin

USE AirTemperatureDailyMean , ONLY : &
    !imported variables:
    dtTemperatureDailyMean, temperatureAirDailyMean

USE AirRelativeHumidity, ONLY : &
    !Imported variables:
    dtRelHumidity, relHumidityAir, hygrometers

USE WindFlux, ONLY : &
    !Imported variables:
    dtWindSpeed, windSpeed, &
    anemometersWS

USE SolarRadiation, ONLY : &
    !Imported routines:
    dtRadiation, netRadiation, netRadiationFAO

USE Plants, ONLY : &
    !Imported variables:
    fvcover, fvcoverLoaded, &
    lai, laiLoaded, &
    plantsHeight, plantsHeightLoaded, &
    rsMinLoaded, rsMin

USE ObservationalNetworks, ONLY : &
  !Imported types:
  ObservationalNetwork, &
  !Imported routines:
  ReadMetadata, ReadData !, &
  !CopyNetwork , &
  !Imported operators:
  !ASSIGNMENT (=)
    
USE Utilities, ONLY : &
    !imported routines:
    GetUnit

USE MeteoUtilities, ONLY: &
  !Imported routines:
  ReadField
    
IMPLICIT NONE 

! Global (i.e. public) Declarations: 
TYPE (grid_real) :: pet !!potential evapotranspiration rate [m/s]
TYPE (grid_real) :: pe !!potential evaporation rate [m/s]
TYPE (grid_real) :: pt !!potential transpiration rate [m/s]
INTEGER (KIND = short) :: dtET !! time step computation [s]



!Private declarations

INTEGER (KIND = short), PRIVATE, PARAMETER :: ET_MODELS             = 6
INTEGER (KIND = short), PRIVATE, PARAMETER :: PENMAN_MONTEITH       = 1
INTEGER (KIND = short), PRIVATE, PARAMETER :: PRIESTLEY_TAYLOR      = 2
INTEGER (KIND = short), PRIVATE, PARAMETER :: HARGREAVES            = 3
INTEGER (KIND = short), PRIVATE, PARAMETER :: HARGREAVES_MOD        = 4
INTEGER (KIND = short), PRIVATE, PARAMETER :: FAO56_PENMAN_MONTEITH = 5
INTEGER (KIND = short), PRIVATE, PARAMETER :: ENERGY_BALANCE        = 6

TYPE (grid_integer),    PRIVATE :: model_map
TYPE (grid_integer),    PRIVATE :: kc_code_map  !!map of crop coefficient code 

TYPE (grid_real),       PRIVATE :: kc_map !! crop coefficient map

INTEGER (KIND = short), PRIVATE :: model_vector ( ET_MODELS ) !!defines active models 
INTEGER (KIND = short), PRIVATE :: model  !!model to compute ET
INTEGER (KIND = short), PRIVATE :: model_assignment  !!method to assign computation
                                                     !!1 = one method for the entire domain, 2 = a map with model codes
INTEGER (KIND = short), PRIVATE :: fileunitKc 
INTEGER (KIND = short), PRIVATE :: dtKc = 0
INTEGER (KIND = short), PRIVATE :: interpolationMethodKc
INTEGER (KIND = short), PRIVATE :: dtGridKc

CHARACTER (LEN = 300), PRIVATE :: filenameKc

TYPE (DateTime),        PRIVATE :: tnewET
TYPE (DateTime),        PRIVATE :: tNewKc

LOGICAL, PRIVATE :: compute_hargreaves = .FALSE.
LOGICAL, PRIVATE :: useCropCoefficient = .FALSE.

TYPE (ObservationalNetwork), PRIVATE :: kc_stations  !! pseudo station to read kc time series values

!Public routines
PUBLIC :: ETinit
PUBLIC :: PETupdate
PUBLIC :: etEnergyBalanceGround


!Local routines
PRIVATE :: HargreavesSamani
PRIVATE :: HargreavesSamaniModified
PRIVATE :: PriestleyTaylor
PRIVATE :: PenmanMonteith
PRIVATE :: FAO56PenmanMonteith
PRIVATE :: KcUpdate

!=======
CONTAINS
!=======
! Define procedures contained in this module.

!==============================================================================
!| Description:
!   Initialize  evapotranspiration computation
SUBROUTINE ETinit &
!
( inifile, time )

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *), INTENT(IN) :: inifile
TYPE (DateTime), INTENT(IN) :: time

!local declarations:
TYPE (IniList) :: etini
INTEGER (KIND = short) :: i, j

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

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

!dt
IF (KeyIsPresent('dt', etini) ) THEN	
   dtET = IniReadInt ('dt', etini)
ELSE
    CALL Catch ('error', 'Evapotranspiration',   &
                'dt missing in configuration file' )
END IF

!initialize PET map
CALL NewGrid (pet, mask)
CALL NewGrid (pe, mask)
CALL NewGrid (pt, mask)


!model assignment method
IF (KeyIsPresent('model-assignment', etini ) ) THEN
    model_assignment = IniReadInt ('model-assignment', etini)
ELSE
   CALL Catch ('error', 'Evapotranspiration',   &
            'model-assignment missing in configuration file')
END IF   

!set model
IF (model_assignment == 1) THEN

    model = IniReadInt ('model', etini )
    
    CALL NewGrid (model_map, mask)
    
    CALL CheckSpecificProperties (model, etini)
    
    model_map = model
  
ELSE 
    model = -1
    
    !check map is present
    IF ( .NOT. SectionIsPresent ('model-map', etini) ) THEN
         CALL Catch ('error', 'Evapotranspiration',   &
            'model-map missing in configuration file')
    END IF
    
    !read map
    CALL GridByIni (etini, model_map, 'model-map')
    
    IF ( .NOT. CRSisEqual (mask, model_map, .TRUE.) ) THEN
        CALL Catch ('error', 'Evapotranspiration',   &
            'something seems wrong in model-map')
    END IF
    
    !scan for et models
    model_vector = -1
    DO i = 1, model_map % idim
        DO j = 1, model_map % jdim
            IF (model_map % mat(i,j) /= model_map % nodata) THEN
                model_vector (model_map % mat (i,j)) = model_map % mat (i,j)
            END IF
        END DO
    END DO

    DO i = 1, ET_MODELS
         CALL CheckSpecificProperties (model_vector (i), etini)
    END DO

END IF

!crop coefficient
IF ( KeyIsPresent ('use-crop-coefficient', etini) ) THEN
    IF (IniReadInt ('use-crop-coefficient', etini) == 1 ) THEN
        useCropCoefficient = .TRUE.
    END IF
END IF

IF (useCropCoefficient) THEN
    IF ( .NOT. SectionIsPresent ('crop-coefficient', etini) ) THEN
        CALL Catch ('error', 'Evapotranspiration',   &
                'crop-coefficient section is missing in configuration file' )
    END IF
    
    !allocate map
    CALL NewGrid (kc_map, mask)
    
    !file
    filenameKc = IniReadString ('file', etini, section = 'crop-coefficient')
    
    !interpolation method
    interpolationMethodKc = IniReadInt ('interpolation', etini, section = 'crop-coefficient')
    
    IF (interpolationMethodKc == 0) THEN !data are stored in net-cdf file
       IF ( KeyIsPresent ('variable', etini, section = 'crop-coefficient') ) THEN
          kc_map % var_name = IniReadString ('variable', etini, section = 'crop-coefficient')
       ELSE IF  (KeyIsPresent ('standard_name', etini, section = 'crop-coefficient') ) THEN
          kc_map % standard_name = IniReadString ('standard_name', etini, section = 'crop-coefficient')
       ELSE
          CALL Catch ('error', 'Evapotranspiration',   &
		       'standard_name or variable missing in section crop-coefficient' )
       END IF
       
       !Get the dt of imported  field. Assume dt is regular	
       dtKc = GetDtGrid (filenameKc, .TRUE.)
       
    ELSE !open file containing local measurements
       fileunitKc = GetUnit ()
	   OPEN ( fileunitKc, file = filenameKc (1:LEN_TRIM(filenameKc)), status = 'old')
       
       !populate  metadata
       CALL ReadMetadata (fileunitKc, kc_stations)
       
       !set dt
       dtKc = kc_stations % timeIncrement
        
       !read map of crop coefficient code
       CALL GridByIni (ini = etini, grid = kc_code_map, &
                section = 'crop-coefficient', subsection = 'code-map')  
    END IF
    
    !set time
    tnewKc = time
    
END IF

!set time
tnewET = time


!  Configuration terminated. Deallocate ini database
CALL IniClose (etini)  

RETURN
END SUBROUTINE ETinit


!==============================================================================
!| Description:
!   Compute  potential evapotranspiration 
SUBROUTINE PETupdate &
!
( time )

IMPLICIT NONE

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

!local declarations:
INTEGER (KIND = short)   :: i, j
INTEGER (KIND = short)   :: jDay !!julian day
REAL (KIND = float)      :: sunsetHourAngle !!sunset hour angle [radians]
REAL (KIND = float)      :: relDist !!relative distance between the earth and the sun
REAL (KIND = float)      :: declination !! declination
REAL (KIND = float)      :: etRad !!extra terrestrial solar radiation [mm/day]


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

!Update Kc
IF (dtKc > 0 ) THEN
	IF( time == tNewKc ) THEN
		CALL  KcUpdate (time)
		tNewKc = tNewKc + dtKc
	END IF
END IF	


IF ( time == tnewET )  THEN !it's time to update ET
    
    !set next time to update PET
    tnewET = tnewET + dtET
    
    IF ( compute_hargreaves ) THEN !compute values for applying Hargreaves
        
        !calculate solar declination
        jDay = DayOfYear(time,'noleap') 
        declination = 0.4093 * SIN( 2. * pi * jDay / 365. - 1.405 )
        
        !calculate sunset hour angle
        sunsetHourAngle = ACOS( - TAN (latCentroid) * TAN (declination) )
        
        !calculate relative distance between the earth and the sun
        relDist = 1. + 0.033 * COS(2. * pi * jDay / 365.)
        
        !calculate extra terrestrial solar radiation (Allen et al. 1998)
        etRad = 15.392 * relDist * ( sunsetHourAngle * SIN(latCentroid) * &
                                     SIN(declination) + COS(latCentroid) * &
                                    COS(declination) * SIN(sunsetHourAngle) )
        
    END IF
    
   
    
    DO i = 1, model_map % idim
        DO j = 1, model_map % jdim
            IF ( model_map % mat (i,j) /= model_map % nodata) THEN
                SELECT CASE (model_map % mat (i,j) )
                    
                CASE (HARGREAVES)
                   CALL HargreavesSamani (time, temperatureAirDailyMean % mat(i,j), &
                                         temperatureAirDailyMax % mat(i,j), &
                                         temperatureAirDailyMin % mat (i,j), &
                                         etRad, pet % mat(i,j) )
                   

                   pe % mat (i,j) = pet % mat(i,j)
                   
                   pt % mat (i,j) = pet % mat(i,j)
                   
                CASE (HARGREAVES_MOD)
                   CALL HargreavesSamaniModified (time, &
                                         temperatureAirDailyMean % mat(i,j), &
                                         temperatureAirDailyMax % mat(i,j), &
                                         temperatureAirDailyMin % mat (i,j), &
                                         etRad, dem % mat (i,j), &
                                         pet % mat(i,j) )
                   

                   pe % mat (i,j) = pet % mat(i,j)
                   
                   pt % mat (i,j) = pet % mat(i,j)
                   
                CASE (PRIESTLEY_TAYLOR)  
                    
                   CALL PriestleyTaylor (temperatureAir % mat (i,j), &
                                         netRadiation % mat (i,j), &
                                         fvcover % mat (i,j), &
                                         dem % mat (i,j), &
                                         pt % mat (i,j), &
                                         pe % mat (i,j), &
                                         pet % mat (i,j) )
                   
                CASE (PENMAN_MONTEITH) 
                    
                    
                    
               !CALL  PetPM (R = netRadiation % mat (i,j), &
               !            T = temperatureAir % mat (i,j), &
               !            W =  windSpeed % mat (i,j), &
               !            rs = rsMin % mat (i,j), &
               !            v = fvcover % mat (i,j), &
               !            lai = lai % mat (i,j), &
               !            z = plantsHeight % mat (i,j), &
               !            um = relHumidityAir % mat (i,j), &
               !            quota = dem % mat (i,j), &
               !            pet = pet % mat (i,j), &
               !            pe = pe % mat (i,j), &
               !            pt = pt % mat (i,j))
                    
                    
                    
                   CALL PenmanMonteith (temperatureAir % mat (i,j), &
                                         netRadiation % mat (i,j), &
                                         relHumidityAir % mat (i,j), &
                                         windSpeed % mat (i,j), &
                                         fvcover % mat (i,j), &
                                         plantsHeight % mat (i,j), &
                                         dem % mat (i,j), &
                                         anemometersWS % offsetZ, &
                                         hygrometers % offsetZ, &
                                         lai % mat (i,j), &
                                         rsMin % mat (i,j), &
                                         pt % mat (i,j), &
                                         pe % mat (i,j), &
                                         pet % mat (i,j) )
                   
                 CASE (FAO56_PENMAN_MONTEITH)  
                    
                   CALL FAO56PenmanMonteith (temperatureAir % mat (i,j), &
                                         netRadiation % mat (i,j), &
                                         netRadiationFAO % mat (i,j), &
                                         relHumidityAir % mat (i,j), &
                                         windSpeed % mat (i,j), &
                                         fvcover % mat (i,j), &
                                         dem % mat (i,j), &
                                         anemometersWS % offsetZ, &
                                         hygrometers % offsetZ, &
                                         lai % mat (i,j), &
                                         pt % mat (i,j), &
                                         pe % mat (i,j), &
                                         pet % mat (i,j) )
  
                 END SELECT
                 
                 
                 IF ( useCropCoefficient) THEN
                       pt % mat (i,j) = pet % mat(i,j) * kc_map % mat (i,j)
                       !update pet
                       pet % mat (i,j) = fvcover % mat (i,j) * pt % mat (i,j) + &
                                        ( 1. - fvcover % mat (i,j) ) *  pe % mat (i,j) 
                 ELSE
                       pt % mat (i,j) = pet % mat(i,j)
                 END IF
                 
                 !check negative values
                 IF ( pt % mat (i,j) < 0.) pt % mat (i,j) = 0.
                 IF ( pe % mat (i,j) < 0.) pe % mat (i,j) = 0.
                 IF ( pet % mat (i,j) < 0.) pet % mat (i,j) = 0.
                 
                 
            END IF
        END DO
    END DO
    
END IF

RETURN
END SUBROUTINE PETupdate



!==============================================================================
!| Description:
!   read new Kc values 
SUBROUTINE KcUpdate &
!
( time )

IMPLICIT NONE

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

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

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

SELECT CASE (interpolationMethodKc)
    
   CASE (0) !update map from netcdf
       CALL ReadField (filenameKc, time, dtKc, dtKc, 'M', kc_map)
    
   CASE (1) !update kc map from map of crop coefficient code 
       !read new data
       CALL ReadData (network = kc_stations, fileunit = fileunitKc, &
                      time = time, aggr_time = dtKc, &
                      aggr_type = 'ave')
	
	    !interpolation
	    DO i = 1, kc_map % idim
	      DO j = 1, kc_map % jdim
		      IF ( kc_map % mat(i,j) == kc_map % nodata ) CYCLE
			
                DO k = 1, kc_stations % countObs
			        IF ( ToString ( kc_code_map % mat(i,j) ) == &
                          kc_stations % obs(k) % id ) THEN
                        
				       kc_map % mat(i,j) = kc_stations % obs(k) % value
                       EXIT
                       
				    END IF
			    END DO
          END DO
	    END DO
END SELECT

RETURN

END SUBROUTINE KcUpdate


!==============================================================================
!| Description:
!   Compute potential evapotranspiration with Hargreaves Samani model
!   at daily time scale
!
! References:
!   Hargreaves, G.H., and Z.A. Samani (1982). Estimating potential 
!   evapotranspiration. J. Irrig. Drain. Eng., ASCE, 108(3), 223-230
!
!   Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop 
!      Evapotranspiration-Guidelines for Computing Crop Water Requirements-FAO 
!      Irrigation and Drainage Paper 56, 9th ed.; Food and Agriculture 
!      Organization of the United Nations: Rome, Italy, 1998; ISBN 92-5-104219-5.
SUBROUTINE HargreavesSamani &
!
(time, Tavg, Tmax, Tmin, etrad, pet)

IMPLICIT NONE

!Arguments with intent(in):
TYPE(DateTime)     , INTENT(in)  :: time !!simulation time
REAL (KIND = float), INTENT(in)  :: Tavg !!average daily temperature [°C]
REAL (KIND = float), INTENT(in)  :: Tmax !!maximum daily temperature [°C]
REAL (KIND = float), INTENT(in)  :: Tmin !!minimum daily temperature [°C]
REAL (KIND = float), INTENT(in)  :: etRad !!extra terrestrial solar radiation [mm/day]

!Arguments with intent(out):
REAL (KIND = float), INTENT(out) :: pet		!! potential evapotranspiration [m/s]

!local declarations.
REAL (KIND = float), PARAMETER   :: HC = 0.0023 !!empirical coefficient
REAL (KIND = float), PARAMETER   :: HE = 0.5  !!empirical exponent
REAL (KIND = float), PARAMETER   :: HT = 17.8  !!empirical temperature coefficient

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

!calculate evapotranspiration
IF ( ( Tmax - Tmin) < 0.) THEN
    CALL Catch ('warning', 'Evapotranspiration', '(Tmax - Tmin) < 0.' )
	pet = 0.
	RETURN
ELSE
	pet = HC * etRad * (Tmax - Tmin) ** HE * (Tavg + HT)
END IF

!conversion from mm/day to m/s
pet = pet * millimeter / day
	
    
RETURN	
END SUBROUTINE HargreavesSamani





!==============================================================================
!| Description:
!   Compute potential evapotranspiration at daily time scale with 
!   Hargreaves Samani model modified by Ravazzani et al. (2012)
!   to include an elevation correction factor
!
! References:
!   Hargreaves, G.H., and Z.A. Samani (1982). Estimating potential 
!   evapotranspiration. J. Irrig. Drain. Eng., ASCE, 108(3), 223-230
!
!   Ravazzani, G., Corbari, C., Morella, S., Gianoli, P., & 
!   Mancini, M.. (2012). Modified Hargreaves-Samani equation 
!   for the assessment of reference evapotranspiration in 
!   Alpine river basins. Journal of irrigation and drainage 
!   engineering, 138(7), 592–599.
SUBROUTINE HargreavesSamaniModified &
!
(time, Tavg, Tmax, Tmin, etrad, elevation, pet)

IMPLICIT NONE

!Arguments with intent(in):
TYPE(DateTime)     , INTENT(in)  :: time !!simulation time
REAL (KIND = float), INTENT(in)  :: Tavg !!average daily temperature [°C]
REAL (KIND = float), INTENT(in)  :: Tmax !!maximum daily temperature [°C]
REAL (KIND = float), INTENT(in)  :: Tmin !!minimum daily temperature [°C]
REAL (KIND = float), INTENT(in)  :: etRad !!extra terrestrial solar radiation [mm/day]
REAL (KIND = float), INTENT(in)  :: elevation !!elevation above sea level [m]


!Arguments with intent(out):
REAL (KIND = float), INTENT(out) :: pet		!! potential evapotranspiration [m/s]

!local declarations.
REAL (KIND = float), PARAMETER   :: HC = 0.0023 !!empirical coefficient
REAL (KIND = float), PARAMETER   :: HE = 0.5  !!empirical exponent
REAL (KIND = float), PARAMETER   :: HT = 17.8  !!empirical temperature coefficient
REAL (KIND = float), PARAMETER   :: c0 = 0.817  !!empirical coefficient for elevation correlation
REAL (KIND = float), PARAMETER   :: c1 = 0.00022   !!empirical coefficient for elevation correlation

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

!calculate evapotranspiration
IF ( ( Tmax - Tmin) < 0.) THEN
    CALL Catch ('warning', 'Evapotranspiration', '(Tmax - Tmin) < 0.' )
	pet = 0.
	RETURN
ELSE
	pet = (c0 + c1 * elevation) * HC * etRad * (Tmax - Tmin) ** HE * (Tavg + HT)
END IF

!conversion from mm/day to m/s
pet = pet * millimeter / day
	
    
RETURN	
END SUBROUTINE HargreavesSamaniModified





!==============================================================================
!| Description:
!   Compute potential evapotranspiration with Priestley-Taylor model
!
! References:
!   Priestley C.H.B., Taylor R.J., 1972. On the assessment
!   of surface heat flux and evaporation using largescale
!   parameters. Monthly Weather Review, 100:81-92.
!
!   Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop 
!      Evapotranspiration-Guidelines for Computing Crop Water Requirements-FAO 
!      Irrigation and Drainage Paper 56, 9th ed.; Food and Agriculture 
!      Organization of the United Nations: Rome, Italy, 1998; ISBN 92-5-104219-5.
!
!   Anderson, M. C., J. M. Norman, J. R. Mecikalski, J. P. Otkin, and W. P. 
!      Kustas (2007), A climatological study of evapotranspiration and moisture
!      stress across the continental U.S. based on the thermal remote sensing: I.
!      Model formulation, J. Geophys. Res., 112, D10117, doi:10.1029/2006JD007506
!
!  ASCE–EWRI. (2005). “The ASCE standardized reference evapotranspiration
!     equation.” ASCE–EWRI Standardization of Reference Evapotranspiration
!     Task Committe Rep., ASCE Reston, Va
!
SUBROUTINE PriestleyTaylor &
!
(airTemp, netRad, fc, elevation, pt, pe, pet)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(in)  :: airTemp !!air temperature [°C]
REAL (KIND = float), INTENT(in)  :: netRad !! net radiawion [W/m2]
REAL (KIND = float), INTENT(in)  :: fc !! fractional coverage by vegetation [0-1]
REAL (KIND = float), INTENT(in)  :: elevation !! terrain elevation [m a.s.l.]

!Arguments with intent(out):
REAL (KIND = float), INTENT(out) :: pt		!! potential transpiration (from vegetation) [m/s]
REAL (KIND = float), INTENT(out) :: pe		!! potential evaporation (from water or saturated soil) [m/s]
REAL (KIND = float), INTENT(out) :: pet		!! potential evapotranspiration [m/s]

!local declarations.
REAL (KIND = float), PARAMETER :: alpha = 1.26 !! priestely taylor advection coefficient 
REAL (KIND = float), PARAMETER :: toMJm2 = 0.000001 !! factor to convert W/me to MJ m-2 s-1
REAL (KIND = float), PARAMETER :: lambda = 2.453 !2.2647 !!latent heat of vaporization (MJ / kg)
REAL (KIND = float) :: netRadMJ !!net radiation flux in MJ / m2 /s 
REAL (KIND = float) :: es !!saturation vapor pressure 
REAL (KIND = float) :: des !!the slope of the relationship between saturation vapour pressure and ]
REAL (KIND = float) :: airPress !!air pressure [KPa]
REAL (KIND = float) :: gamma !!psychrometric constant [kPa °C-1]
REAL (KIND = float) :: groundHeat !!ground heat flux [MJ m-2 s-1]

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

! compute saturation vapor pressure
es = 0.6108 * EXP ( ( 17.27 * airTemp ) / ( airTemp + 237.3 ) )

!compute the slope of saturation vapour pressure curve (Allen. et al. 1998)
des = ( 4098. * es ) / ( ( airTemp + 237.3 )**2. )

!compute atmospheric pressure (ASCE-EWRI, 2005)
airPress = 101.3 * ( ( 293. - 0.0065 * elevation ) / 293. )** 5.26

!compute psychrometric constant
gamma = 0.665 * 0.001 * airPress

!compute neat radiation in MJ / m2 / s
netRadMJ = netRad * toMJm2

! compute ground heat flux as a fraction of net radiation Anderson et al. (2007)
! groundHeat is assumed to affect evaporation from bare soil
groundHeat = 0.1  * netRadMJ
 

!compute potential evaporation (from saturated bare soil or water)
pe = alpha * ( des / ( des + gamma ) ) * ( netRadMJ - groundHeat) / lambda * millimeter ! (m/s) 

!compute potential transpiration from vegetation
pt = alpha * ( des / ( des + gamma ) ) * netRadMJ / lambda * millimeter ! (m/s) 

!compute potential as a weighted average of bare soil and vegetated fluxes
pet = ( 1. - fc) * pe + fc * pt

RETURN
END SUBROUTINE PriestleyTaylor




!==============================================================================
!| Description:
!   Compute potential evapotranspiration with Penman-Monteith model
!
!   Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop 
!      Evapotranspiration-Guidelines for Computing Crop Water Requirements-FAO 
!      Irrigation and Drainage Paper 56, 9th ed.; Food and Agriculture 
!      Organization of the United Nations: Rome, Italy, 1998; ISBN 92-5-104219-5.
!
!     ASCE–EWRI. 
2005. “The ASCE standardized reference evapotranspiration
!     equation.” ASCE–EWRI Standardization of Reference Evapotranspiration
!     Task Committe Rep., ASCE Reston, Va
!
! https://wetlandscapes.github.io/blog/blog/penman-monteith-and-priestley-taylor/
!
SUBROUTINE PenmanMonteith &
!
(airTemp, netRad, rh, wind, fc, z, elevation, zws, zrhum, lai, srmin, pt, pe, pet)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(in)  :: airTemp !!air temperature [°C]
REAL (KIND = float), INTENT(in)  :: netRad !! net radiawion [W/m2]
REAL (KIND = float), INTENT(in)  :: rh !! air relative humidity [0-100]
REAL (KIND = float), INTENT(in)  :: wind !! wind speed [m/s]
REAL (KIND = float), INTENT(in)  :: fc !! fractional coverage by vegetation [0-1]
REAL (KIND = float), INTENT(in)  :: z !! vegetation height [m]
REAL (KIND = float), INTENT(in)  :: elevation !! terrain elevation [m a.s.l.]
REAL (KIND = float), INTENT(in)  :: zws !! wind speed measurement heigth [m]
REAL (KIND = float), INTENT(in)  :: zrhum !! relative humidity measurement heigth [m]
REAL (KIND = float), INTENT(in)  :: lai !! leaf area index [m2/m2]
REAL (KIND = float), INTENT(in)  :: srmin !! minimum stomatal resistance

!Arguments with intent(out):
REAL (KIND = float), INTENT(out) :: pt		!! potential transpiration (from vegetation) [m/s]
REAL (KIND = float), INTENT(out) :: pe		!! potential evaporation (from water or saturated soil) [m/s]
REAL (KIND = float), INTENT(out) :: pet		!! potential evapotranspiration [m/s]

!local declarations:
REAL (KIND = float) :: groundHeat !!ground heat flux [MJ m-2 s-1]
REAL (KIND = float) :: es !!saturation vapor pressure (Pa)
REAL (KIND = float) :: ea !!actual vapor pressure (Pa)
REAL (KIND = float) :: des !!the slope of the relationship between saturation vapour pressure and temperature (kPa°C-1)
REAL (KIND = float) :: airPress !!air pressure [KPa]
REAL (KIND = float) :: gamma !!psychrometric constant [kPa °C-1]
REAL (KIND = float) :: ws !!wind speed (m/s)
REAL (KIND = float) :: z1  !!bare soil heigth (m)
REAL (KIND = float) :: zm1  !! wind speed measurement heigth + z1 (m)
REAL (KIND = float) :: zh1  !! relative humidity measurement heigth + z1 (m)
REAL (KIND = float) :: zd1   !! bare soil zero plane displacement height (m)
REAL (KIND = float) :: zrw1  !! Roughness length governing momentum transfer above zd1 (m)
REAL (KIND = float) :: zrh1   !! Roughness length governing transfer of heat and vapor above zd1 (m)
REAL (KIND = float) :: rabs !! aerodynamic resistance of bare soil (s/m)
REAL (KIND = float) :: zm   !! wind speed measurement heigth + z (m)
REAL (KIND = float) :: zd    !! vegetated soil zero plane displacement height (m)
REAL (KIND = float) :: zrw   !! Roughness length governing momentum transfer above zd (m)
REAL (KIND = float) :: zh  !! relative humidity measurement heigth + z (m)
REAL (KIND = float) :: zrh   !! Roughness length governing transfer of heat and vapor above zd (m)
REAL (KIND = float) :: ra !! aerodynamic resistance of vegetated soil (s/m)
REAL (KIND = float) :: cp !!  humid air specific heat (MJkg^-1°C^-1)
REAL (KIND = float) :: Tkv  !! virtual air temperature (K)
REAL (KIND = float) :: rhoAir !!air density (kgm^-3)
REAL (KIND = float) :: netRadMJ !! net radiation in MJm^-2s^-1
REAL (KIND = float) :: rc !! total vegetation  resistance (s/m)
REAL (KIND = float), PARAMETER :: k = 0.41 !!von Karman’s constant
REAL (KIND = float), PARAMETER :: epsilon = 0.622   ! ratio molecular weight of water vapour/dry air
REAL (KIND = float), PARAMETER :: lambda = 2.453 !2.2647 !!latent heat of vaporization (MJ / kg)
REAL (KIND = float), PARAMETER :: R1 = 0.287	!!gas specific constant	(kJkg^-1K^-1)
!-----------------------------------------end of declarations------------------

! compute saturation vapor pressure
es = 0.6108 * EXP ( ( 17.27 * airTemp ) / ( airTemp + 237.3 ) )

!actual vapor pressure
ea = rh / 100. * es

!compute the slope of saturation vapour pressure curve (Allen. et al. 1998)
des = ( 4098. * es ) / ( ( airTemp + 237.3 )**2. )

!compute atmospheric pressure (ASCE-EWRI, 2005)
airPress = 101.3 * ( ( 293. - 0.0065 * elevation ) / 293. )** 5.26

!compute psychrometric constant
gamma = 0.665 * 0.001 * airPress

!check wind speed: set minimum to 0.01 m/s
IF ( wind <= 0 ) THEN
	ws = 0.01
ELSE
	ws = wind
END IF

!aerodynamic variables of bare soil
z1 = 0.1 
zm1 = zws + z1
zd1 = 2. / 3. * z1
zrw1 = 0.123 * z1
zh1 = zrhum + z1
zrh1 = 0.1 * zrw1

!compute aerodynamic resistance of bare soil
rabs = LOG ( (zm1 - zd1) / zrw1 ) * LOG ( (zh1 - zd1) / zrh1 ) / ( k**2 * ws ) 

!aerodynamic variables of vegetated soil
zm = zws + z
zd = 2. / 3. * z
zrw = 0.123 * z
zh = zrhum + z
zrh = 0.1 * zrw

!compute aerodynamic resistance of bare soil
IF (z == 0 ) THEN  !vegetation not present
	ra = rabs
ELSE
	ra = LOG ( ( zm - zd ) / zrw ) * LOG ( ( zh - zd ) / zrh ) / ( k**2 * ws)
END IF

! humid air specific heat (MJkg^-1°C^-1)
cp = gamma * epsilon * lambda / airPress

!virtual air temperature (K)
Tkv = 1.01 * ( airTemp + 273. )

!air density (kgm^-3)
rhoAir = airPress / ( Tkv * R1 )


!check and convert net radiation
IF ( netRad < 0 ) THEN
	netRadMJ = 0.
ELSE
	netRadMJ = netRad * 0.000001
END IF 

! compute ground heat flux as a fraction of net radiation Anderson et al. (2007)
! groundHeat is assumed to affect evaporation from bare soil
!groundHeat = 0.31 * ( 1. - fc ) * netRadMJ
groundHeat = 0.1 * netRadMJ 


!compute potential evaporation (from saturated bare soil or water)
pe = ( des * ( netRadMJ - groundHeat ) + rhoAir * cp * ( es - ea ) / rabs  )  / ( lambda *  ( des + gamma )  ) * millimeter
    

!compute potential transpiration from vegetation
IF (fc > 0.) THEN  
    IF ( lai == 0.) THEN
        rc = srmin / 1.
    ELSE
        rc = srmin / lai
    END IF
    
    pt = ( des * netRadMJ + rhoAir * cp * ( es - ea ) / ra  )  / ( lambda *  ( des + gamma * (1. + rc / ra) )  ) * millimeter

ELSE
    
    pt = 0.

END IF

!compute potential as a weighted average of bare soil and vegetated fluxes
pet = ( 1. - fc) * pe + fc * pt

!IF (isnan (pet)) then
!    write(*,*) pe, pt, des, ra, rc, lai, srmin
!    pause
!end if

RETURN
END SUBROUTINE PenmanMonteith


!==============================================================================
!| Description:
!   Compute potential evapotranspiration with FAO56 Penman-Monteith model
!   for a reference grass
!   "A hypothetical reference crop with an assumed crop height of 0.12 m, 
!    a fixed surface resistance of 70 s m-1 and an albedo of 0.23."
!
! References:
!
!   Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop 
!      Evapotranspiration-Guidelines for Computing Crop Water Requirements-FAO 
!      Irrigation and Drainage Paper 56, 9th ed.; Food and Agriculture 
!      Organization of the United Nations: Rome, Italy, 1998; ISBN 92-5-104219-5.
!
!   Lincoln Zotarelli, Michael D. Dukes, Consuelo C. Romero, 
!     Kati W. Migliaccio, and Kelly T. Morgan, Step by Step Calculation 
!     of the Penman-Monteith Evapotranspiration (FAO-56 Method), 
!     Agricultural and Biological Engineering Department, University of Florida,
!     2009. https://edis.ifas.ufl.edu/pdf/AE/AE45900.pdf
!
!   https://wetlandscapes.github.io/blog/blog/penman-monteith-and-priestley-taylor/
!
SUBROUTINE FAO56PenmanMonteith &
!
(airTemp, netRad, netRadFAO, rh, wind, fc, elevation, zws, zrhum, lai, pt, pe, pet)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(in)  :: airTemp !!air temperature [°C]
REAL (KIND = float), INTENT(in)  :: netRad !! net radiation [W/m2]
REAL (KIND = float), INTENT(in)  :: netRadFAO !! net radiation for reference vegetation, with albedo = 0.23 [W/m2]
REAL (KIND = float), INTENT(in)  :: rh !! air relative humidity [0-100]
REAL (KIND = float), INTENT(in)  :: wind !! wind speed [m/s]
REAL (KIND = float), INTENT(in)  :: fc !! fractional coverage by vegetation [0-1]

REAL (KIND = float), INTENT(in)  :: elevation !! terrain elevation [m a.s.l.]
REAL (KIND = float), INTENT(in)  :: zws !! wind speed measurement heigth [m]
REAL (KIND = float), INTENT(in)  :: zrhum !! relative humidity measurement heigth [m]
REAL (KIND = float), INTENT(in)  :: lai !! leaf area index [m2/m2]


!Arguments with intent(out):
REAL (KIND = float), INTENT(out) :: pt		!! potential transpiration (from vegetation) [m/s]
REAL (KIND = float), INTENT(out) :: pe		!! potential evaporation (from water or saturated soil) [m/s]
REAL (KIND = float), INTENT(out) :: pet		!! potential evapotranspiration [m/s]

!local declarations:
REAL (KIND = float) :: z = 0.12 !! reference grass  height [m]
REAL (KIND = float) :: srmin = 70 !! minimum stomatal resistance of reference grass [s/m]

REAL (KIND = float) :: groundHeat !!ground heat flux [MJ m-2 s-1]
REAL (KIND = float) :: es !!saturation vapor pressure (Pa)
REAL (KIND = float) :: ea !!actual vapor pressure (Pa)
REAL (KIND = float) :: des !!the slope of the relationship between saturation vapour pressure and temperature
REAL (KIND = float) :: airPress !!air pressure [KPa]
REAL (KIND = float) :: gamma !!psychrometric constant [kPa °C-1]
REAL (KIND = float) :: ws2m !!wind speed at 2m (m/s)
REAL (KIND = float) :: z1  !!bare soil heigth (m)
REAL (KIND = float) :: zm1  !! wind speed measurement heigth + z1 (m)
REAL (KIND = float) :: zh1  !! relative humidity measurement heigth + z1 (m)
REAL (KIND = float) :: zd1   !! bare soil zero plane displacement height (m)
REAL (KIND = float) :: zrw1  !! Roughness length governing momentum transfer (m)
REAL (KIND = float) :: zrh1   !! Roughness length governing transfer of heat and vapor above zd1 (m)
REAL (KIND = float) :: rabs !! aerodynamic resistance of bare soil (s/m)
REAL (KIND = float) :: cp !!  humid air specific heat (MJkg^-1°C^-1)
REAL (KIND = float) :: Tkv  !! virtual air temperature (K)
REAL (KIND = float) :: rhoAir !!air density (kgm^-3)
REAL (KIND = float) :: netRadMJ !! net radiation in MJm^-2s^-1
REAL (KIND = float) :: netRadFAO_MJ !! FAO net radiation in MJm^-2s^-1
REAL (KIND = float), PARAMETER :: k = 0.41 !!von Karman’s constant
REAL (KIND = float), PARAMETER :: epsilon = 0.622   ! ratio molecular weight of water vapour/dry air
REAL (KIND = float), PARAMETER :: lambda = 2.453 !2.2647 !!latent heat of vaporization (MJ / kg)
REAL (KIND = float), PARAMETER :: R1 = 0.287	!!gas specific constant	(kJkg^-1K^-1)
!-----------------------------------------end of declarations------------------

! compute saturation vapor pressure
es = 0.6108 * EXP ( ( 17.27 * airTemp ) / ( airTemp + 237.3 ) )

!actual vapor pressure
ea = rh / 100. * es

!compute the slope of saturation vapour pressure curve (Allen. et al. 1998)
des = ( 4098. * es ) / ( ( airTemp + 237.3 )**2. )

!compute atmospheric pressure
airPress = 101.3 * ( ( 293. - 0.0065 * elevation ) / 293. )** 5.26

!compute psychrometric constant
gamma = 0.665 * 0.001 * airPress


!move wind speed to 2 m reference height (Zotarelli et al., 2009)
IF ( zws /= 2. ) THEN
	ws2m = wind
ELSE
	ws2m = wind * 4.87 / LOG ( 67.8 * zws - 5.42 )
END IF

!aerodynamic variables of bare soil
z1 = 0.1 
zm1 = zws + z1
zd1 = 2. / 3. * z1
zrw1 = 0.123 * z1
zh1 = zrhum + z1
zrh1 = 0.1 * zrw1

!compute aerodynamic resistance of bare soil
rabs = LOG ( (zm1 - zd1) / zrw1 ) * LOG ( (zh1 - zd1) / zrh1 ) / ( k**2 * wind ) 

! humid air specific heat (MJkg^-1°C^-1)
cp = gamma * epsilon * lambda / airPress

!virtual air temperature (K)
Tkv = 1.01 * ( airTemp + 273. )

!air density (kgm^-3)
rhoAir = airPress / ( Tkv * R1 )


!check and convert net radiation
IF ( netRad < 0 ) THEN
	netRadMJ = 0.
ELSE
	netRadMJ = netRad * 0.000001
END IF 

IF ( netRadFAO < 0 ) THEN
	netRadFAO_MJ = 0.
ELSE
	netRadFAO_MJ = netRadFAO * 0.000001
END IF 

! compute ground heat flux as a fraction of net radiation Anderson et al. (2007)
! groundHeat is assumed to affect evaporation from bare soil only
groundHeat = 0.31 * netRadMJ

!compute evaporation (from saturated bare soil or water)
    
pe = ( des * ( netRadMJ - groundHeat ) + rhoAir * cp * ( es - ea ) / rabs  )  / &
     ( lambda *  ( des + gamma )  ) * millimeter
    

!compute transpiration from vegetation
IF (fc > 0.) THEN  

    pt =  ( 0.408 * des * ( netRadMJ - groundHeat ) + gamma *  (900. / day ) / &
          ( airTemp + 273. ) * ws2m * (es - ea) ) / &
          (des + gamma * ( 1. + 0.34 * ws2m ) ) * millimeter

ELSE
    
    pt = 0.

END IF

!compute evapotranspiration

pet = ( 1. - fc) * pe + fc * pt

RETURN
END SUBROUTINE FAO56PenmanMonteith








!==============================================================================
!| Description:
!   Compute actual evapotranspiration from ground 
!   by solving energy balance
!   Original code written by Chiara Corbari and Jacopo Martinelli
!
SUBROUTINE etEnergyBalanceGround &
!
(airTemp, netRad, rh, wind, fc, z, elevation, zws, zrhum, lai, srmin, pt, pe, pet)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(in)  :: airTemp !!air temperature (°C)
REAL (KIND = float), INTENT(in)  :: netRad !! net radiawion (W/m2)
REAL (KIND = float), INTENT(in)  :: rh !! air relative humidity (0-100)
REAL (KIND = float), INTENT(in)  :: wind !! wind speed(m/s)
REAL (KIND = float), INTENT(in)  :: fc !! fractional coverage by vegetation (0-1)
REAL (KIND = float), INTENT(in)  :: z !! vegetation height (m)
REAL (KIND = float), INTENT(in)  :: elevation !! terrain elevation (m a.s.l.)
REAL (KIND = float), INTENT(in)  :: zws !! wind speed measurement heigth (m)
REAL (KIND = float), INTENT(in)  :: zrhum !! relative humidity measurement heigth (m)
REAL (KIND = float), INTENT(in)  :: lai !! leaf area index (m2/m2)
REAL (KIND = float), INTENT(in)  :: srmin !! minimum stomatal resistance

!Arguments with intent(out):
REAL (KIND = float), INTENT(out) :: pt		!! potential transpiration (from vegetation) (m/s)
REAL (KIND = float), INTENT(out) :: pe		!! potential evaporation (from water or saturated soil) (m/s)
REAL (KIND = float), INTENT(out) :: pet		!! potential evapotranspiration (m/s)

!local declarations:
REAL (KIND = float) :: groundHeat !!ground heat flux (MJ m-2 s-1)
REAL (KIND = float) :: es !!saturation vapor pressure (Pa)
REAL (KIND = float) :: ea !!actual vapor pressure (Pa)
REAL (KIND = float) :: des !!the slope of the relationship between saturation vapour pressure and temperature
REAL (KIND = float) :: airPress !!air pressure (KPa)
REAL (KIND = float) :: gamma !!psychrometric constant (kPa °C-1)
REAL (KIND = float) :: ws !!wind speed (m/s)
REAL (KIND = float) :: z1  !!bare soil heigth (m)
REAL (KIND = float) :: zm1  !! wind speed measurement heigth + z1 (m)
REAL (KIND = float) :: zh1  !! relative humidity measurement heigth + z1 (m)
REAL (KIND = float) :: zd1   !! bare soil zero plane displacement height (m)
REAL (KIND = float) :: zrw1  !! Roughness length governing momentum transfer above zd1 (m)
REAL (KIND = float) :: zrh1   !! Roughness length governing transfer of heat and vapor above zd1 (m)
REAL (KIND = float) :: rabs !! aerodynamic resistance of bare soil (s/m)
REAL (KIND = float) :: zm   !! wind speed measurement heigth + z (m)
REAL (KIND = float) :: zd    !! vegetated soil zero plane displacement height (m)
REAL (KIND = float) :: zrw   !! Roughness length governing momentum transfer above zd (m)
REAL (KIND = float) :: zh  !! relative humidity measurement heigth + z (m)
REAL (KIND = float) :: zrh   !! Roughness length governing transfer of heat and vapor above zd (m)
REAL (KIND = float) :: ra !! aerodynamic resistance of vegetated soil (s/m)
REAL (KIND = float) :: cp !!  humid air specific heat (MJkg^-1°C^-1)
REAL (KIND = float) :: Tkv  !! virtual air temperature (K)
REAL (KIND = float) :: rhoAir !!air density (kgm^-3)
REAL (KIND = float) :: netRadMJ !! net radiation in MJm^-2s^-1
REAL (KIND = float) :: rc !! total vegetation  resistance (s/m)
REAL (KIND = float), PARAMETER :: k = 0.41 !!von Karman’s constant
REAL (KIND = float), PARAMETER :: epsilon = 0.622   ! ratio molecular weight of water vapour/dry air
REAL (KIND = float), PARAMETER :: lambda = 2.453 !2.2647 !!latent heat of vaporization (MJ / kg)
REAL (KIND = float), PARAMETER :: R1 = 0.287	!!gas specific constant	(kJkg^-1K^-1)
!-----------------------------------------end of declarations------------------

! compute saturation vapor pressure
es = 0.6108 * EXP ( ( 17.27 * airTemp ) / ( airTemp + 237.3 ) )

!actual vapor pressure
ea = rh / 100. * es

!compute the slope of saturation vapour pressure curve (Allen. et al. 1998)
des = ( 4098. * es ) / ( ( airTemp + 237.3 )**2. )

!compute atmospheric pressure (ASCE-EWRI, 2005)
airPress = 101.3 * ( ( 293. - 0.0065 * elevation ) / 293. )** 5.26

!compute psychrometric constant
gamma = 0.665 * 0.001 * airPress

!check wind speed: set minimum to 0.01 m/s
IF ( wind <= 0 ) THEN
	ws = 0.01
ELSE
	ws = wind
END IF

!aerodynamic variables of bare soil
z1 = 0.1 
zm1 = zws + z1
zd1 = 2. / 3. * z1
zrw1 = 0.123 * z1
zh1 = zrhum + z1
zrh1 = 0.1 * zrw1

!compute aerodynamic resistance of bare soil
rabs = LOG ( (zm1 - zd1) / zrw1 ) * LOG ( (zh1 - zd1) / zrh1 ) / ( k**2 * ws ) 

!aerodynamic variables of vegetated soil
zm = zws + z
zd = 2. / 3. * z
zrw = 0.123 * z
zh = zrhum + z
zrh = 0.1 * zrw

!compute aerodynamic resistance of bare soil
IF (z == 0 ) THEN  !vegetation not present
	ra = rabs
ELSE
	ra = LOG ( ( zm - zd ) / zrw ) * LOG ( ( zh - zd ) / zrh ) / ( k**2 * ws)
END IF

! humid air specific heat (MJkg^-1°C^-1)
cp = gamma * epsilon * lambda / airPress

!virtual air temperature (K)
Tkv = 1.01 * ( airTemp + 273. )

!air density (kgm^-3)
rhoAir = airPress / ( Tkv * R1 )


!check and convert net radiation
IF ( netRad < 0 ) THEN
	netRadMJ = 0.
ELSE
	netRadMJ = netRad * 0.000001
END IF 

! compute ground heat flux as a fraction of net radiation Anderson et al. (2007)
! groundHeat is assumed to affect evaporation from bare soil
!groundHeat = 0.31 * ( 1. - fc ) * netRadMJ
groundHeat = 0.1 * netRadMJ 


!compute potential evaporation (from saturated bare soil or water)
pe = ( des * ( netRadMJ - groundHeat ) + rhoAir * cp * ( es - ea ) / rabs  )  / ( lambda *  ( des + gamma )  ) * millimeter
    

!compute potential transpiration from vegetation
IF (fc > 0.) THEN  
    IF ( lai == 0.) THEN
        rc = srmin / 1.
    ELSE
        rc = srmin / lai
    END IF
    
    pt = ( des * netRadMJ + rhoAir * cp * ( es - ea ) / ra  )  / ( lambda *  ( des + gamma * (1. + rc / ra) )  ) * millimeter

ELSE
    
    pt = 0.

END IF

!compute potential as a weighted average of bare soil and vegetated fluxes
pet = ( 1. - fc) * pe + fc * pt

!IF (isnan (pet)) then
!    write(*,*) pe, pt, des, ra, rc, lai, srmin
!    pause
!end if

RETURN
END SUBROUTINE etEnergyBalanceGround









!==============================================================================
!| Description:
!   check properties  for each ET model
SUBROUTINE CheckSpecificProperties &
!
( model, ini )

IMPLICIT NONE

!Arguments with intent(in):
INTEGER (KIND = short), INTENT(in) :: model
TYPE (IniList)        , INTENT(in) :: ini

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

!check data specific for each ET model
SELECT CASE (model)
   CASE ( PENMAN_MONTEITH ) 
      IF ( dtTemperature <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires air temperature' )
      END IF
      
      IF ( dtRelHumidity <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires air relative humidity' )
      END IF
      
      IF ( dtWindSpeed <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires wind speed' )
      END IF
      
      IF ( dtRadiation <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires net radiation' )
      END IF
      
      IF ( .NOT. fvcoverLoaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires vegetation fraction coverage' )
      END IF
      
      IF ( .NOT. laiLoaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires leaf area index' )
      END IF
      
      IF ( .NOT. dem_loaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires digital elevation model' )
      END IF
      
      IF ( .NOT. rsMinloaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires minimum stomatal resistance' )
      END IF
      
      
   CASE ( FAO56_PENMAN_MONTEITH ) 
      IF ( dtTemperature <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires air temperature' )
      END IF
      
      IF ( dtRelHumidity <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires air relative humidity' )
      END IF
      
      IF ( dtWindSpeed <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires wind speed' )
      END IF
      
      IF ( dtRadiation <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires net radiation' )
      END IF
      
      IF ( .NOT. fvcoverLoaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires vegetation fraction coverage' )
      END IF
      
      
      IF ( .NOT. dem_loaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires digital elevation model' )
      END IF
      
      
   CASE (PRIESTLEY_TAYLOR) 
       IF ( dtTemperature <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires air temperature' )
      END IF
      
      IF ( dtRadiation <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires net radiation' )
      END IF
      
      IF ( .NOT. fvcoverLoaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires vegetation fraction coverage' )
      END IF
      
      IF ( .NOT. dem_loaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires digital elevation model' )
      END IF
      
      !IF ( .NOT. plantsHeightLoaded ) THEN
      !     CALL Catch ('error', 'Evapotranspiration',   &
      !          'ET model requires vegetation height' )
      !END IF
 
   CASE (HARGREAVES) 
       
       compute_hargreaves = .TRUE.
      
       
      IF ( dtET /= day ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires daily dt' ) 
      END IF
       
      IF ( dtTemperatureDailyMean <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires daily air temperature' )
      END IF
      
      IF ( dtTemperatureDailyMin <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires daily minimum air temperature' )
      END IF
      
      IF ( dtTemperatureDailyMax <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires daily maximum air temperature' )
      END IF
      
   CASE (HARGREAVES_MOD) 
      
      compute_hargreaves = .TRUE. 
       
      IF ( dtET /= day ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires daily dt' ) 
      END IF
       
      IF ( dtTemperatureDailyMean <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires daily air temperature' )
      END IF
      
      IF ( dtTemperatureDailyMin <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires daily minimum air temperature' )
      END IF
      
      IF ( dtTemperatureDailyMax <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires daily maximum air temperature' )
      END IF
      
      IF ( .NOT. dem_loaded ) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires digital elevation model' )
      END IF
      
      
      
        
   CASE (ENERGY_BALANCE)
       
       IF ( dtTemperature <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires air temperature' )
      END IF
      
      IF ( dtRelHumidity <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires air relative humidity' )
      END IF
      
      IF ( dtWindSpeed <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires wind speed' )
      END IF
      
      IF ( dtRadiation <= 0) THEN
          CALL Catch ('error', 'Evapotranspiration',   &
                'ET model requires net radiation' )
      END IF
      
	
END SELECT


RETURN
END SUBROUTINE CheckSpecificProperties



END MODULE Evapotranspiration