Hello all,

The NBS/NRC Steam Table book published in 1984 contains a FORTRAN program that will calculate the various properties of steam. My employer has used this FORTRAN program for many years and has asked me to write it in C++.

The most basic quantity is the saturation temperature for a given pressure. It is calculated through an interation loop that converges within 1-8 iterations.

My problem is that when C++ calculates the saturation temperature the result is incorrect. There are no typos or errors in the C++ code. However, a coworker wrote the same code in VB (with corrected syntax for VB) and the correct results are obtained. Neither of us can explain why this is occuring.

Has anyone else ever encountered this type of behavior between languages? If so do you have any reasons why this would occur?

I believe that the FORTRAN code is still copyrighted so I can't post my code.

Any help on this matter would be greatly appreciated.

Thanks
Matt

...There are no typos or errors in the C++By that statement, are you saying that your code compiles with no complaints from the compiler? This is quite different from saying that there are no errors. I mean, if there were no errors, the results would be whatever you expect them to be, right?

Anyhow...

,,,any reasons...

One possibility: In "traditional" Fortran, and, therefore by default in most (all?) current Fortran compilers, all subprogram arguments in are passed by reference. For any of the Fortran subroutines that modify any values of arguments back in the calling program, you must pass your C++ arguments by reference.

Other possibilities: There are probably a million (maybe more). It's pretty hard to list all of them.

Maybe inconsistent data types of program variables.

Maybe different treatment of arrays in subroutines:
---Fortran arrays are 1-based; C++ are zero-based.
---If there are any 2-D arrays, some Fortran programmers use tricks that are not applicable to C++.

Maybe something "clever" in the original programmer's use of Fortran common data that was subtle enough to escape the C++ programmer's attention.

Maybe something else.

Regards,

Dave

By that statement, are you saying that your code compiles with no complaints from the compiler? This is quite different from saying that there are no errors. I mean, if there were no errors, the results would be whatever you expect them to be, right?

There are no compilier complaints. The code executes but converges on the wrong result.

Maybe different treatment of arrays in subroutines:
---Fortran arrays are 1-based; C++ are zero-based.

Yes, I corrected the two for loops to start at zero instead of 1. There are no 2D arrays used.

Maybe something "clever" in the original programmer's use of Fortran common data that was subtle enough to escape the C++ programmer's attention.


No data is comming from a common block (I hate those. I have spent hours trying to find where information is comming from in them.) As far as clever, I am ruling that out because the FORTRAN code is verbatim from the NBS/NRC book. (By the way there are errors in NBS/NRC book that the original programmer corrected and noted but those errors are not in the saturation temperature routines.)

I will keeping working on it.

Dave, the actual FORTAN code is posted here http://orion.math.iastate.edu/burkardt/f_src/steam/steam.f90 Attention moderator(s) If the code is posted here is it OK for me to post my code?

The routine PSAT_EST and TDPSDT are used to calculate TSAT.

If it is OK to post my code I will.

Thanks for your help.

Matt

All,

The problem has been solved. The cause of the problem was the wrong variable typedef being used in the (infamous) pow function. The pow function was using an int typedef when I should have been passing a double typedef.

All is right with the world again.

Thanks alot.

Matt

There are no compilier complaints. The code executes but converges on the wrong result.



Yes, I corrected the two for loops to start at zero instead of 1.
Two loops? There are a ton of them that deal with arrays, and not all of them start with 1. Also note that all array index values must be adjusted, not just the ones involving loop counters!

For example, if you have something like

double precision, parameter, dimension ( 10 ) :: bp = (/ &
0.7478629D+00, -0.3540782D+00, 0.0D+00, 0.0D+00, &
0.007159876D+00, 0.0D+00, -0.003528426D+00, 0.0D+00, &
0.0D+00, 0.0D+00 /)

double precision v(10)

.
.
.
!
! Set V(I) = ( T_REF / T )**(I-1).
!
v(1) = 1.0D+00
do i = 2, 10
v(i) = v(i-1) * t_ref / t
end do
!
! Set B1, B1T, B1TT.
!
b1 = bp(1) + bp(2) * log ( 1.0D+00 / v(2) )
b1t = bp(2) * v(2) / t_ref



What does your C++ look like?

One way to translate: Make the sizes of all C++ arrays one larger than the corresponding Fortran ones. Change the initialization lists to put a "dummy" value (zero) into the zero position of each C++ array. Then all references remain the same. (Including all loop counters, etc.)



So: The loop counters can stay the same, but the subscript values must be changed everywhere.

OK, OK, this wastes some memory since each array has an unused element. My reaction to this: So what! I can't imagine that this will be a problem. Especially when you consider some alternatives

For example you could...


Make the C++ arrays the same size as the Fortran arrays. Then every single reference to an array must have its index value decreased by one.

Change values of loop counters and keep track of what has to be changed inside each loop. Also keep track of whatever must be changed outside of loops.

or you could do...



Something else.

Regards,

Dave

For completeness, here is my code.

All I am calculating is the saturation temperature, TSAT().

// Function to calculate the vapor pressure.
double PS(double T){
double V,W,B,Q,PL,zz;

double A[8]={-7.8889166, 2.5514255, -6.716169, 33.239495, -105.38479, 174.35319, -148.39348, 48.631602};

if(T>314){
PL=6.3573118-8858.843/T+607.56335*pow(T,-0.6);
return 0.1 * exp(PL);
}
V=T/647.25;
W=abs(1-V);
B=0;
for(int z=0;z<8;z++){
zz=z;
B=B+A[z]*pow(W,(zz+2)/2);
}
Q=B/V;
return 22.093*exp(Q);
}

// Function to calculate TDPSDT
double TDPSDT(double T){
double C,V,W,Y,Q,B,zz;//,TDPSDT;

double A[8]={-7.8889166, 2.5514255, -6.716169, 33.239495, -105.38479, 174.35319, -148.39348, 48.631602};

V=T/647.25;
W=1-V;
B=0;
C=0;
for(int z=1;z<9;z++){
zz=z;
Y=A[z-1]*pow(W,(zz+1)/2);
C=C+Y/W*(0.5-0.5*z-1/V);
B=B+Y;
}
Q=B/V;
return 22.093*exp(Q)*C;
}


// Function to calculate TSAT
double TSAT(double P){

double TG,PP,DP,PL;
int k;

if(P > 22.05) return 0;

PL=2.302585+log(P);
TG=372.83+PL*(27.7589+PL*(2.3819+PL*(0.24834+PL*0. 0193855)));
for(k=0;k<8;k++){
if(TG<273.15) TG=273.15;
if(TG>647) TG=647;
PP = PS(TG);
DP = TDPSDT(TG);
if(abs(1-(PP/P)) < .00001) break;
TG=TG*(1+(P-PP)/DP);
}
return TG;
}



int main()
{
double P,tsat;
int pause;

P=2.0;
tsat=TSAT(P);
cout << tsat-273.15 << endl;
cin >> pause;
return 0;
}


This code works correctly and produces the correct results. The problem was in the two seperate pow functions. We added the variable zz of type double and passed it the pow function. Before the variable z was passed incorrectly as type int.

Thanks
Matt

For completeness, here is my code.
...So you did it the "Right Way," and not the simple-minded "easy way" of my post. I like the "Right Way" better for the long run, but in my experience the "easy way" gets me up and running faster so that I can compare results with known good output values. Then I can refine things one function at a time to make it cleaner.

...We added the variable zz of type double and passed it the pow function....

Instead of creating a new "dummy" variable, a nicer way might be to use floating point constants in integer expressions where it makes a difference (not a bad habit to get into, I'm thinking, even if it does not make an obvious difference). For this particular expression, do this with either or both of the integer literals:

for(int z=1;z<9;z++){
Y=A[z-1]*pow(W,(z+1.0)/2.0);
.
.
.



Regards,

Dave

...
This code works correctly and produces the correct results....

I don't know what you are using to decide what you mean by "correct results," so I tried to compare the results of your code with the results from the original Fortran that you referenced.

I mad a minor change to your main() to print more digits.

Here's the main() function that I used with your code:

int main()
{
double p, t;

p = 2.0;
t = TSAT(p);
cout << scientific << setprecision(15);
cout << "t-273.15 = " << t-273.15 << endl;
return 0;
}


Here is a main program that I used with the original Fortran:

! Program calls tsat in the file steam.f90
!
! http://orion.math.iastate.edu/burkardt/f_src/steam/steam.f90
!
program testtstat
double precision p, t
double precision rhol, rhov

p = 2.0
call tsat(p,t,rhol,rhov)
write(*,'("t-273.15 = ", 1pd21.15)') t-273.15

end program testtstat



Here's the output that I got from your code:


t-273.15 = 2.139205071338439e+02


Here's the output that I got from the Fortran:

t-273.15 = 2.124166097105207D+02


If I were the judge, I wouldn't accept that. I mean, with any floating point operations, there might be slight differences in output results, but I would certainly expect more than two significant decimal digits accuracy for problems of this nature. (If you round off these two output values, they differ by two digits in the third significant decimal digit.)

Now, doing a little detective work, I notice what I believe to be an error in your function PS. I know that you said that there are no typographical errors in your code, but here is what I see in PS():

if(T>314){
PL=6.3573118-8858.843/T+607.56335*pow(T,-0.6);
return 0.1 * exp(PL);
}


But I think that these lines correspond to the following lines in the psat_est subroutine of the Fortran code:

if ( t <= 314.0D+00 ) then

p = 0.1D+00 * exp ( 6.3573118D+00 - 8858.843D+00 / t &
+ 607.56335D+00 * t**( -0.6D+00 ) )


If I change your code to

if(T <= 314){
PL=6.3573118-8858.843/T+607.56335*pow(T,-0.6);
return 0.1 * exp(PL);
}

I get


t-273.15 = 2.124171635956101e+02



Now, there may be other problems too.

I mean, this is closer to the Fortran output value. It agrees with the Fortran result to six significant decimal digits, which might be reasonable for single precision arithmetic, but I'm wondering if it could be closer.


For example, in the Fortran tsat subroutine I see

errtol = sqrt ( epsilon ( errtol ) )
.
.
.
do it = 1, it_max
.
.
if ( abs ( delg ) < errtol ) then
return
end if

end do


But in your TSAT() function, I see

for(k=0;k<8;k++){
.
.
if(abs(1-(PP/P)) < .00001) break;
TG=TG*(1+(P-PP)/DP);
}
return TG;

If the goal is to get as close to the original as possible, I would try to use the same error tolerance:


double errtol = sqrt(numeric_limits<double>::epsilon());
.
.
.
if(abs(1-(PP/P)) < errtol) break;
.
.
.


It may very well turn out that it won't make much difference for this particular value of p, but is it possible that it could affect results for other input values? How many values will you use to compare results? Are there any other "shortcuts" in your code that don't attempt to implement the exact functionality of the original?

The point is, if I were going to replace an existing program with a modernized one, I would compare results. Not just for one value of p, but for lots of them. If I found any discrepancies anywhere I would grab them by the throat and wrestle them to the ground. See Footnote.

Regards,

Dave

Footnote:
From a great old ballad. A love song between me and one of my programs:

"I grab you by the throat
And we come crashing down through the window
On the dirt ground below
And we wrestle in the mud and the blood and the beer"

---Blindside
Fell In Love With the Game


Regards,

Dave

The point is, if I were going to replace an existing program with a modernized one, I would compare results. Not just for one value of p, but for lots of them.


Yes, I agree. I am testing hundreds of points. The TSAT(P) function will start to produce appreciable errors around the 2000 psig or 13.7 MPa pressures. This was taken into account in the original FORTRAN program.

As far as the accuracy is concerned, three digits is plenty. The pressures we are supplied with from the end user/purchaser of our boilers and heat exchangers are no greater then +/- 10%. The same goes for the temperatures. The temperatures can be wildly off if the thermocouples are poorly shielded from radiation effects or the controller isn't compensating for the radiation load.

I have a long way to go on this program still. I have to calculate the entropy, enthalpy, of both saturated liquid and vapor, the inverse of TSAT, and the superheat routines.

The code I linked to in the previous post is a much larger version of the code in the NBS/NRC Steam Table book, it just contained TSAT() and the functions needed by it.

Thanks alot for help on this. I am sure by now you have gathered that I am a pretty green programmer, but I am getting better.

Thanks
Matt

Good song. Here is the youtube video. http://www.youtube.com/watch?v=fxkBwpejig4

Matt

...three digits is plenty....

Maybe I didn't make my point.

Let's get hypothetical. Instead of "me," I'll talk about "my boss."

Suppose "my boss" gave me an assignment to write a replacement for a program that the company has been "using for years."

I can guarandamtee that she would insist that results from a number of test cases agree with the output from the previous program to as many digits as the program is capable of producing.

If the old program produced, say 1.99643321123, and the new program produced 1.98210029432, I know that she would reject it, even though an output value of 2.00 would what we actually use in practical cases. I mean, given that double precision arithmetic should be "good for" something like 15 significant decimal digits, the outputs from the two programs don't even agree to three rounded decimal digits. Therefore, I wouldn't even think of proceeding until I could reconcile (correct) the difference.

Now, if the old program produced, say 1.99999999743 and the new one produced 1.99999999812 I still might have to defend it even though in practical terms 2.00 would be "OK." She's funny that way.

It very well may be that differences in the tenth or eleventh rounded decimal digit can be explained, given certain differences in the floating point implementation, and are not due to any mistakes or "shortcuts" in my transmogrification of the old relic of a program into a modern work or art. It would be my job to show that the observed differences are innocuous, and will not lead to problems with other input values.

Regards,

Dave

Dave,

Maybe I didn't make my point.

Yes, I totally agree. I don't have any way to know what the FORTRAN values are after the 3rd decimal place because that is all that there is on the output. I don't have a FORTRAN compilier (yes I know I can get one for free) and the steam table book only uses at most 5 decimal places and on average 3 decimal places are used. For instance, entropy is at 5, internal energy is at 2, enthalpy is at 2 (sometimes 1), density is at 2, and specific volume is at 5 decimal places.

I have the same values to the 3 decimal places as the FORTRAN program.

That is about as good as I can get right now.

In my previous position outstanding accuracy was required in my CFD (Computational Fluid Dynamics) models but where I am now 3 - 5 decimal places is acceptable.

Thanks Matt. (Your boss sounds pretty tough.)

...and the steam table book only uses at most 5 decimal places A lot of old tables were done with single precision arithmetic on 36-bit computers (IBM 7090 for example). Modern (32-bit) computers tend to give less accurate estimates at single precision and more accurate estimates at double precision. My point is not whether three or four or five significant decimal digits is "good enough," or whether output results are exactly the same, rounded off to the nearest gnat's eyebrow, for old and new versions, but to make sure the expectations are reasonable and are agreed upon before the big revelation of the new whizbang. See Footnote.



...Your boss sounds pretty toughWell, I was speaking hypothetically. Consider it a parable. (Not a fable, however.)

Regards,

Dave

Footnote:

"My recipe for satisfaction in life: Expectations Management."
---davekw7x

(You can't fall out of bed if you sleep on the floor.)

Dave,

I spoke with my managers about accuracy on this program and was assured by them that as long as my program reproduced the FORTRAN results to the same accuracy that everything was good.

In short, "we ain't build'n no watch".

All though I prefer the one that I tell my wife, "Good enough for the girls we go out with".

LOL.

Thanks
Matt

I spoke with my managers...

I like to start with the most strict requirements and then back off if necessary and/or appropriate.



Regards,

Dave

"Say, 'No,' then negotiate."
---Helga's Rule

The function that I am programming will calculate the density corresponding to an input pressure.

This FORTRAN function uses the following data from the common block.

COMMON /QQQQ/ Q0, Q5
COMMON /ACONST/ WM,GASCON,TZ,AA,Z,DZ,Y,UREF,SREF

This is the full COMMON block

IMPLICIT REAL*8 (A-H, O-Z)
COMMON /ACONST/ WM,GASCON,TZ,AA,ZB,DZB,YB,UREF,SREF
COMMON /NCONST/ G(40),II(40),JJ(40),NC
COMMON /ELLCON/ G1,G2,GF,B1,B2,B1T,B2T,B1TT,B2TT
COMMON /BCONST/ BP(10),BQ(10)
COMMON /ADDCON/ ATZ(4),ADZ(4),AAT(4),AAD(4)
DATA ATZ /2*64.D1,641.6D0,27.D1/, ADZ /3*.319D0,1.55D0/,
Z AAT/2*2.D4,4.D4,25.D0/, AAD/34.D0,4.D1,3.D1,1.05D3/
DATA WM /18.0152D0/, GASCON /.461522D0/, TZ /647.073D0/,
Z AA /1.D0/, NC /36/
DATA UREF, SREF /-4328.455039D0, 7.6180802D0/
DATA G1, G2, GF /11.D0,44.333333333333D0,3.5D0/
DATA BP /.7478629D0,-.3540782D0,2*0.D0,.7159876D-2,0.D0,
Z -.3528426D-2,3*0.D0/, BQ /1.1278334D0,0.D0,-.5944001D0,
Z -5.010996D0,0.D0,.63684256D0,4*0.D0/
DATA G /-.53062968529023D3,.22744901424408D4,.7877933302068 7D3,
Z -.69830527374994D2,.17863832875422D5,-.39514731563338D5,
Z .33803884280753D5,-.13855050202703D5,-.25637436613260D6,
Z .48212575981415D6,-.34183016969660D6,.12223156417448D6,
Z .11797433655832D7,-.21734810110373D7,.10829952168620D7,
Z -.25441998064049D6,-.31377774947767D7,.52911910757704D7,
Z -.13802577177877D7,-.25109914369001D6,.46561826115608D7,
Z -.72752773275387D7,.41774246148294D6,.1401635824461 4D7,
Z -.31555231392127D7,.47929666384584D7,.4091266478120 9D6,
Z -.13626369388386D7,.69625220862664D6,-.10834900096447D7,
Z -.22722827401688D6,.38365486000660D6,.6883325794433 2D4,
Z .21757245522644D5,-.26627944829770D4,-.70730418082074D5,
Z -.225D0,-1.68D0,.055D0,-93.0D0/
DATA II /4*0,4*1,4*2,4*3,4*4,4*5,4*6,4*8,2*2,0,4,3*2,4/
DATA JJ / 2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2, 3,5,7,
Z 2,3,5,7,1,3*4,0,2,0,0/
END


My question is this;

What information will these variables (shown below) contain?

COMMON /QQQQ/ Q0, Q5
COMMON /ACONST/ WM,GASCON,TZ,AA,Z,DZ,Y,UREF,SREF

Thanks
Matt

...
What information will these variables (shown below) contain?

COMMON /QQQQ/ Q0, Q5
Here's how I would try to find what happens with the variables in the QQQQ block:

From bottom-up:

The common block QQQQ appears in the following:
Subroutine CORR: COMMON /QQQQ/ Q00,Q11

Subroutine DFIND: COMMON /QQQQ/ Q0,Q5

Subroutine QQ: COMMON /QQQQ/ Q,Q5

Subroutine THERM: COMMON /QQQQ/ QP, QDP

Note that the only one of those that changes values of the common block variables is Subroutine QQ. The others use the values.

Now, I look at things from the top down. Start with a top-level sequence of calculations and see where we run into QQQQ.

To find saturation temperature as a function of pressure, I look at steam_nbs_interact.f

The actual calculation is done here:

C
C CALCULATE SATURATION PROPERTIES AS A FUNCTION OF P
C
104 P=X/FP
IF(P .LE. 0D0 .OR. P .GT. PCRIT) GO TO 18
DLL=0.0
DVV= 0.0
CALL TCORR(TSAT,P,DLL,DVV)
T = TSAT
RT = GASCON*T
CALL BB(T)
GOTO 108
...108 shown below...


Now, back in steam_nbs.f we see that subroutine TCORR calls TSAT and CORR

TSAT invokes functions PS and TDPSTD. Neither of these affects the values of the QQQQ variables (I think).

On the other hand, CORR calls DFIND and THERM.

Now, DFIND calls QQ (which affects values of the QQQQ variables). Then THERM uses those values.

Another path in steam_nbs_interact.f that refers to the QQQQ common block is the one to callculate saturation properties as a function of T.


C
C CALCULATE SATURATION PROPERTIES AS A FUNCTION OF T
C
103 IF (T .GE. TCRIT) GO TO 16
DLL = 0.0
DVV = 0.0
CALL PCORR(T, PSAT, DLL, DVV)
C
C
C CALCULATE SATURATED LIQUID VALUES
C
108 D = DLL
DD = D/FD
TT = TTI(T)
CALL QQ(T,D)
Z = BASE(D,T)
PL = FP*(RT*D*Z + Q0)
CALL THERM(D,T)
.
.
.



Subroutine PCORR calls CORR, which affects the values of the QQQQ variables but they aren't used here (I think) since the subsequent explicit call to QQ calculates them and then THERM uses them.

Look for places in steam_nbs_interact where PCORR, TCORR, QQ, DFIND, and THERM are called. (See if there are any other paths that lead to QQ, CORR, DFIND, and THERM.)


COMMON /ACONST/ WM,GASCON,TZ,AA,Z,DZ,Y,UREF,SREF

The variables ACONST are defined in Subroutine BLOCKDATA.

Maybe their names or maybe the context of their usage can give a clue to the meanings.
For example: The Specific Gas Constant
GASCON = 0.461522

(Value derived as indicated in comments in the gascon() function in the f90 program file):

The specific gas constant R is equal R-bar/M.

Where R-bar is the universal gas constant and M is the molar mass:

R-bar = 8.31441 J/(mol degrees Kelvin)
M = 18.0152 g/mol


Etc.


Note that I might have missed some things and I might even have got some things wrong (gasp), but that's the way I would try to get to the bottom of things.


Regards,

Dave

This project is still ongoing. This week we encountered this FORTRAN statement.

EQUIVALENCE

I am told that it allows an array to be called with a different index. At least that is how this program uses it.

Why would you need this functionallity? To me, it seems that it would only add confusion. I am not sure if the FORTRAN NR routines use this statement.

Can anyone explain the reason/benefit for this EQUIVALENCE keyword/function?

Thanks
Matt

...
Can anyone explain the reason/benefit for this EQUIVALENCE...
I think that nowadays, equivalence is frowned upon for a number of reasons. See, for example, FORBIDDEN / ESOTERIC / OBSOLETE STATEMENTS (http://www.ibiblio.org/pub/languages/fortran/ch1-5.html)


However, let's take a trip back in time: Or as they used to say on the radio, "Return with us now to those thrilling days of yesteryear..." (If that doesn't sound familiar, you can google it --- and that's not something you could do way back when.)

Back in the dawn of the modern computer age, when Fortran and Cobol were the king and queen of computing languages (and PL/I was the crown prince who, like Prince Charles, never got to wear the crown, even though it's still hanging around after all of these years):

In Fortran, variables were static, which meant that there were only a fixed number of words of memory that the entire program and all of its data could occupy at a given time. There was no dynamic memory allocation, and mainframe computers, costing multiple hundreds of thousands of 1960 dollars might have only a few tens of thousands of words of core memory. Accessing more memory required external storage. (The Job Control Language that controlled the running of each program might have provisions to notify the sleepy grad-student operator to mount a "scratch tape" on a 9-track tape drive, and the program had to explicitly swap stuff onto and off of the tape.)

Sometimes equivalence statements were used so that an array declared, say, as 1,000 integers could be used as scratch memory in different parts of the program. Typically this was in a common block that could be made visible to whatever subprograms needed some temporary memory for whatever...

Here's a little program that looks like something out of Fortran II. Back in those days, they used roman numerals (like world wars) for Fortran versions. The IBM 026 keypunches didn't have lower case characters, which was OK, since line printers didn't have lower case characters either. In those days, to use Hollerith fields in print statements the programmer had to count the number of characters and hard-code into the FORMAT statement. Do statements had to have statement numbers, and I/O always required a FORMAT statement. Variables starting with I, J, K, L, M, and N were integers; all others were reals. Columns after 72 on the punch cards were sometimes used for sequence numbers (not seen by the compiler, but useful to humans in case you or the sleepy grad student dropped your card deck). See Footnote.


C 00001
C AN ARRAY OF 1000 INTEGERS USED FOR DIFFERENT 00002
C THINGS AT DIFFERENT POINTS IN THE PROGRAM 00003
C 00004
DIMENSION ISCRATCH(1000) 00005
00006
DIMENSION ISCR(2) 00007
DIMENSION FSCR(2) 00008
00009
EQUIVALENCE (ISCR(1),FSCR(1),ISCRATCH(1)) 00010
00011
ISCR(1)=1 00012
ISCR(2)=2 00013
00014
ISUM=0 00015
DO 10 I=1,2 00016
10 ISUM=ISUM+ISCR(I) 00017
WRITE(6,1000)ISUM 00018
1000 FORMAT(1X,6HISUM = ,I10) 00019
00020
FSCR(1)=1.1 00021
FSCR(2)=2.2 00022
00023
XSUM=0.0 00024
DO 20 I=1,2 00025
20 XSUM = XSUM+FSCR(I) 00026
WRITE(6,2000)XSUM 00027
2000 FORMAT(1X,6HXSUM = ,F10.6) 00028

ISUM=0 00015
DO 30 I=1,2 00016
30 ISUM=ISUM+ISCR(I) 00017
WRITE(6,1000)ISUM 00018
00029
STOP 00020
00021
END 00022


If you save that as "something.for" and compile with GNU g77, you will get an output like

ISUM = 3
XSUM = 3.300000
ISUM =2140772762

The second loop for ISUM shows that the memory previously used for ISCR has been overwritten by the FSCR loop.

(Comment out the EQUIVALENCE statement and you will get the same answer after the second ISCR loop as the first, since in that case each of the ISCR and FSCR arrays occupies its own memory space, separate from the other.)




Another use of EQUIVALENCE, would be to declare an array of something (maybe integers) that could be used as a pool of memory that subprograms could dip into to allocate storage for arrays. I won't bore you with the details, but it's kind of like a C programmer having to write his own malloc() that keeps track of what blocks of the array can be given out to requestors. There were no pointer variables, so integer index variables were used to get into the memory pool, and functions were defined to make the programming a little easier, but it still looked kind of messy.

The index value of the beginning of the next available block was typically put into a common block along with the pool array itself.

Nowadays, with Fortran 90 and beyond, people just use dynamic memory allocation if they have to worry about accessing vast amounts of memory efficiently. In fact, Fortran compilers, with the cooperation of modern operating systems, are able to access very large amounts of memory with no special effort from the programmer.

Finally, here's a little program that uses an equivalence statement to print out the bytes in memory of various program variables. It works for me with g77 or gfortran:

!
! Fortran program to see the bytes in memory
!
program showbytes
implicit none

integer ix, i, iac
real fx
character c(4)
equivalence (ix, c(1), fx)

ix = 123456789
write(*, '("Enter an integer: ")')
read *, ix

write(*,'(i10," decimal = ", Z8.8," hex")') ix, ix

write(*,'("Here are the bytes stored in memory:")')
do i=1,4
iac = iachar(c(i))
write(*,'(i10," decimal = ", z2.2, " hex")') iac, iac
end do

write(*, '(/"Enter a real number: ")')
read *, fx

write(*,'("You entered ",1pg12.6)') fx
write(*,'("Here are the bytes stored in memory:")')
do i=1,4
write(*,'(i10," decimal = ", z2.2, " hex")') &
+ iachar(c(i)), iachar(c(i))
end do

stop

end program showbytes


Here's a run from my Intel/AMD x86-based (little-endian architecture) workstation:


Enter an integer:
12345678
12345678 decimal = 00BC614E hex
Here are the bytes stored in memory:
78 decimal = 4E hex
97 decimal = 61 hex
188 decimal = BC hex
0 decimal = 00 hex

Enter a real number:
12.3
You entered 12.3000
Here are the bytes stored in memory:
205 decimal = CD hex
204 decimal = CC hex
68 decimal = 44 hex
65 decimal = 41 hex




Regards,

Dave

Footnote: Back in those days, it was properly spelled FORTRAN (as an acronym for FORmula TRANslating system). It was totally managed by its originators (IBM) before the international standards guys got their hands on it. IBM used all-caps, and so did the rest of the world. Afterwords, it was declared to be a proper noun and not an acronym, so we get the more modern-looking Fortran. On the other hand they didn't see fit to do the same for, COBOL (COmmon Business-Oriented Language. I never figured out why. (One of those things that I wished I had been able to ask Grace Murray Hopper, but I never got the chance.)

Also, computers were made of wood and they ran on kerosene.

Actually, I just made up that last part. As my dad used to say, nostalgia just ain't what it used to be. Oh, well...

(like world wars)

LOL

We are making good progress on this project now.

Once it is done I will very happy.

Good writeup on the history of FORTRAN and EQIVALENCE.

Thanks Dave

Matt

The re-write of this program is now complete.

If anyone is trying to write their own and has any questions, please ask me.

Thanks
Matt