//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 CSimCro::TotAss(float fDayL, float fAMax, float fEff, float fLAI,
                     float fKDif, float fAvgRad, float fDSinBE,
                     float fSinLD, float fCosLD, float fDiffPP,
                                         float *pfDTGA, float *pfFaPar)
{
float fHour;
float fSinB;          // Solar height
float fPar;           // Photosyntetically active radiation flux
float fParDiff;       // Diffuse part of the photosyntetically active radiation flux at top of the canopy
float fParDir;        // Direct part of the photosyntetically active radiation flux at top of the canopy

// Modifications for FUL, request by Katty Wouters, 8-Jan-1999
float fParDayTop;          // Incoming PAR, top of canopy, daily total
float fParDayAbs;          // Absorbed PAR, daily total
float fVISSLAVG;           // Average absorbtion intensity on sunlit leaf area
float fFaPar = 999.0f;     // Fraction of absorbed photosyntetically active radiation flux

float fRefH;
float fRefS;       // Reflection coefficient of a green leaf canopy
float fKDirBL;     // Extinction coefficient for the direct component of direct light
float fKDirT;      // Extinction coefficient for the total direct radiation flux
float fLAIC;       // Leaf area index at relative distance L in the canopy (L= 0 at top)
float fVISDF;      // Absorbed amount of diffuse radiation flux
float fVIST;       // Absorbed amount of the total direct radiation flux
float fVISD;       // Absorbed amount of direct component of the direct radiation flux
float fVISShD;     // Absorbed amount of the total direct radiation flux by shaded leaves
float fFGRSh;      // Instantaneous gross assimilation rate for shaded leaves
float fVISPP;      // Absorbed amount of the direct radiation flux by leaves perpendicular to the direct beam
float fFGRSun;     // Instantaneous gross assimilation rate for sunlit leaves
float fFSLLA;      // Fraction sunlit leaf area
float fFGL;        // Total instantaneous gross assimilation rate for the whole canopy

float fDTGA = 0.0f;

fParDayTop = 0.0f;
fParDayAbs = 0.0f;

if (fAMax > 0.0f && fLAI > 0.0f)
    {
    for (int nIndex1 = 0; nIndex1 < 3; nIndex1++)
        {
        fHour    = 12.0f + 0.5f * fDayL * m_AXGauss[nIndex1];
        fSinB    = max(0.0f, fSinLD + fCosLD * (float)cos(2.0f * PI * (fHour + 12.0f) / 24.0f)); // Supit 4.24

        // Multiply fAvgRad by 1000 to convert from kJ to J
        fPar     = 0.5f * (fAvgRad * 1000.0f) * fSinB * (1.0f + 0.4f * fSinB) / fDSinBE;
        fParDiff = min(fPar, (fSinB * fDiffPP));
        fParDir  = fPar - fParDiff;

        // --- Was ASSIM code -------------------------------------------------------
        // Extinction coefficients KDIF, KDIRBL, KDIRT
        fRefH   = (1.0f - (float)sqrt(1.0f - m_fSCV)) / (1.0f + (float)sqrt(1.0f - m_fSCV));
        fRefS   = fRefH * 2.0f / (1.0f + 1.6f * fSinB);  // Supit formula 5.12, p.47
        fKDirBL = (0.5f / fSinB) * fKDif / (float)(0.8 * sqrt(1.0f - m_fSCV));
        fKDirT  = fKDirBL * (float)sqrt(1.0f - m_fSCV);

        float fFGross = 0.0f;
        float fParAbs = 0.0f;

        // Three-point Gaussian integration over LAI
        for (int nIndex2 = 0; nIndex2 < 3; nIndex2++)
            {
            fLAIC   = fLAI * m_AXGauss[nIndex2];

            // Absorbed diffuse radiation (VISDF),light from direct
            // origine (VIST) and direct light(VISD)
            fVISDF  = (1.0f - fRefS)  * fParDiff * fKDif   * (float)exp(-fKDif   * fLAIC);
            fVIST   = (1.0f - fRefS)  * fParDir  * fKDirT  * (float)exp(-fKDirT  * fLAIC);
            fVISD   = (1.0f - m_fSCV) * fParDir  * fKDirBL * (float)exp(-fKDirBL * fLAIC);

            // Absorbed flux in W/m2 for shaded leaves and assimilation
            fVISShD = fVISDF + fVIST - fVISD;
            fFGRSh  = fAMax * (1.0f - (float)exp(-fVISShD * fEff / max(2.0f, fAMax)));

            // Direct light absorbed by leaves perpendicular on direct
            // beam and assimilation of sunlit leaf area
            fVISPP  = (1.0f - m_fSCV) * fParDir / fSinB;

            if (fVISPP <= 0.0f)
               fFGRSun = fFGRSh;
            else
               fFGRSun = fAMax * (1.0f - (fAMax - fFGRSh) *
                         (1.0f - (float)exp(-fVISPP * fEff / max(2.0f, fAMax))) / (fEff * fVISPP));

            // Fraction of sunlit leaf area (FSLLA) and local
            // assimilation rate (FGL)
            fFSLLA  = (float)exp(-fKDirBL * fLAIC);
            fFGL    = fFSLLA * fFGRSun + (1.0f - fFSLLA) * fFGRSh;

            // Integration
            fFGross += (fFGL * m_AWGauss[nIndex2]);

            // Average absorbtion intensity on sunlit leaf area
            // equation 14 in Spitters, 1986, part II
            fVISSLAVG = fVISShD + (1.0f - m_fSCV) * fKDirBL * fParDir;

            // PAR radiation absorbed by canopy (Gauss weighted sum of values at 3 depths)
            fParAbs += ((fFSLLA * fVISSLAVG + (1.0f - fFSLLA) * fVISShD)*m_AWGauss[nIndex2]);

            }

        fFGross *= fLAI;
        fParAbs *= fLAI;

        // --- End of ASSIM code ----------------------------------------------------

        fDTGA += (fFGross * m_AWGauss[nIndex1]);

        fParDayTop += (fPar * m_AWGauss[nIndex1]);
        fParDayAbs += (fParAbs * m_AWGauss[nIndex1]);
        }

    fDTGA *= fDayL;
    fFaPar = 100.0f * (fParDayAbs / fParDayTop);
    }

// and the return values
*pfDTGA = fDTGA;
*pfFaPar = fFaPar;
}