Проблема с операторами EQUIVALENCE в коде Fortran 77

Я работаю над тем, чтобы заставить работать код трассировки лучей, и я думаю, что, возможно, изолировал проблему. Я новичок в работе с Fortran 77, но хотел бы получить больше опыта использования этого языка (даже если он устарел). У меня есть несколько операторов EQUIVALENCE в одной из подпрограмм, которые можно использовать для передачи переменных в подпрограмму (здесь это может быть половиной проблемы).

Подпрограмма:

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-

Непосредственно перед вызовом подпрограммы или записи ELECTX я поместил несколько операторов печати в подпрограмму/запись RINDEX.
Я проверяю несколько входных данных непосредственно перед вызовом 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

Что я получаю из этого небольшого участка кода:

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

Так что это многословный способ спросить - что происходит? Я ожидал, что такие переменные, как f, например, будут переданы в подпрограмму QPARAB, поэтому, когда я печатаю в подпрограмме, я ожидаю увидеть F = 7. Вероятно, я в корне неправильно понимаю что-то простое. Как я уже упоминал, тот факт, что я не могу получить такие переменные, как F, в подпрограмме QPARAB, на самом деле является большой проблемой, потому что следующие вычисления выходят на 0 или NaN. Я ожидаю, что это будет иметь какую-то ценность. Так может данные как-то не попадают? Все остальное (на данный момент), кажется, в какой-то степени работает.

Откуда этот код:

И я использую небольшой сценарий оболочки (это может быть полный беспорядок):

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

Подпрограммы g77, которые я использую, были загружены по адресу: http://hpc.sourceforge.net/, и, наконец, я получаю та же ошибка при использовании 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)