//writen by: E. van der Goot and Giorgio Cifani
//JRC, I-21020 Ispra, Italy
//date: 1996

//based on the original FORTRAN/SQL CGMS and WOFOST code
//produced by Tamme van der Wal, Kees van Diepen and Daniel van Kralingen
//Alterra, P.O. Box 6700 AA, Wageningen, The Netherlands

void PotEvTransp(SWSData WS, SMetData MetData, long globrad, float angot, float daylength, SPenman *pPenmanResults)
{

// The formulae used are according to SUPIT et al,
// System description of the WOFOST 6.0 Crop simulation model
// implemented in CGMS
// The comment about gamma0 at [4.6] is wrong, 0.67 is in mbar/C, not kPa/C

// All calculations are performed in mbar, rather than kPa.

const double Boltzmann = 4.9e-3;
const double gamma = 0.67;                 // psychometric constant, but in mbar/C [4.6]
const float Goudr = 238.102f;                 // from Goudriaan (1977) [4.9]

// NOTE in POTEVT the coeffs are as follows:
// PARAMETER (PSYCON=0.67, REFCFW=0.05, REFCFS=0.15, REFCFC=0.20)
// BUT in documentation for WOFOST.... Canopy is set at .25..

const double ReflCoeffWater = 0.05;             // reflection coeff of water [4.18]
const double ReflCoeffSoil = 0.15;                // reflection coeff of bare soil [4.18]
const double ReflCoeffCanopy = 0.20;          // reflection coeff of canopy [4.18]
const double LatentHeat = 2.45e6;              // Latent heat of vaporization of water

const double Ten2Two = log(2.0/0.033)/log(10.0/0.033);

const double BU = 0.54;

float TG;                                                // Corrected temperature, [4.9]
float TempAvg = (MetData.TempMin + MetData.TempMax) / 2.0f;
double SatVapPress;                             // mbar
double VapPress;                                  // relative sunshine duration
double delta;                                         // slope of saturated vapour pressure/temp curve
double BruntE = 0.1f;                             // Constants used in Brunt
double BruntF = 0.9f;
double Rnl;                                            // net outgoing longwave
double Rnw;                                           // net absorbed
double Rns;
double Rnc;
double Ea;                                             // evaporative demand of atmosphere
double Eac;

float E0;                                                // Penman potential evaporation from a free water surface
float ES0;                                              // Penman potential evaporation from a moist bare soil surface
float ET0;                                               // Penman potential transpiration from a crop canopy

TG = TempAvg + Goudr;
SatVapPress = 6.10588 * exp(17.32491*TempAvg/TG); // mbar [4.9]
VapPress = min(MetData.VapourPressure, (SatVapPress)); // mbar

delta = (Goudr * 17.32491* SatVapPress) / (TG*TG); // mbar/C [4.10]

// We need the relative sunshine duration for the calculation of the net outgoing
// long wave radiation, but check for daylength
if (daylength == 0.0f)
{
RelSunDur = 0.0;
}
else
{
if (MetData.Sunshine != FloatPseudoNull)
{
// We've got it, so use it
RelSunDur = MetData.Sunshine / daylength;
}
else
{
// Estimate the relative sunshine duration from the Angstrom formula, using
// the already calculated value of the global radiation.
// globrad = (long)(angot*(a+b*sun/daylength));
RelSunDur = ((globrad/angot)-WS.AngstromA)/WS.AngstromB;
}
}
 

// Brunt uses two constants (e = 0.1 and f = 0.9)
// In the CGMS source code 22-3-1996, this is made dependent on the latitude value
if (WS.Latitude < 45.0)
{
BruntE = 0.3f;
BruntF = 0.7f;
}

// Calculate net outgoing long wave radiation (J/m2/d according to Brunt)
Rnl = Boltzmann * pow((TempAvg + 273.0), 4) *
(0.56 - 0.08 * sqrt(VapPress)) *
(BruntE + BruntF * RelSunDur); // [4.16]
 

// Net absorbed radiation in mm/d
// BEWARE, globrad is in kJ, and needed in J.
Rnw = ((1.0-ReflCoeffWater)*(globrad*1000.0) - Rnl) / LatentHeat; // [4.18]
Rns = ((1.0-ReflCoeffSoil)*(globrad*1000.0) - Rnl) / LatentHeat; // [4.18]
Rnc = ((1.0-ReflCoeffCanopy)*(globrad*1000.0) - Rnl) / LatentHeat; // [4.18]
 

// Now calculate the evaporative demand of the atmosphere (isothermal evaporation)
// in mm/d. Assume BU constant...
Ea = 0.26 * (SatVapPress - VapPress) * (0.5 + BU * MetData.Wind10*Ten2Two); // [4.19]
Eac = 0.26 * (SatVapPress - VapPress) * (1.0 + BU * MetData.Wind10*Ten2Two);

// Penman formula (1948), but correct for negative Es
E0 = (float)((delta * Rnw + gamma * Ea) / (delta + gamma));
if (E0 < 0.0f) E0 = 0.0f;

ES0 = (float)((delta * Rns + gamma * Ea) / (delta + gamma));
if (ES0 < 0.0f) ES0 = 0.0f;

ET0 = (float)((delta * Rnc + gamma * Eac) / (delta + gamma));
if (ET0 < 0.0f) ET0 = 0.0f;

pPenmanResults ->E0 = E0;
pPenmanResults ->ES0 = ES0;
pPenmanResults ->ET0 = ET0;
}