Skip to content

Commit

Permalink
Applied formatting with miss-hit. Runs on Octave (#147)
Browse files Browse the repository at this point in the history 
* Applied formatting with miss-hit. Runs on Octave

* Update EXE
  • Loading branch information
@BSchilperoort
BSchilperoort authored Mar 1, 2023
1 parent 63fc611 commit 6832a30
Show file tree
Hide file tree
Showing 119 changed files with 8,768 additions and 8,679 deletions.
Binary file modified run_model_on_snellius/exe/STEMMUS_SCOPE
View file
Binary file not shown.
12 changes: 6 additions & 6 deletions src/+equations/Planck.m
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
function Lb = Planck(wl,Tb,em)
function Lb = Planck(wl, Tb, em)

c1 = 1.191066e-22;
c2 = 14388.33;
if nargin<3
if nargin < 3
em = ones(size(Tb));
end
Lb = em.* c1*(wl*1e-9).^(-5)./(exp(c2./(wl*1e-3*Tb))-1);
end

Lb = em .* c1 * (wl * 1e-9).^(-5) ./ (exp(c2 ./ (wl * 1e-3 * Tb)) - 1);

end
24 changes: 12 additions & 12 deletions src/+equations/Root_properties.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,15 @@
% Subfunction - Root - Properties %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%REFERENCES
function[Rl]=Root_properties(Rl,rroot,frac,BR)
%%% INPUTS
% BR = 10:1:650; %% [gC /m^2 PFT]
% rroot = 0.5*1e-3 ; % 3.3*1e-4 ;%% [0.5-6 *10^-3] [m] root radius
%%% OUTPUTS
root_den = 250*1000; %% [gDM / m^3] Root density Jackson et al., 1997
R_C = 0.488; %% [gC/gDM] Ratio Carbon-Dry Matter in root Jackson et al., 1997
Delta_Rltot = BR/R_C/root_den/(pi*(rroot^2)); %% %% root length index [m root / m^2 PFT]
Delta_Rl = Delta_Rltot*frac;
Rl = Rl + Delta_Rl;
end
%% REFERENCES
function [Rl] = Root_properties(Rl, rroot, frac, BR)
%%% INPUTS
% BR = 10:1:650; %% [gC /m^2 PFT]
% rroot = 0.5*1e-3 ; % 3.3*1e-4 ;%% [0.5-6 *10^-3] [m] root radius
%%% OUTPUTS
root_den = 250 * 1000; %% [gDM / m^3] Root density Jackson et al., 1997
R_C = 0.488; %% [gC/gDM] Ratio Carbon-Dry Matter in root Jackson et al., 1997
Delta_Rltot = BR / R_C / root_den / (pi * (rroot^2)); %% %% root length index [m root / m^2 PFT]
Delta_Rl = Delta_Rltot * frac;
Rl = Rl + Delta_Rl;
end
6 changes: 3 additions & 3 deletions src/+equations/Soil_Inertia0.m
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
function [GAM] = Soil_Inertia0(cs,rhos,lambdas)
% soil thermal inertia
GAM = sqrt(cs*rhos*lambdas); % soil thermal intertia
function [GAM] = Soil_Inertia0(cs, rhos, lambdas)
% soil thermal inertia
GAM = sqrt(cs * rhos * lambdas); % soil thermal intertia
43 changes: 21 additions & 22 deletions src/+equations/Soil_Inertia1.m
View file
Original file line number Diff line number Diff line change
@@ -1,35 +1,34 @@
function [GAM] = Soil_Inertia1(SMC)
global theta_s0
%soil inertia method by Murray and Verhoef (
global theta_s0
% soil inertia method by Murray and Verhoef (

% % parameters
% % parameters

theta_s = theta_s0; %(saturated water content, m3/m3)
Sr = SMC/theta_s;
theta_s = theta_s0; % (saturated water content, m3/m3)
Sr = SMC / theta_s;

%fss = 0.58; %(sand fraction)
gamma_s = 0.27; %(soil texture dependent parameter)
dels = 1.33; %(shape parameter)
% fss = 0.58; %(sand fraction)
gamma_s = 0.27; % (soil texture dependent parameter)
dels = 1.33; % (shape parameter)

ke = exp(gamma_s * (1 - power(Sr, gamma_s - dels)));

ke = exp(gamma_s*(1- power(Sr,(gamma_s - dels))));
phis = theta_s0; % (phis == theta_s)
lambda_d = -0.56 * phis + 0.51;

phis = theta_s0; %(phis == theta_s)
lambda_d = -0.56*phis + 0.51;
QC = 0.20; % (quartz content)
lambda_qc = 7.7; % (thermal conductivity of quartz, constant)

QC = 0.20; %(quartz content)
lambda_qc = 7.7; %(thermal conductivity of quartz, constant)
lambda_s = (lambda_qc^(QC)) * lambda_d^(1 - QC);
lambda_wtr = 0.57; % (thermal conductivity of water, W/m.K, constant)

lambda_s = (lambda_qc^(QC))*lambda_d^(1-QC);
lambda_wtr = 0.57; %(thermal conductivity of water, W/m.K, constant)
lambda_w = (lambda_s^(1 - phis)) * lambda_wtr^(phis);

lambda_w = (lambda_s^(1-phis))*lambda_wtr^(phis);
lambdas = ke * (lambda_w - lambda_d) + lambda_d;

lambdas = ke*(lambda_w - lambda_d) + lambda_d;
Hcs = 2.0 * 10^6;
Hcw = 4.2 * 10^6;

Hcs = 2.0*10^6;
Hcw = 4.2*10^6;
Hc = (Hcw * SMC) + (1 - theta_s) * Hcs;

Hc = (Hcw * SMC)+ (1-theta_s)*Hcs;

GAM = sqrt(lambdas.*Hc);
GAM = sqrt(lambdas .* Hc);
4 changes: 2 additions & 2 deletions src/+equations/calc_msoc_fraction.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,5 +5,5 @@ soil organic carbon (MSOC). Density of organic matter is 1300 kg m-3,
and the particle density of mineral material is 2700 kg m-3.
%}

fraction = MSOC.*2700./((MSOC.*2700)+(1-MSOC).*1300);
end
fraction = MSOC .* 2700 ./ ((MSOC .* 2700) + (1 - MSOC) .* 1300);
end
10 changes: 5 additions & 5 deletions src/+equations/calc_rssrbs.m
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
function [rss,rbs] = calc_rssrbs(SMC,LAI,rbs)
global SaturatedMC ResidualMC fieldMC
aa=3.8;
rss = exp((aa+4.1)-aa*(SMC-ResidualMC(1))/(fieldMC(1)-ResidualMC(1)));
rbs = rbs*LAI/4.3;
function [rss, rbs] = calc_rssrbs(SMC, LAI, rbs)
global SaturatedMC ResidualMC fieldMC
aa = 3.8;
rss = exp((aa + 4.1) - aa * (SMC - ResidualMC(1)) / (fieldMC(1) - ResidualMC(1)));
rbs = rbs * LAI / 4.3;
110 changes: 55 additions & 55 deletions src/+equations/calczenithangle.m
View file
Original file line number Diff line number Diff line change
@@ -1,64 +1,64 @@
function [Fi_s,Fi_gs,Fi_g,Omega_s] = calczenithangle(Doy,t,Omega_g,Fi_gm,Long,Lat)
%
% author: Christiaan van der Tol ([email protected])
% date: Jan 2003
% update: Oct 2008 by Joris Timmermans ([email protected]):
% - corrected equation of time
% Oct 2012 (CvdT) comment: input time is GMT, not local time!
%
% function [Fi_s,Fi_gs,Fi_g]= calczenithangle(Doy,t,Omega_g,Fi_gm,Long,Lat)
%
% calculates pi/2-the angle of the sun with the slope of the surface.
%
% input:
% Doy day of the year
% t time of the day (hours, GMT)
% Omega_g slope azimuth angle (deg)
% Fi_gm slope of the surface (deg)
% Long Longitude (decimal)
% Lat Latitude (decimal)
%
% output:
% Fi_s 'classic' zenith angle: perpendicular to horizontal plane
% Fi_gs solar angle perpendicular to surface slope
% Fi_g projected slope of the surface in the plane through the solar beam and the vertical
%
function [Fi_s, Fi_gs, Fi_g, Omega_s] = calczenithangle(Doy, t, Omega_g, Fi_gm, Long, Lat)
%
% author: Christiaan van der Tol ([email protected])
% date: Jan 2003
% update: Oct 2008 by Joris Timmermans ([email protected]):
% - corrected equation of time
% Oct 2012 (CvdT) comment: input time is GMT, not local time!
%
% function [Fi_s,Fi_gs,Fi_g]= calczenithangle(Doy,t,Omega_g,Fi_gm,Long,Lat)
%
% calculates pi/2-the angle of the sun with the slope of the surface.
%
% input:
% Doy day of the year
% t time of the day (hours, GMT)
% Omega_g slope azimuth angle (deg)
% Fi_gm slope of the surface (deg)
% Long Longitude (decimal)
% Lat Latitude (decimal)
%
% output:
% Fi_s 'classic' zenith angle: perpendicular to horizontal plane
% Fi_gs solar angle perpendicular to surface slope
% Fi_g projected slope of the surface in the plane through the solar beam and the vertical
%

%parameters (if not already supplied)
if nargin<6
Long = 13.75; % longitude
Lat = 45.5; % latitude
if (nargin<4)
Omega_g = 210; % aspect
Fi_gm = 30; % slope angle
% parameters (if not already supplied)
if nargin < 6
Long = 13.75; % longitude
Lat = 45.5; % latitude
if nargin < 4
Omega_g = 210; % aspect
Fi_gm = 30; % slope angle
end
end
end

%convert angles into radials
G = (Doy-1)/365*2*pi; % converts day of year to radials
Omega_g = Omega_g/180*pi; % converts direction of slope to radials
Fi_gm = Fi_gm/180*pi; % converts maximum slope to radials
Lat = Lat/180*pi; % converts latitude to radials
% convert angles into radials
G = (Doy - 1) / 365 * 2 * pi; % converts day of year to radials
Omega_g = Omega_g / 180 * pi; % converts direction of slope to radials
Fi_gm = Fi_gm / 180 * pi; % converts maximum slope to radials
Lat = Lat / 180 * pi; % converts latitude to radials

%computes the declination of the sun
d = 0.006918-0.399912*cos(G )+ 0.070247*sin(G )- ...
0.006758*cos(2*G)+ 0.000907*sin(2*G)- ...
0.002697*cos(3*G)+ 0.00148*sin(3*G);

%Equation of time
Et = 0.017 + .4281 * cos(G) - 7.351 * sin(G) - 3.349 * cos(2*G) - 9.731 * sin(2*G);
% computes the declination of the sun
d = 0.006918 - 0.399912 * cos(G) + 0.070247 * sin(G) - ...
0.006758 * cos(2 * G) + 0.000907 * sin(2 * G) - ...
0.002697 * cos(3 * G) + 0.00148 * sin(3 * G);

%computes the time of the day when the sun reaches its highest angle
tm = 12+(4*(-Long)-Et)/60; % de Pury and Farquhar (1997), Iqbal (1983)
% Equation of time
Et = 0.017 + .4281 * cos(G) - 7.351 * sin(G) - 3.349 * cos(2 * G) - 9.731 * sin(2 * G);

%computes the hour angle of the sun
Omega_s = (t-tm)/12*pi;
% computes the time of the day when the sun reaches its highest angle
tm = 12 + (4 * (-Long) - Et) / 60; % de Pury and Farquhar (1997), Iqbal (1983)

%computes the zenithangle (equation 3.28 in De Bruin)
Fi_s = acos(sin(d)*sin(Lat)+cos(d)*cos(Lat).*cos(Omega_s));
% computes the hour angle of the sun
Omega_s = (t - tm) / 12 * pi;

%computes the slope of the surface Fi_g in the same plane as the solar beam
Fi_g = atan(tan(Fi_gm).*cos(Omega_s-Omega_g));
% computes the zenithangle (equation 3.28 in De Bruin)
Fi_s = acos(sin(d) * sin(Lat) + cos(d) * cos(Lat) .* cos(Omega_s));

%computes the angle of the sun with the vector perpendicular to the surface
Fi_gs = Fi_s + Fi_g;
% computes the slope of the surface Fi_g in the same plane as the solar beam
Fi_g = atan(tan(Fi_gm) .* cos(Omega_s - Omega_g));

% computes the angle of the sun with the vector perpendicular to the surface
Fi_gs = Fi_s + Fi_g;
Loading

0 comments on commit 6832a30

Please sign in to comment.