HomeProgrammingOctave

logo

gnu5.2-IC-v0.3 Hydrometeorological method application

Application of the hydrometeorological method for calculating flow rates in basins. Adapted to the Spanish road drainage instruction 5.2-IC of 2016 (pdf)

#############################################################
Application of the rational method according to 5.2-IC (2016)
Pere Hernàndez i Casellas
version 0.3
30 April 2019
GNU General Public License 3.0
http://www.gnu.org
#############################################################

Files:

gnu52IC.m > source code file
data.txt > data file
results.txt > results file

The program gnuIC52 is written in the Octave numerical computing language and must be run with the GNU Octave interpreter, which can be downloaded from: https://www.gnu.org/software/octave/download.html or installed directly on GNU/Linux systems with the command:
sudo apt-get install octave

To see the source code and the details of the formulation used, open the gnu52IC.m code with a text editor. The input data must be written in the 'data.txt' file as explained in the file itself. When the program is run, it reads the data file and generates the results file 'results.txt' with a detailed report of the results.

Source code:

######################################################
# Application of the rational method according to 5.2-IC (2016)
# Pere Hernandez i Casellas
# perehc@geocities.com
# Date 20190429
# Version 0.3
# GNU General Public License 3.0
######################################################

######################################################
# Open files
######################################################

fid = fopen ("data.txt");
fid2 = fopen ("results.txt",'w');

######################################################
# Data reading functions. Component 2.
######################################################

function read1lines(fid)
txt = fgetl (fid); line000 = strsplit (txt);
endfunction

function read3lines(fid)
for i=1:3;
txt = fgetl (fid); line000 = strsplit (txt);
endfor;
endfunction

function read7lines(fid)
for i=1:7;
txt = fgetl (fid); line000 = strsplit (txt);
endfor;
endfunction

function nint8=readint8(fid)
txt = fgetl (fid); line000 = strsplit (txt,":");
nint8=(str2double((line000(2))));
endfunction

function double=readdouble(fid)
txt = fgetl (fid); line000 = strsplit (txt,":");
double=str2double(line000(2));
endfunction

function tchar=readchar(fid)
txt = fgetl (fid); line000 = strsplit (txt,":");
tchar=char(line000(2));
endfunction

######################################################
# Data reading functions. Component 1.
# OPTIONAL
######################################################

function tchart=readchart(fid)
txt = fgetl (fid); line000 = strsplit (txt,":");
tchart=char(line000(1));
endfunction

######################################################
# Read data
######################################################

read7lines(fid);

Nom_projecte=readchar(fid);
Data_calcul=readchar(fid);
Nom_conca=readchar(fid);
Nom_OD=readchar(fid);
Area_conca=readdouble(fid);
Lon_llera=readdouble(fid);
Des_llera=readdouble(fid);
Per_retorn=readint8(fid);

read3lines(fid);

Pd=readdouble(fid);
I1_Id=readdouble(fid);
#IDFtc=readdouble(fid);
#IDF24=readdouble(fid);
kb=readdouble(fid);

read3lines(fid);

Poi=readdouble(fid);
betam=readdouble(fid);
delta50=readdouble(fid);
FT=readdouble(fid);
Tip_OD=readchar(fid);


######################################################
# tc = time of concentration (h) (according to 2.2.2.5)
######################################################

pen_m=0.001*Des_llera/Lon_llera;

tc=0.3*realpow(Lon_llera,0.76)*realpow(pen_m,-0.19);

######################################################
# Calculation of I(T,t)=Id*Fint (according to 2.2.2)
######################################################

# KA = precipitation reduction factor for basin area

if (Area_conca < 1)
  KA=1;
else
  KA=1-((log10(Area_conca))/15);
endif

# Id = mean corrected precipitation intensity mm/h

Id=Pd*KA/24;

######################################################
# Fint = max(Fa,Fb) intensity factor (according to 2.2.2.4)
######################################################

# Fa = factor obtained from the torrent index

Fa=realpow(I1_Id,3.5287-2.5287*realpow(tc,0.1));

# Fb = factor obtained from IDF curves

IDFtc=realpow(I1_Id,(realpow(28,0.1)-realpow(tc,0.1))/(realpow(28,0.1)-1));
IDF24=realpow(I1_Id,(realpow(28,0.1)-realpow(24,0.1))/(realpow(28,0.1)-1));

Fb=kb*IDFtc/IDF24;

Fint=max(Fa,Fb);

ITt=Id*Fint;

######################################################
# Po = Runoff threshold (according to 2.2.3.2)
######################################################

betaPM=betam*FT;
betaDT=(betam-delta50)*FT;

if (Tip_OD==' PM')
Po=Poi*betaPM;
else
Po=Poi*betaDT;
endif


######################################################
# Coe = runoff coefficient (according to 2.2.3)
######################################################

if (Pd*KA > Po)
Coe=(((Pd*KA/Po)-1)*((Pd*KA/Po)+23))/realpow((Pd*KA/Po)+11,2);
else
Coe=0.0;
endif

######################################################
# Kt = uniformity coefficient (according to 2.2.5)
######################################################

Kt=1+((realpow(tc,1.25))/(realpow(tc,1.25)+14));

######################################################
# QT = annual maximum flow for T (according to 2.2.1)
######################################################

QT=ITt*Coe*Area_conca*Kt/3.6;


######################################################
# Write variables (verify correct reading)

######################################################

Nom_projecte;
Data_calcul;
Nom_conca;
Nom_OD;
Area_conca;
Lon_llera;
Des_llera;
Per_retorn;

Pd;
I1_Id;
IDFtc;
IDF24;
kb;

Poi;
betam;
delta50;
FT;
Tip_OD;

pen_m;
tc;

KA;
Id;

Fa;
Fb;
Fint;
ITt;

betaPM;
betaDT;
Po;

Coe;
Kt;

QT;

######################################################
# Write results file
######################################################

fprintf(fid2,"######################################################\n");
fprintf(fid2,"gnu5.2-IC \n");
fprintf(fid2,"Application of the rational method according to 5.2-IC (2016)\n");
fprintf(fid2,"######################################################\n");

fprintf(fid2,"\n")
fprintf(fid2,"1. General data:\n")
fprintf(fid2,"Project: %s\n",Nom_projecte);
fprintf(fid2,"Calculation date: %s\n",Data_calcul);
fprintf(fid2,"Basin name: %s\n",Nom_conca);
fprintf(fid2,"Outlet name: %s\n",Nom_OD);

fprintf(fid2,"\n")
fprintf(fid2,"2. Geometric data:\n")
fprintf(fid2,"Basin area %6.3f km2 \n",Area_conca);
fprintf(fid2,"Watercourse length %5.3f km \n",Lon_llera);
fprintf(fid2,"Average watercourse drop %4.0f m \n",Des_llera);
fprintf(fid2,"Return period %5.0f years \n",Per_retorn);
fprintf(fid2,"Average watercourse slope %5.3f \n",pen_m);

fprintf(fid2,"\n")
fprintf(fid2,"3. Time of concentration: \n")
fprintf(fid2,"tc Time of concentration %4.2f hours\n",tc);

fprintf(fid2,"\n")
fprintf(fid2,"4. Corrected precipitation intensity: \n")
fprintf(fid2,"Pd daily precipitation corresponding to T %7.2f mm \n",Pd);
fprintf(fid2,"KA precipitation reduction factor %5.3f \n",KA);
fprintf(fid2,"Id=Pd*KA/24 mean corrected precipitation intensity %5.1f mm \n", Id);

fprintf(fid2,"\n")
fprintf(fid2,"5. Precipitation intensity for T and tc: \n")
fprintf(fid2,"I1/Id torrent index %4.1f \n", I1_Id);
fprintf(fid2,"Fa factor depending on IDF curves %5.2f \n", Fa);
fprintf(fid2,"Fb factor depending on IDF curves %5.2f \n", Fb);
fprintf(fid2,"IDF(T,tc) %5.2f \n", IDFtc);
fprintf(fid2,"IDF(T,24) %5.2f \n", IDF24);
fprintf(fid2,"Fint=max(Fa,Fb) %5.2f \n", Fint);
fprintf(fid2,"I(T,tc)=Id*Fint precipitation intensity %5.1f mm/h \n", ITt);

fprintf(fid2,"\n")
fprintf(fid2,"6. Corrected runoff threshold: \n")
fprintf(fid2,"betam %4.2f \n", betam);
fprintf(fid2,"betaPM %4.2f \n", betaPM);
fprintf(fid2,"betaDT %4.2f \n", betaDT);
fprintf(fid2,"Poi initial runoff threshold %5.2f \n", Poi);
fprintf(fid2,"Po runoff threshold %5.2f \n", Po);

fprintf(fid2,"\n")
fprintf(fid2,"7. Application of the rational formula: \n");
fprintf(fid2,"QT=I(T,tc)*C*A*Kt/3.6 \n");
fprintf(fid2,"C runoff coefficient %5.3f \n", Coe);
fprintf(fid2,"Kt uniformity coefficient  %5.3f \n", Kt);
fprintf(fid2,"QT annual maximum flow for T  %6.2f m3/s \n", QT);


######################################################
# Close files
######################################################

fclose (fid);
fclose (fid2);