Troubles compiling an old Fortran program (macOS 12 / Intel)

I'm no Fortran programmer, and I'm trying to help a friend compiling some of his old Fortran programs, which he doesn't remember how to do.

I'm trying to get this to work on macOS 12 on an Intel processor using gfortran.

I know there are different revisions of Fortran, and what struck me about the programs he gave me to try and compile is that none of them defines a program symbol at the top, unlike the Fortran hello world examples you can find on the internet.

According to this friend, each and every source code he gave me wans't meant to be a library and could compile and run just fine alone. So I believe these files might be compiled and run as they are.

Researching on this, I found out that maybe, just maybe the fact that no program symbol is needed might be related to this "fixed form'" which I read about here. Which checks with the source code files' extension as well, all ending in .f.

Anyways, I've tried compiling one of the source files using the most plain gfortran command I could think of:

gfortran -o test BL2CaK.f

which yields:

Undefined symbols for architecture x86_64:
  "_main", referenced from:
     implicit entry/start for main executable
ld: symbol(s) not found for architecture x86_64
collect2: error: ld returned 1 exit status

And I understand that the linker cannot resolve the symbol _main. But that's pretty much it.

No other variation of the gfortran command I've tried to run ever produced something different than this, and I have no clue how I should go about solving this. I tried to compile the file against the legacy version of Fortran (F77?):

gfortran -std=legacy -o test BL2CaK.f

to the same result.

Here's the source file for reference (the original source file is indented using tabs as opposed to spaces as shown below, the file just lost formatting while copy-pasting):

c   program: nuclcarfilm1:it assumes nucleation and growth of an adsorbed film of monomers, 
c   whose initial value is given by theta0;
c   X1=initial time in s (it cannot be exactly equal to zero!; if it is zero, set X1=0.0001);
c   X2=final time; NSTEP=number of steps between X1 ed X2 (for a good accuracy, 
c   NSTEP should not be <100);
c   the  number of molecules composing a nucleus is always set equal to 1
c   Er=reversal transmembrane potential in V 
c   E= potenziale transmembrana misurato rispetto a un  resting potential infinitamente negativo, in V
c   z=charge number della gating charge q=ze, con e carica del protone

    COMMON /PATH/ XX(10000),Y(10,10000),DYDX(10,10000)
    DIMENSION VSTART(5),V(5),DV(5),XX(10000),Y(10,10000),DYDX(10,10000)
    COMMON/VALORI/S0,cappaN,cappaR,p
    write(*,*) 'BL2CaK.f'
    write(*,*) 'X2,E,Er,z,a,S0,cappaN,cappaR,smorz,cost'
    read(*,*) X2,E,Er,z,a,S0,cappaN,cappaR,smorz,cost
    write(*,*)'X2=',X2,'E=',E,'Er=',Er,'z=',z,'a=',a,'S0=',S0
    write(*,*)'cappaN=',cappaN,'cappaR=
 ',cappaR,'smorz=',smorz,'cost=',cost
    X1=0.00000
    NSTEP=20000
    Kount=0
c   Fourth-order Runge-Kutta method
    NVAR=5
    VSTART(1)=S0  ! Parte ad un potenziale dove vi sono solo subunits down
    VSTART(2)=0.00000
    VSTART(3)=0.00000
    VSTART(4)=0.00000
    VSTART(5)=0.00000

    DO 11 I=1,NVAR
    V(I)=VSTART(I)
    Y(I,1)=V(I)
11  CONTINUE
    XX(1)=X1
    X=X1
    H=(X2-X1)/NSTEP

    DO 13 K=1,NSTEP
    Kount=Kount+1 
    CALL DERIVS(X,V,DV)
    CALL RK4(V,DV,NVAR,X,H,V,DERIVS)
    X=X+H
    XX(K+1)=X
    DO 12 I=1,NVAR
    Y(I,K+1)=V(I)
    DYDX(I,K+1)=DV(I)
12  CONTINUE
    p=1/(1+exp(-z*E/0.0257)/a) !ze=gating charge;0.0257=kT/e a 25�C
    curr=cost*S0*p*Y(2,K+1)*(1-exp(smorz*(1-Y(2,K+1))))*(Er-E) ! ON current density
    if(K.eq.19900) goto15
    if(Kount.eq.100) goto 14
    goto13
14  write (*,*) X,curr,p,DYDX(2,K+1),Y(2,K+1),Y(1,K+1)/S0
    Kount=0
13  CONTINUE
15  PAUSE
    END
c----------------------------------------------------------
    SUBROUTINE DERIVS(X,Y,DYDX)
    DIMENSION Y(5),DYDX(5)
    COMMON/VALORI/S0,cappaN,cappaR,p 
c   teta=surface coverage by monome
 rs=Y(1);S=surface coverage by channels=Y(2);
c   Sx= "extended" area=Y(3); p=primo integrale=Y(4);
c   q=secondo integrale=Y(5);X=time in s
    DYDX(3)=6.283*cappaR*p*Y(1)*Y(4)  ! time derivative of Sx
    DYDX(2)=(1-Y(2))*DYDX(3) ! time derivative of S
    Y(1)=S0*(1-Y(2))
    DYDX(4)=cappaR*p*Y(1)*Y(5) ! time derivative of p 
    DYDX(5)=cappaN*p*Y(1)
    RETURN
    END
c----------------------------------------------------------
    SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS)
    PARAMETER (NMAX=10)
    DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX)
    HH=H*0.5
    H6=H/6.
    XH=X+HH
    DO 11 I=1,N
    YT(I)=Y(I)+HH*DYDX(I)
11  CONTINUE
    CALL DERIVS(XH,YT,DYT)
    DO 12 I=1,N
    YT(I)=Y(I)+HH*DYT(I)
12  CONTINUE
    CALL DERIVS(XH,YT,DYM)
    DO 13 I=1,N
    YT(I)=Y(I)+H*DYM(I)        
    DYM(I)=DYT(I)+DYM(I)
13  CONTINUE
    CALL DERIVS(X+H,YT,DYT)
    DO 14 I=1,N
    YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I))
14  CONTINUE
    RETURN
    END

It'd be wonderful if someone could advise. Thanks!