      Subroutine type224       
!**********************************************************************************************
! Object: Evapotranspiration model
! Simulation Studio Model: type224
! 

! Author: S.Alkhatib
! Editor: 
! Date:	 June 23, 2018
! last modified: June 23, 2018
! 
! 
! *** 
! *** Model Parameters 
! *** 
!			Area	m^2 [0;+Inf]
!			Plant's height	m [0;+Inf]
!			Minimum Leaf resistance	any [0;+Inf]
!			Height of wind measurements	m [0;+Inf]
!			Height of humidity measurements	m [0;+Inf]
!			Coef. J-S for radiation. S1	W/m^2 [0;+Inf]
!			Coef. J-S for radiation. S2	W/m^2 [0;+Inf]
!			Coef. J-S for temperature. t1	C [0;+Inf]
!			Coef. J-S for temperature. t2	C [0;+Inf]
!			Coef. J-S for temperature. t3	C [0;+Inf]
!			Coef. J-S for vapour pressure deficit. V1	kPa [0;+Inf]
!			Coef. J-S for vapour pressure deficit. V2	kPa [0;+Inf]
!			Coef. J-S for vapour pressure deficit. V3	kPa [0;+Inf]

! *** 
! *** Model Inputs 
! *** 
!			Air temperature	C [-Inf;+Inf]
!			Relative Humidity	any [0;1]
!			Air pressure	kPa [0;+Inf]
!			Air density	kg/m^3 [0;+Inf]
!			Air cp	J/kg.K [0;+Inf]
!			Air velocity	m/s [0;+Inf]
!			Irradiance	W/m^2 [0;+Inf]
!			Interception factor	any [0;1]

! *** 
! *** Model Outputs 
! *** 
!			Evapotranspiration Rate	kg/hr [0;+Inf]
!			Aerodynamic resistance	any [-Inf;+Inf]
!			Plant canopy resistance	any [-Inf;+Inf]
!			QEvapotranspiration	kJ/hr [-Inf;+Inf]

! *** 
! *** Model Derivatives 
! *** 

! (Comments and routine interface generated by TRNSYS Studio)
!************************************************************************

!export this subroutine for its use in external DLLs.
!DEC$ATTRIBUTES DLLEXPORT :: TYPE224

!-----------------------------------------------------------------------------------------------------------------------
      Use TrnsysConstants
      Use TrnsysFunctions
!-----------------------------------------------------------------------------------------------------------------------
!Variable Declarations
      Implicit None  !force explicit declaration of local variables
      Double Precision  Timestep,Time
      Integer CurrentUnit,CurrentType
      DOUBLE PRECISION area !Area [m^2]
      DOUBLE PRECISION h_c  !Plant_s_height [m]
      DOUBLE PRECISION r_l_min !Minimum_Leaf_resistance [-]
      DOUBLE PRECISION z_m     !Height_of_wind_measurements[m]
      DOUBLE PRECISION z_h     !Height_of_humidity_measurements [m]
      DOUBLE PRECISION Coef_S1     !Coef__J_S_for_radiation__S1 [W/m^2]
      DOUBLE PRECISION Coef_S2     !Coef__J_S_for_radiation__S2 [W/m^2]
      DOUBLE PRECISION Coef_t1     !Coef__J_S_for_temperature__t1 [C]
      DOUBLE PRECISION Coef_t2     !Coef__J_S_for_temperature__t2 [C]
      DOUBLE PRECISION Coef_t3     !Coef__J_S_for_temperature__t3 [C]
      DOUBLE PRECISION Coef_V1     !Coef__J_S_for_vapour_pressure_deficit__V1 [KPa]
      DOUBLE PRECISION Coef_V2     !Coef__J_S_for_vapour_pressure_deficit__V2 [KPa]
      DOUBLE PRECISION Coef_V3     !Coef__J_S_for_vapour_pressure_deficit__V3 [KPa]
      DOUBLE PRECISION e_sat       ! Mean saturation Vapor Pressure for a given temperature [KPa]
      DOUBLE PRECISION e_i         ! Actual vapor pressure derived from RH data [KPa]
      DOUBLE PRECISION delta_e     ! Vapor pressure defict VPD
      DOUBLE PRECISION DELTA       !rate of change (slope) of the saturation vapor pressure - temperature
      DOUBLE PRECISION lambda      !latent heat
      DOUBLE PRECISION gamma       !phychrometric constant
      DOUBLE PRECISION z_om        !roughness lenght governing momentum transfer
      DOUBLE PRECISION z_oh        !roughness lenght governing transfer of heat and vapour
      DOUBLE PRECISION k           !Karman's constant
      DOUBLE PRECISION v_a_in
      DOUBLE PRECISION r_a         !Aerodynamic resistance
      DOUBLE PRECISION LAI         !leaf area index [m^2]
      DOUBLE PRECISION LAI_act     !activate leaf area index [m^2]
      DOUBLE PRECISION F1
      DOUBLE PRECISION F2
      DOUBLE PRECISION F3
      DOUBLE PRECISION F4
      DOUBLE PRECISION F5
      DOUBLE PRECISION F6
      DOUBLE PRECISION d        !zero plane displacement height
      DOUBLE PRECISION r_l      !leaf resistance
      DOUBLE PRECISION r_c      !plant canopy resistance
      DOUBLE PRECISION ET0 !Massflow rate of Evapotranspiration [kg/hr] 
      DOUBLE PRECISION OUT1
      DOUBLE PRECISION OUT2
      DOUBLE PRECISION OUT3
      DOUBLE PRECISION QEVAP ! Evapotranspiration [kJ/hr]
      DOUBLE PRECISION T_a   !Air_temperature [C]
      DOUBLE PRECISION rh_a  !Relative_Humidity [%]
      DOUBLE PRECISION P_a   !Air_pressure [KPa]
      DOUBLE PRECISION rho_a !Air_density [kg/m^3]
      DOUBLE PRECISION cp_a  !Air_cp [J/kg.K]
      DOUBLE PRECISION v_a   !Air_velocity [m/s]
      DOUBLE PRECISION G     !Irradiance [W/m^2]
      DOUBLE PRECISION f     !Interception_factor [-]
!-----------------------------------------------------------------------------------------------------------------------
!Get the Global Trnsys Simulation Variables
      Time=getSimulationTime()
      Timestep=getSimulationTimeStep()
      CurrentUnit=getCurrentUnit()
      CurrentType=getCurrentType()
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Set the Version Number for This Type
      If(getIsVersionSigningTime()) Then
		Call SetTypeVersion(17)
		Return
      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Do Any Last Call Manipulations Here
      If(getIsLastCallofSimulation()) Then
		Return
      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Perform Any "After Convergence" Manipulations That May Be Required 
      If(getIsEndOfTimestep()) Then
		Return
      EndIf
!-----------------------------------------------------------------------------------------------------------------------

      
!-----------------------------------------------------------------------------------------------------------------------
!Do All of the "Very First Call of the Simulation Manipulations" Here
      If(getIsFirstCallofSimulation()) Then

		!Tell the TRNSYS Engine How This Type Works
		Call SetNumberofParameters(13)           !The number of parameters that the the model wants
		Call SetNumberofInputs(8)                   !The number of inputs that the the model wants
		Call SetNumberofDerivatives(0)         !The number of derivatives that the the model wants
		Call SetNumberofOutputs(4)                 !The number of outputs that the the model produces
		Call SetIterationMode(1)                             !An indicator for the iteration mode (default=1).  Refer to section 8.4.3.5 of the documentation for more details.
		Call SetNumberStoredVariables(0,0)                   !The number of static variables that the model wants stored in the global storage array and the number of dynamic variables that the model wants stored in the global storage array
		Call SetNumberofDiscreteControls(0)               !The number of discrete control functions set by this model (a value greater than zero requires the user to use Solver 1: Powell's method)

		Return

      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Do All of the First Timestep Manipulations Here - There Are No Iterations at the Intial Time
      If (getIsStartTime()) Then
      area = getParameterValue(1)
      h_c = getParameterValue(2)
      r_l_min = getParameterValue(3)
      z_m = getParameterValue(4)
      z_h = getParameterValue(5)
      Coef_S1 = getParameterValue(6)
      Coef_S2 = getParameterValue(7)
      Coef_t1 = getParameterValue(8)
      Coef_t2 = getParameterValue(9)
      Coef_t3 = getParameterValue(10)
      Coef_V1 = getParameterValue(11)
      Coef_V2 = getParameterValue(12)
      Coef_V3 = getParameterValue(13)
      
      T_a = GetInputValue(1)
      rh_a = GetInputValue(2)
      P_a = GetInputValue(3)
      rho_a = GetInputValue(4)
      cp_a = GetInputValue(5)
      v_a = GetInputValue(6)
      G = GetInputValue(7)
      f = GetInputValue(8)
   !Set the Initial Values of the Outputs (#,Value)
		Call SetOutputValue(1, 0) ! Evapotranspiration Rate
		Call SetOutputValue(2, 0) ! Aerodynamic resistance
		Call SetOutputValue(3, 0) ! Plant canopy resistance
		Call SetOutputValue(4, 0) ! QEvapotranspiration
		Return

      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!ReRead the Parameters if Another Unit of This Type Has Been Called Last
      If(getIsReReadParameters()) Then
          !Read in the Values of the Parameters from the Input File
      area = getParameterValue(1)
      h_c  = getParameterValue(2)
      r_l_min = getParameterValue(3)
      z_m = getParameterValue(4)
      z_h = getParameterValue(5)
      Coef_S1 = getParameterValue(6)
      Coef_S2 = getParameterValue(7)
      Coef_t1 = getParameterValue(8)
      Coef_t2 = getParameterValue(9)
      Coef_t3 = getParameterValue(10)
      Coef_V1 = getParameterValue(11)
      Coef_V2 = getParameterValue(12)
      Coef_V3 = getParameterValue(13)
      EndIf
!-----------------------------------------------------------------------------------------------------------------------
!Read the Inputs
      T_a = GetInputValue(1)
      rh_a = GetInputValue(2)
      P_a = GetInputValue(3)
      rho_a = GetInputValue(4)
      cp_a = GetInputValue(5)
      v_a = GetInputValue(6)
      G = GetInputValue(7)
      f = GetInputValue(8)
      If(ErrorFound()) Return
!-----------------------------------------------------------------------------------------------------------------------
!    *** PERFORM ALL THE CALCULATION HERE FOR THIS MODEL. ***
!-----------------------------------------------------------------------------------------------------------------------
 !Penman-Monteith Parameters
	 e_sat = exp((16.78*T_a-116.9)/(T_a+237.3)) 
       e_i = rh_a*e_sat
       delta_e = e_sat-e_i 
       DELTA = 4098*e_sat/(T_a+237.3)**2 
       lambda = (2.501-2.361E-3*T_a)*1E6 
       gamma = cp_a*1E-6*P_a/(0.622*lambda*1E-6) 
       d = 2/3*h_c 
       z_om = 0.123*h_c 
       z_oh = 0.1*z_om 
       k = 0.41 
       v_a_in = v_a 
       LAI = 1.5*LOG(h_c*100)-1.4 
       LAI_act = LAI *0.5 
       r_l = 1/(1/r_l_min) 
       r_c = r_l/LAI_act 
       F1 = LOG(max(1E-10,(z_m-d)/z_om))
       F2 = LOG(max(1E-10,(z_h-d)/z_oh))
       F3 = k**2*max(1E-10,v_a_in)
       r_a = (F1*F2)/F3
       F4 = (DELTA*(f*G)+rho_a*cp_a*delta_e/r_a)
       F5 = ((DELTA+gamma*(1+r_c/r_a))*lambda)
       F6 = 3600*area
       
       ET0 = F4/F5*F6
       
       QEVAP = (ET0*lambda)/1000
!-----------------------------------------------------------------------------------------------------------------------
!Set the Outputs from this Model (#,Value)
		Call SetOutputValue(1, ET0) ! Evapotranspiration Rate
		Call SetOutputValue(2, r_a) ! Aerodynamic resistance
		Call SetOutputValue(3, r_c) ! Plant canopy resistance
		Call SetOutputValue(4, QEVAP) ! QEvapotranspiration

!If Needed, Store the Desired Disceret Control Signal Values for this Iteration (#,State)
!Sample Code:  Call SetDesiredDiscreteControlState(1,1)
!-----------------------------------------------------------------------------------------------------------------------
!If Needed, Store the Final value of the Dynamic Variables in the Global Storage Array (#,Value)
!Sample Code:  Call SetValueThisIteration(1,T_FINAL_1)
!-----------------------------------------------------------------------------------------------------------------------
 
      Return
      End 
!-----------------------------------------------------------------------------------------------------------------------

