Fortan
The first thing that you will need to do is extract your model results from the Selafin file. The Selafin file is in Binary format which makes things a little complex. The structure of a Selafin file is well described in Appendix 3 of the TELEMAC-2D User Manual. I have not programmed in Fortran since 1990, but I thought I would give it a try.
PROGRAM ReadSelaphin
! http://www.opentelemac.org/index.php?option=com_kunena&func=view&catid=20&id=1162&Itemid=62&lang=en
INTEGER IERR
CHARACTER(LEN=80) :: strTitle
CHARACTER(LEN=32) :: strVar
INTEGER I, NBV1, NBV2, NELEM,NPOIN,NDP, TEMP, ERROR
INTEGER IPARAM(10)
INTEGER, allocatable :: IKLE(:,:)
INTEGER, allocatable :: IPOBO(:)
REAL*4, allocatable :: X(:)
REAL*4, allocatable :: Y(:)
REAL*4 TIME
REAL*4 , allocatable :: RES(:)
OPEN(1, FILE="C:\Telemac\bin\Yu_3d_sm_240940_5.slf", FORM="UNFORMATTED",STATUS='OLD',IOSTAT=IERR)
IF(IERR.NE.0) PRINT*,'ERROR ',IERR
OPEN(2, FILE="C:\Telemac\bin\Sta2.txt", FORM="FORMATTED")
READ(1, END=5, ERR=7) strTitle(1:80)
WRITE(2,*) strTitle(1:80)
!1 record containing the two integers NBV(1) and NBV(2) (NBV(1) the number of variables, NBV(2) with the value of 0),
READ(1, END=5, ERR=7) NBV1, NBV2
WRITE(2,*) NBV1, NBV2
! NBV(1) records containing the names and units of each variable (over 32 characters),
DO 20 j=1, NBV1
READ(1, END=5, ERR=7) strVar(1:32)
WRITE(2,*) strVar(1:32)
20 continue
! 1 record containing the integers table IPARAM (10 integers, of which only 4 are currently being used).
READ(1, END=5, ERR=7) IPARAM(1:10)
WRITE(2,*) IPARAM(1:10)
READ(1, END=5, ERR=7) NELEM,NPOIN,NDP, TEMP
WRITE(2,*) NELEM,NPOIN,NDP, TEMP
WRITE(2,*) "=================IKLE================"
!1 record containing table IKLE (integer array of dimension (NDP,NELEM) which is the connectivity table
allocate(IKLE(NDP,NELEM),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IKLE"
GO TO 1000
end if
READ(1, END=5, ERR=7) IKLE(1:NDP,1:NELEM)
WRITE(2,*) IKLE(1:NDP,1:NELEM)
WRITE(2,*) "=================IPOBO================"
!1 record containing table IPOBO (integer array of dimension NPOIN);
allocate(IPOBO(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) IPOBO(1:NPOIN)
WRITE(2,*) IPOBO(1:NPOIN)
WRITE(2,*) "=================X================"
!1 record containing table X (real array of dimension NPOIN containing the abscissas of the points),
allocate(X(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) X(1:NPOIN)
WRITE(2,*) X(1:NPOIN)
WRITE(2,*) "=================Y================"
!1 record containing table Y (real array of dimension NPOIN containing the ordinates of the points),
allocate(Y(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) Y(1:NPOIN)
WRITE(2,*) Y(1:NPOIN)
! FOR EACH TIMESTEP
allocate(RES(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory RES"
GO TO 1000
end if
do
READ(1, END=5, ERR=7) TIME
WRITE(2,*) TIME
! For each record
DO 30 j=1, NBV1
READ(1, END=5, ERR=7) RES(1:NPOIN)
WRITE(2,*) RES(1:NPOIN)
30 continue
end do
5 PRINT*,'END OF fILE 1'
GO TO 1000
7 PRINT*,'ERROR IN FILE 1'
GO TO 1000
999 WRITE(2,*) strTitle(1:80)
1000 CONTINUE
CLOSE(1)
CLOSE(2)
STOP
END
Please excuse any errors with this file - I have run it successfully but have not had the chance to verify the results. This is because I ended up using MATLAB. I found a nice set of scripts for reading and writing Selafin files on the MATLAB website. It has a set of m-files used to read and write the Selafin file. The main ones that I was using were:
- telheadr.m - reads the header info from a Selafin file.
- telstepr.m - reads a specified time step
- telmean.m - calculates the mean from all time-steps.
MATLAB
I needed to generate the mean flow velocities over a specific timespan so I modified telmean.m to take in a start and end timestep. Next I developed my own script to parse the m.RESULT variable and plot the profiles.
The variables for each time step are stored in an (NPOIN * NBV x 1) array- where NPOIN is the number of points and NBV is the number of variables (e.g. elevation, U, V, W for my model runs). This is a very long array and parsing it into the variables is a somewhat complex task.
In the MATAB scripts, the array is accessible through the variable m.RESULT. The RESULT file for a 3D Selafin file is written out by the vertical layers. Each layer has NPOIN/NPLAN points, where NPLAN is the number of vertical layers. The parser that I wrote ends up creating and populating X,Y,Z,U,V,W arrays which are much more easy to access. It then generates a velocity profile along one of my "channels". There is the code:
function m = telplotprofile(INFILE,CHANNEL,NCCHAN,vectorScale,INTERVAL)
%****M* Telemac/telplotprofile.m
% NAME
% telplotprofile.m
%
% PURPOSE
% Reads a processed Telemac file (e.g. mean of set of time steps) and plots
% the velocity profile. Assumes that the channel mesher has been used.
%
% USAGE
% m = telplotprofile(INFILE,CHANNEL,vectorScale,INTERVAL)
%
% INPUTS
% INFILE - the input file
% CHANNEL - the channel number
% vectorScale - scaling factor for the
% INTERVAL - spacing of the velocity profiles (1 = all, 2 = every
% second)
% NCCHAN - number of cross-channel elements
%
% Example Usage:
% inputFile = 'C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5';
% meanFile = 'C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5_mean.slf';
% telmean2(inputFile,41, 100, meanFile);
% telplotprofile(meanFile,2,9,0.02,3);
%
% author: Patrick Grover
% email: patrick.grover@queensu.ca
% release date: 31-Aug-2011
% For testing run telplotprofile('') from the command line
if isempty(INFILE)
INFILE ='C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5_mean.slf';
CHANNEL = 9; % The channel to plot INPUT
vectorScale = 0.02;
INTERVAL = 1;
NCCHAN = 9; % number of cross-channel elements
end
% Read the output file which now has the averaged values
m = telheadr(INFILE);
m.NPOINCHAN = m.NPOIN / (m.NPLAN*NCCHAN); % the number of points along a profile (downstream
m.NPOINLAY = m.NPOIN / (m.NPLAN); % the number of points on a vertical layer
m.X = zeros(m.NPOINCHAN,1);
m.Y = zeros(m.NPOINCHAN,1);
m.Z = zeros(m.NPOINCHAN,m.NPLAN);
m.U = zeros(m.NPOINCHAN,m.NPLAN);
m.V = zeros(m.NPOINCHAN,m.NPLAN);
m.W = zeros(m.NPOINCHAN,m.NPLAN);
m = telstepr(m,1); %since is an average file - read only the first timestep
% Read in the X and Y values
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
m.X(i) =m.XYZ(POINT,1);
m.Y(i) =m.XYZ(POINT,2);
%fprintf('\nX = %d\n',m.X(i));
%fprintf('\nY = %d\n',m.Y(i));
end
% now read in the z values
for layer=1:m.NPLAN
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
%fprintf('\nPOINT = %d\n',POINT);
POINT = POINT + ((layer-1)*m.NPOINLAY);
%fprintf('\nPOINT = %d\n',POINT);
m.Z(i,layer) =m.RESULT(POINT);
%fprintf('\nZ = %d\n',m.Z(i));
end
end
% now read in the U,V,W values
for layer=1:m.NPLAN
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
%fprintf('\nPOINT = %d\n',POINT);
POINT = POINT + ((layer-1)*m.NPOINLAY);
POINTU = POINT + m.NPOIN;
POINTV = POINT + m.NPOIN*2;
POINTW = POINT + m.NPOIN*3;
m.U(i,layer) = m.RESULT(POINTU);
m.V(i,layer) = m.RESULT(POINTV);
m.W(i,layer) = m.RESULT(POINTW);
end
end
umax = max(max(m.U()));
m.U = m.U ./ umax;
figure; hold on;
for i=1 : INTERVAL: m.NPOINCHAN
bFlag = 1; % a simple flag that indicates if it is the first time through
for layer=1:m.NPLAN
X0 = m.X(i);
Z0 = m.Z(i,layer);
X1 = X0 + (m.U(i,layer))*vectorScale;
Z1 = Z0;
VX = [X0 X1];
VY = [Z0 Z1];
if m.U(i,layer) > 0
plot(VX,VY,'k');
else
plot(VX,VY,'r');
end
if bFlag~=1
BackX = [X0OLD0 X0];
BackZ = [Z0OLD0 Z0];
plot(BackX,BackZ,'k');
FrontX = [X0OLD1 X1];
FrontZ = [Z0OLD1 Z1];
plot(FrontX,FrontZ,'k');
end
bFlag = 0;
X0OLD0 = X0;
X0OLD1 = X1;
Z0OLD0 = Z0;
Z0OLD1 = Z1;
end
end
% Finally plot the surface.
bFlag = 1; % a simple flag that indicates if it is the first time through
for i=1: m.NPOINCHAN
X1 = m.X(i);
Z1 = m.Z(i,1);
if bFlag~=1;
BottomX = [X0 X1];
BottomZ = [Z0 Z1];
fprintf('\nZ = %d\n',BottomZ);
plot(BottomX,BottomZ,'k');
end
X0 =X1;
Z0 = Z1;
bFlag = 0;
end
hold off;
This generates the following velocity profile graph:
