<commit_before>
/*****************************************************************************
 * File		: initialize.c							
 * Function	: Initialization of elemental attributes using relational	
 *		  database							
 ****************************************************************************/

#include "pihm.h"

void Initialize (pihm_struct pihm, N_Vector CV_Y)
{
#ifdef _NOAH_
    int             i;
#endif

    if (verbose_mode)
    {
        printf ("\n\nInitialize data structure\n");
    }

    pihm->elem = (elem_struct *)malloc (pihm->numele * sizeof (elem_struct));
    pihm->riv = (river_struct *)malloc (pihm->numriv * sizeof (river_struct));

    InitMeshStruct (pihm->elem, pihm->numele, pihm->meshtbl);

    InitTopo (pihm->elem, pihm->numele, pihm->meshtbl);

#ifdef _NOAH_
    /* Calculate average elevation of model domain */
    pihm->elevation = AvgElev (pihm->elem, pihm->numele);
#endif

    InitSoil (pihm->elem, pihm->numele, pihm->atttbl, pihm->soiltbl,
#ifdef _NOAH_
        pihm->noahtbl,
#endif
        pihm->cal);

    InitLC (pihm->elem, pihm->numele, pihm->atttbl, pihm->lctbl,
        pihm->cal);

    InitForcing (pihm->elem, pihm->numele, pihm->riv, pihm->numriv,
        pihm->atttbl, pihm->rivtbl, &pihm->forc, pihm->cal);

    InitRiver (pihm->riv, pihm->numriv, pihm->elem, pihm->rivtbl,
        pihm->shptbl, pihm->matltbl, pihm->meshtbl, pihm->cal);

    if (debug_mode)
    {
        CorrectElevation (pihm->elem, pihm->numele, pihm->riv, pihm->numriv);
    }

    InitSurfL (pihm->elem, pihm->numele, pihm->riv, pihm->meshtbl);

#ifdef _NOAH_
    InitLsm (pihm->elem, pihm->numele, pihm->ctrl, pihm->noahtbl, pihm->cal);
#endif

    if (pihm->ctrl.init_type == RELAX)
    {
#ifdef _NOAH_
        ApplyForcing (&pihm->forc, pihm->ctrl.starttime);
#endif
        SaturationIC (pihm->elem, pihm->numele, pihm->riv, pihm->numriv);
    }
    else if (pihm->ctrl.init_type == RST_FILE)
    {
        ReadInit (pihm->filename.init, pihm->elem, pihm->numele, pihm->riv,
            pihm->numriv);
    }

    InitVar (pihm->elem, pihm->numele, pihm->riv, pihm->numriv, CV_Y);

    CalcModelStep (&pihm->ctrl);

//#ifdef _DAILY_
//    InitDailyStruct (pihm);
//#endif
}

void InitMeshStruct (elem_struct *elem, int numele, meshtbl_struct meshtbl)
{
    int             i, j;

    for (i = 0; i < numele; i++)
    {
        elem[i].ind = i + 1;

        for (j = 0; j < 3; j++)
        {
            elem[i].node[j] = meshtbl.node[i][j];
            elem[i].nabr[j] = meshtbl.nabr[i][j];
        }
    }
}

void InitTopo (elem_struct *elem, int numele, meshtbl_struct meshtbl)
{
    int             i, j;
    double          x[3];
    double          y[3];
    double          zmin[3];
    double          zmax[3];

    for (i = 0; i < numele; i++)
    {
        for (j = 0; j < 3; j++)
        {
            x[j] = meshtbl.x[elem[i].node[j] - 1];
            y[j] = meshtbl.y[elem[i].node[j] - 1];
            zmin[j] = meshtbl.zmin[elem[i].node[j] - 1];
            zmax[j] = meshtbl.zmax[elem[i].node[j] - 1];
        }

        elem[i].topo.area =
            0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (y[1] - y[0]) * (x[2] -
                x[0]));
        /* Calculate centroid of triangle */
        elem[i].topo.x = (x[0] + x[1] + x[2]) / 3.0;
        elem[i].topo.y = (y[0] + y[1] + y[2]) / 3.0;

        elem[i].topo.zmin = (zmin[0] + zmin[1] + zmin[2]) / 3.0;
        elem[i].topo.zmax = (zmax[0] + zmax[1] + zmax[2]) / 3.0;
        elem[i].topo.edge[0] =
            sqrt (pow ((x[1] - x[2]), 2) + pow ((y[1] - y[2]), 2));
        elem[i].topo.edge[1] =
            sqrt (pow ((x[2] - x[0]), 2) + pow ((y[2] - y[0]), 2));
        elem[i].topo.edge[2] =
            sqrt (pow ((x[0] - x[1]), 2) + pow ((y[0] - y[1]), 2));
    }

#ifdef _NOAH_
    CalcSlopeAspect (elem, numele, meshtbl);
#endif
}

void InitSoil (elem_struct *elem, int numele, atttbl_struct atttbl,
    soiltbl_struct soiltbl,
#ifdef _NOAH_
    noahtbl_struct noahtbl,
#endif
    calib_struct cal)
{
    int             i;
    int             soil_ind;

    for (i = 0; i < numele; i++)
    {
        soil_ind = atttbl.soil[i] - 1;

        elem[i].soil.depth = elem[i].topo.zmax - elem[i].topo.zmin;
        elem[i].soil.ksath = cal.ksath * soiltbl.ksath[soil_ind];
        elem[i].soil.ksatv = cal.ksatv * soiltbl.ksatv[soil_ind];
        elem[i].soil.kinfv = cal.kinfv * soiltbl.kinfv[soil_ind];
        elem[i].soil.porosity =
            cal.porosity * (soiltbl.smcmax[soil_ind] -
            soiltbl.smcmin[soil_ind]);
        if (elem[i].soil.porosity > 1.0 || elem[i].soil.porosity <= 0.0)
        {
            printf ("Warning: Porosity value out of bounds for element %d",
                i + 1);
            PihmExit (1);
        }
        elem[i].soil.dinf = cal.dinf * soiltbl.dinf[soil_ind];
        elem[i].soil.alpha = cal.alpha * soiltbl.alpha[soil_ind];
        elem[i].soil.beta = cal.beta * soiltbl.beta[soil_ind];

        elem[i].soil.smcmin = soiltbl.smcmin[soil_ind];
        elem[i].soil.smcmax = elem[i].soil.smcmin + elem[i].soil.porosity;

        /* Calculate field capacity and wilting point following Chan and Dudhia 2001 MWR, but replacing Campbell with van Genuchten */
        elem[i].soil.smcwlt =
            WiltingPoint (soiltbl.smcmax[soil_ind],
            soiltbl.smcmin[soil_ind], soiltbl.alpha[soil_ind],
            soiltbl.beta[soil_ind]);
#ifdef _NOAH_
        elem[i].soil.smcwlt =
            (elem[i].soil.smcwlt - elem[i].soil.smcmin) * cal.thetaw +
            elem[i].soil.smcmin;
#endif
        elem[i].soil.smcref =
            FieldCapacity (soiltbl.alpha[soil_ind], soiltbl.beta[soil_ind],
            soiltbl.ksatv[soil_ind], soiltbl.smcmax[soil_ind],
            soiltbl.smcmin[soil_ind]);
#ifdef _NOAH_
        elem[i].soil.smcref =
            (elem[i].soil.smcref - elem[i].soil.smcmin) * cal.thetaref +
            elem[i].soil.smcmin;
#endif
<commit_msg>Apply porosity calibration to smcmin as well<commit_after>
/*****************************************************************************
 * File		: initialize.c							
 * Function	: Initialization of elemental attributes using relational	
 *		  database							
 ****************************************************************************/

#include "pihm.h"

void Initialize (pihm_struct pihm, N_Vector CV_Y)
{
#ifdef _NOAH_
    int             i;
#endif

    if (verbose_mode)
    {
        printf ("\n\nInitialize data structure\n");
    }

    pihm->elem = (elem_struct *)malloc (pihm->numele * sizeof (elem_struct));
    pihm->riv = (river_struct *)malloc (pihm->numriv * sizeof (river_struct));

    InitMeshStruct (pihm->elem, pihm->numele, pihm->meshtbl);

    InitTopo (pihm->elem, pihm->numele, pihm->meshtbl);

#ifdef _NOAH_
    /* Calculate average elevation of model domain */
    pihm->elevation = AvgElev (pihm->elem, pihm->numele);
#endif

    InitSoil (pihm->elem, pihm->numele, pihm->atttbl, pihm->soiltbl,
#ifdef _NOAH_
        pihm->noahtbl,
#endif
        pihm->cal);

    InitLC (pihm->elem, pihm->numele, pihm->atttbl, pihm->lctbl,
        pihm->cal);

    InitForcing (pihm->elem, pihm->numele, pihm->riv, pihm->numriv,
        pihm->atttbl, pihm->rivtbl, &pihm->forc, pihm->cal);

    InitRiver (pihm->riv, pihm->numriv, pihm->elem, pihm->rivtbl,
        pihm->shptbl, pihm->matltbl, pihm->meshtbl, pihm->cal);

    if (debug_mode)
    {
        CorrectElevation (pihm->elem, pihm->numele, pihm->riv, pihm->numriv);
    }

    InitSurfL (pihm->elem, pihm->numele, pihm->riv, pihm->meshtbl);

#ifdef _NOAH_
    InitLsm (pihm->elem, pihm->numele, pihm->ctrl, pihm->noahtbl, pihm->cal);
#endif

    if (pihm->ctrl.init_type == RELAX)
    {
#ifdef _NOAH_
        ApplyForcing (&pihm->forc, pihm->ctrl.starttime);
#endif
        SaturationIC (pihm->elem, pihm->numele, pihm->riv, pihm->numriv);
    }
    else if (pihm->ctrl.init_type == RST_FILE)
    {
        ReadInit (pihm->filename.init, pihm->elem, pihm->numele, pihm->riv,
            pihm->numriv);
    }

    InitVar (pihm->elem, pihm->numele, pihm->riv, pihm->numriv, CV_Y);

    CalcModelStep (&pihm->ctrl);

//#ifdef _DAILY_
//    InitDailyStruct (pihm);
//#endif
}

void InitMeshStruct (elem_struct *elem, int numele, meshtbl_struct meshtbl)
{
    int             i, j;

    for (i = 0; i < numele; i++)
    {
        elem[i].ind = i + 1;

        for (j = 0; j < 3; j++)
        {
            elem[i].node[j] = meshtbl.node[i][j];
            elem[i].nabr[j] = meshtbl.nabr[i][j];
        }
    }
}

void InitTopo (elem_struct *elem, int numele, meshtbl_struct meshtbl)
{
    int             i, j;
    double          x[3];
    double          y[3];
    double          zmin[3];
    double          zmax[3];

    for (i = 0; i < numele; i++)
    {
        for (j = 0; j < 3; j++)
        {
            x[j] = meshtbl.x[elem[i].node[j] - 1];
            y[j] = meshtbl.y[elem[i].node[j] - 1];
            zmin[j] = meshtbl.zmin[elem[i].node[j] - 1];
            zmax[j] = meshtbl.zmax[elem[i].node[j] - 1];
        }

        elem[i].topo.area =
            0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (y[1] - y[0]) * (x[2] -
                x[0]));
        /* Calculate centroid of triangle */
        elem[i].topo.x = (x[0] + x[1] + x[2]) / 3.0;
        elem[i].topo.y = (y[0] + y[1] + y[2]) / 3.0;

        elem[i].topo.zmin = (zmin[0] + zmin[1] + zmin[2]) / 3.0;
        elem[i].topo.zmax = (zmax[0] + zmax[1] + zmax[2]) / 3.0;
        elem[i].topo.edge[0] =
            sqrt (pow ((x[1] - x[2]), 2) + pow ((y[1] - y[2]), 2));
        elem[i].topo.edge[1] =
            sqrt (pow ((x[2] - x[0]), 2) + pow ((y[2] - y[0]), 2));
        elem[i].topo.edge[2] =
            sqrt (pow ((x[0] - x[1]), 2) + pow ((y[0] - y[1]), 2));
    }

#ifdef _NOAH_
    CalcSlopeAspect (elem, numele, meshtbl);
#endif
}

void InitSoil (elem_struct *elem, int numele, atttbl_struct atttbl,
    soiltbl_struct soiltbl,
#ifdef _NOAH_
    noahtbl_struct noahtbl,
#endif
    calib_struct cal)
{
    int             i;
    int             soil_ind;

    for (i = 0; i < numele; i++)
    {
        soil_ind = atttbl.soil[i] - 1;

        elem[i].soil.depth = elem[i].topo.zmax - elem[i].topo.zmin;
        elem[i].soil.ksath = cal.ksath * soiltbl.ksath[soil_ind];
        elem[i].soil.ksatv = cal.ksatv * soiltbl.ksatv[soil_ind];
        elem[i].soil.kinfv = cal.kinfv * soiltbl.kinfv[soil_ind];
        elem[i].soil.smcmin = cal.porosity * soiltbl.smcmin[soil_ind];
        elem[i].soil.smcmax = cal.porosity * soiltbl.smcmax[soil_ind]; 
        elem[i].soil.porosity = elem[i].soil.smcmax - elem[i].soil.smcmin;
        if (elem[i].soil.porosity > 1.0 || elem[i].soil.porosity <= 0.0)
        {
            printf ("Warning: Porosity value out of bounds for element %d",
                i + 1);
            PihmExit (1);
        }
        elem[i].soil.dinf = cal.dinf * soiltbl.dinf[soil_ind];
        elem[i].soil.alpha = cal.alpha * soiltbl.alpha[soil_ind];
        elem[i].soil.beta = cal.beta * soiltbl.beta[soil_ind];


        /* Calculate field capacity and wilting point following Chan and Dudhia 2001 MWR, but replacing Campbell with van Genuchten */
        elem[i].soil.smcwlt =
            WiltingPoint (soiltbl.smcmax[soil_ind],
            soiltbl.smcmin[soil_ind], soiltbl.alpha[soil_ind],
            soiltbl.beta[soil_ind]);
#ifdef _NOAH_
        elem[i].soil.smcwlt =
            (elem[i].soil.smcwlt - elem[i].soil.smcmin) * cal.thetaw +
            elem[i].soil.smcmin;
#endif
        elem[i].soil.smcref =
            FieldCapacity (soiltbl.alpha[soil_ind], soiltbl.beta[soil_ind],
            soiltbl.ksatv[soil_ind], soiltbl.smcmax[soil_ind],
            soiltbl.smcmin[soil_ind]);
#ifdef _NOAH_
        elem[i].soil.smcref =
            (elem[i].soil.smcref - elem[i].soil.smcmin) * cal.thetaref +
            elem[i].soil.smcmin;
#endif

{'fraction_non_alphanumeric': 0.19361917623227548, 'fraction_numerical': 0.013166779203241054, 'mean_word_length': 2.2445235487404163, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 147, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '183001', 'n_tokens_mistral': 4965, 'n_tokens_neox': 4796, 'n_words': 856}