AD0AE
AD0AE

Reputation: 111

Trouble with EQUIVALENCE statements in Fortran 77 Code

I am working on getting a raytracing code working and I think I may have isolated the problem. I am new to working with Fortran 77, but would like to gain more experience using this language (even if it is dated). I have some EQUIVALENCE statements in one of the subroutines that may be used to pipe variables into the subroutine (this could be half the problem right here).

The subroutine:

c   s/r qparxdp
  SUBROUTINE QPARAB                                                 PARA001 
implicit real*8 (a-h, o-z)
character*8 modx, id
C        PLAIN PARABOLIC OR QUASI-PARABOLIC PROFILE                     PARA002 
C        W(104) = 0. FOR A PLAIN PARABOLIC PROFILE                      PARA003 
C               = 1. FOR A QUASI-PARABOLIC PROFILE                      PARA004 
  COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  PARA005 
  COMMON R(6) /WW/ ID(10),W0,W(400)                                 PARA006 

 EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FC,W(101)),(HM,W(102)),       PARA007 
 1 (YM,W(103)),(QUASI,W(104)),(PERT,W(150))                         PARA008 
 data ipass / 0 /

 ENTRY ELECTX                                                      PARA010 
 print*, W(2), W(6), W(101), W(102), W(103), W(104), W(150)
 print*, ' Electx W(6), f ', F, EARTHR, FC, HM, YM, QUASI, PERT, ipass
ipass = ipass + 1
if(ipass.gt.10000) ipass = 2
if(ipass.eq.1) return
  modx(1) = 'qparab'
  HMAX=HM                                                           PARA011 
x = 0.d0
pxpr = 0.d0
pxpth = 0.d0
pxpph = 0.d0
  H=R(1)-EARTHR                                                     PARA013 
if(f.le.0.d0) print*, ' YM', YM
  FCF2=(FC/F)**2                                                    PARA014 
  CONST=1.d0                                                        PARA015 
  IF (QUASI.EQ.1.d0) CONST=(EARTHR+HM-YM)/R(1)                      PARA016 
  Z=(H-HM)/YM*CONST                                                 PARA017 
  X=dMAX1(0.d0,FCF2*(1.d0-Z*Z))                                     PARA018 
  print*, 'X in qparx', X, Z
  IF (X.EQ.0.d0) GO TO 50                                           PARA019 
  IF (QUASI.EQ.1.d0) CONST=(EARTHR+HM)*(EARTHR+HM-YM)/R(1)**2       PARA020 
  PXPR=-2.d0*Z*FCF2/YM*CONST                                        PARA021 
  50 IF (PERT.NE.0.d0) CALL ELECT1                                     PARA022 
  RETURN                                                            PARA023 
     END                                                            PARA024-

Immediately before the subroutine or entry ELECTX is called I placed some print statements in the RINDEX Subroutine/Entry.
I check a few of the inputs immediately before the call of RINDEX

      ENTRY RINDEX
  write(*,*), 'Starting Rindex in ahnwfnc', F
if(ray.eq.0.d0.and.ipass.eq.0) print*, '  no magnetic field'
ipass = 1
  OM=PIT2*1.d6*F
  C2=C*C
  K2=KR*KR+KTH*KTH+KPH*KPH
  OM2=OM*OM
  VR =C/OM*KR
  VTH=C/OM*KTH
  VPH=C/OM*KPH
  write(*,*), OM, C2, K2, OM2, VR, VTH, VPH, F
  CALL ELECTX

What I get out of this little section of code is:

fstep,fbeg,fend   1.  7.  8.
fbeg,fstep,f   7.  1.  0.
f  7.  7.
f before Rindex  7.
Starting Rindex in ahnwfnc  7.
43982297.2  8.98755431E+10  1.  1.93444246E+15  0.00640514066  0.00231408417
0.000282636641  7.
0.  0.  0.  0.  0.  0.  0.
Electx W(6),f   0.  0.  0.  0.  0.  0.  0. 1

So this is a longwinded way of asking - what is going on? I expected the variables like f, for example, to be passed into the subroutine QPARAB, so when I print in the subroutine, I'd expect to see F = 7. I am probably fundamentally misunderstanding something simple. As I have mentioned, the fact that I can't seem to get variables like F into the subroutine QPARAB is actually a big issue because the following calculations come out to 0s or NaNs. I would expect it to have some value. So maybe the data isn't getting in somehow? Everything else (at this point) seems to be working, to some degree.

Where this code comes from:

And I am using a small shell script (this could be a total mess):

g77 -c -O3  raytr_dp.for readw_dp.for trace_dp.for reach_dp.for backup_d.for dummy.for  \
polcar_d.for printr_d.for rkam_dp.for hamltn_d.for ahwfnc_d.for \
qparxdp.for dipoly_d.for spoints.for ggm_dp.for secnds.for

g77 -o main -O3 raytr_dp.o readw_dp.o trace_dp.o reach_dp.o backup_d.o dummy.o \
polcar_d.o printr_d.o rkam_dp.o hamltn_d.o ahwfnc_d.o \
qparxdp.o dipoly_d.o spoints.o ggm_dp.o secnds.o

The g77 routines I am using were downloaded at: http://hpc.sourceforge.net/ and finally I get the same error using gfortran,

Using built-in specs.
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=/usr/local/gfortran/libexec/gcc/x86_64-apple-darwin13/4.9.0/lto-wrapper
Target: x86_64-apple-darwin13
Thread model: posix
gcc version 4.9.0 (GCC)

Upvotes: 1

Views: 345

Answers (1)

casey
casey

Reputation: 6915

Subroutine QPARAB takes no arguments, e.g. nothing is passed to it. It loads the following variables from common blocks (variables shared between scope) MODX, X, PXPR, PXPTH, PXPPH, PXPT, HMAX, ID, W0, and W. Additionally it declares local scope variables modx and id and then assigns implicit typing to all undeclared variables (which are local in scope).

Your variable of interest, F is equivalent to writing W(6). This says that implicit variable F (type real*8) must have the same memory location as W(6). F isn't passed into this subroutine, it is a name local to the subroutine that is really a specific array member of W. If you want to pass a value into the subroutine into F, you need to set the variable W(6) prior to calling the subroutine. Note that in order to do this you will need W in scope and thus you will need the /WW/ common block referenced in the subroutine you are calling from.

Upvotes: 2

Related Questions