EcoExtreML / STEMMUS_SCOPE

Integrated code of SCOPE and STEMMUS
GNU General Public License v3.0
14 stars 2 forks source link

Rm unused variables and equations #239

Open SarahAlidoost opened 4 days ago

SarahAlidoost commented 4 days ago

Description

closes #240

:red_circle: exe should be updated by merging the main into this branch.

Checklist

SarahAlidoost commented 4 days ago

@MostafaGomaa93 in this PR, unused variables and equations are commented out, only botmSoilLevel is removed. see Files changed.

@BSchilperoort the exe file is re-genertaed. please let me know if more fixes are needed.

BSchilperoort commented 8 hours ago

I found it difficult to understand the current code.

Why not reduce the initialization to:

function GroundwaterSettings = readGroundwaterSettings()
    % Activate/deactivate Groundwater coupling
    GroundwaterSettings.GroundwaterCoupling = 0; % (value = 0 -> deactivate coupling, or = 1 -> activate coupling); default = 0, update value to = 1 -> through BMI

    % Initialize the variables (head, temperature, air pressure) at the bottom boundary (start of saturated zone)
    GroundwaterSettings.tempBotm = 0; % groundwater temperature (C), received from MODFLOW through BMI

    % Call MODFLOW layers information (number of aquifer layers and their elevations, etc)
    GroundwaterSettings.gw_Dep = 1.0; % elevation of the top surface aquifer layer
    %     to avoid model crashing, assign minimum value of 1 cm
end

Move the soil thickness calculation to a more appropriate place. It does nothing with the groundwater inputs:

    % Calculate soil layers thickness (cumulative layers thickness; e.g. 1, 2, 3, 5, 10, ......., 480, total soil depth)
    soilThick = zeros(NN, 1); % cumulative soil layers thickness
    soilThick(1) = 0;
    DeltZ = ModelSettings.DeltZ;
    DeltZ_R = ModelSettings.DeltZ_R;

    % Load model settings
    ModelSettings = io.getModelSettings();
    NN = ModelSettings.NN; % Number of nodes
    NL = ModelSettings.NL; % Number of layers

    for i = 2:NL
        soilThick(i) = soilThick(i - 1) + DeltZ_R(i - 1);
    end
    soilThick(NN) = ModelSettings.Tot_Depth; % total soil depth

Also, the indxBotmLayer thing is a bit problematic. It is calculated in readGroundwaterSettings.m based on a groundwater depth value, but this is only run once at initialization, even though it's supposed to skip saturated layers:

        % index of bottom layer after neglecting saturated layers (from bottom to top)
        indxBotm = GroundwaterSettings.indxBotmLayer;

So, shouldn't the following code be run every time step...?

    for i = 1:NL
        midThick = (soilThick(i) + soilThick(i + 1)) / 2;
        if gw_Dep >= soilThick(i) && gw_Dep < soilThick(i + 1)
            if gw_Dep < midThick
                indxBotmLayer_R = i;
            elseif gw_Dep >= midThick
                indxBotmLayer_R = i + 1;
            end
            break
        elseif gw_Dep >= soilThick(i + 1)
            continue
        end
    end

If these things are fixed, then the BMI implementation will be simple: only three variables (coupling toggle, temperature of the groundwater, groundwater head at top aquifer) have to be implemented.

SarahAlidoost commented 7 hours ago

@MostafaGomaa93 can you check the calculation of indxBotmLayer and if it should be calculated every time step?

MostafaGomaa93 commented 6 hours ago

I found it difficult to understand the current code.

Why not reduce the initialization to

Move the soil thickness calculation to a more appropriate place. It does nothing with the groundwater inputs:

Hi @BSchilperoort , its fine for me to split this function into two as you suggested if that will increase the readability of the code

MostafaGomaa93 commented 6 hours ago

So, shouldn't the following code be run every time step...?

    for i = 1:NL
        midThick = (soilThick(i) + soilThick(i + 1)) / 2;
        if gw_Dep >= soilThick(i) && gw_Dep < soilThick(i + 1)
            if gw_Dep < midThick
                indxBotmLayer_R = i;
            elseif gw_Dep >= midThick
                indxBotmLayer_R = i + 1;
            end
            break
        elseif gw_Dep >= soilThick(i + 1)
            continue
        end
    end

Yes, gw_Dep has to be calculated everytime step based on the updated GroundwaterSettings.headBotmLayer

MostafaGomaa93 commented 6 hours ago

@MostafaGomaa93 can you check the calculation of indxBotmLayer and if it should be calculated every time step?

Yes, indxBotmLayer has to be updated everytime step of MODFLOW