//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

// Water limited calculation (was WATFD)

void CSoil::Calculate(float fTra, SSoilPhysicalGroup& SPG, float fRainfall,
                      SSiteParameters& Site, BOOL bAirDu,
                      float fRootD, float fRootDMax,
                      float fEvWMax, float fEvSMax)
{
// NOTE: convert units
fRainfall /= 10.0f;

// Local variables without persistence
float fEVW  = 0.0f;
float fEVS  = 0.0f;

float fPERC;
float fLOSS;
float fWE;
float fPERC1;
float fPERC2;
float fWELOW;
float fDW;
float fDWLOW;
float fWDR;
float fSSPRE;
float fRINPRE = 0.0f;
 

// Actual evaporation rates from surface water if surface storage more than
1 cm, ...
if (m_fSS > 1.0f)
    {
    fEVW = fEvWMax;
    }
//... else from soil surface
else
    {
    if (m_fRIN >= 1.0f)
        {
        fEVS = fEvSMax;
        m_fDSLR = 1.0f;
        }
    else
        {
        float fEVSMXT;

        m_fDSLR += 1.0f;
        fEVSMXT  = fEvSMax * (float)(sqrt(m_fDSLR) - sqrt(m_fDSLR - 1.0f));
        fEVS     = min(fEvSMax, fEVSMXT + m_fRIN);
        }
    }

// Preliminary infiltration rate
if (m_fSS <= 0.1f)
    {
    if (Site.nIFunRN == 0)
        // Without surface storage
        fRINPRE = (float)(1.0 - Site.dNotInf) * fRainfall + m_fSS;

    else
        if (Site.nIFunRN == 1)
            fRINPRE = (float)(1.0 - Site.dNotInf * AFGen(m_NINFTB, 10,
fRainfall)) * fRainfall + m_fSS;
    }
else
    {
    // With surface storage, infiltration limited by SOPE
    float fAvail  = m_fSS + (fRainfall - fEVW);
 
    fRINPRE = min(SPG.fSope, fAvail);
    }

// PERCOLATION

// Equilibrium amount of soil moisture in rooted zone
fWE = SPG.fSoilMoistureContentFC * fRootD;

// Percolation from rooted zone to subsoil equals amount of
// excess moisture in rooted zone (not to exceed conductivity)
fPERC1 = Limit(0.0f, SPG.fSope, (m_fW - fWE) - fTra - fEVS);

// Loss of water at the lower end of the maximum root zone
// equilibrium amount of soil moisture below rooted zone
fWELOW = SPG.fSoilMoistureContentFC * (fRootDMax - fRootD);
fLOSS  = Limit(0.0f, SPG.fKSub, (m_fWLOW - fWELOW) + fPERC1);

// For rice water losses are limited to K0/20
if (bAirDu == TRUE)
    fLOSS = min(fLOSS, 0.05f * SPG.fKO);

// Percolation not to exceed uptake capacity of subsoil
fPERC2 = ((fRootDMax - fRootD) * SPG.fSoilMoistureContentSAT - m_fWLOW) +
fLOSS;
fPERC  = min(fPERC1, fPERC2);

// Adjustment of infiltration rate
m_fRIN = min(fRINPRE, (SPG.fSoilMoistureContentSAT - m_fSM) * fRootD + fTra
+ fEVS + fPERC);

// Rates of change in amounts of moisture W and WLOW
fDW    = - fTra - fEVS - fPERC + m_fRIN;
fDWLOW = fPERC - fLOSS;

// This is calculated once in ITASK == 2 (rate only, not integrated).

m_fRelSM = (m_fSM - SPG.fSoilMoistureContentWP) /
           (SPG.fSoilMoistureContentFC - SPG.fSoilMoistureContentWP) *
100.0f;

// Surface storage and runoff
fSSPRE = m_fSS + (fRainfall - fEVW - m_fRIN);
m_fSS  = (float)min(fSSPRE, Site.dMaxSurfaceStorage);

// Amount of water in rooted zone
m_fW += fDW;

// Amount of water in unrooted, lower part of rootable zone
m_fWLOW += fDWLOW;

m_fTraJWL += fTra;

// Calculation of amount of soil moisture in new rootzone
if ((fRootD - m_fRootDOld) > 0.001f)
    {
    // Water added to root zone by root growth, in cm
    fWDR     = m_fWLOW * (fRootD - m_fRootDOld) / (fRootDMax -
m_fRootDOld);
    m_fWLOW -= fWDR;

    // Total water addition to rootzone by root growth
    m_fWDRT += fWDR;

    // Amount of soil moisture in extended root zone
    m_fW    += fWDR;
    }

// Mean soil moisture content in rooted zone
m_fSM = m_fW / fRootD;

// Save rooting depth
m_fRootDOld = fRootD;
}