Structure of subroutine subprograms

Subroutines are similar to functions, yet differ from them in several ways. The most important one is that there is no value associated with the name of the subroutine or, in other words, subroutine has no return statement. Subroutine is invoked using a call statement from anywhere else in your F77 code. Keywords subroutine and end are used to define the beginning and end of a subroutine. Compared to functions which usually have at least one argument, subroutines are often used without any arguments.

The best way to contrast the difference between a function and a subroutine is to rewrite the function for the calculation of the average of three numbers as a subroutine. This might look as follows:

      subroutine AVERAGE(a,b,c,avg) 
      implicit none

      real a, b, c, avg

      avg = (a+b+c)/3.0

      end

The difference between this expression and the function in the previous section is that the subroutine contains one more parameter, avg. This serves as a variable in which the subroutine returns the computed average of the three numbers a, b and c. To call this subroutine, you do not have to declare the type of AVERAGE, as it was needed if it were a function. Simply write:

      :
      CALL average(num1,num2,num3,navg)
      :

The calculated average of the three numbers num1, num2, num3 is then returned in the variable navg.

The most striking difference between functions and subroutines becomes obvious once you have to write a subprogram that returns more than one value. If you still doubt that you will ever need subroutines, check with the following example.

Problem: An interaction between two atoms can be described using the Lennard-Jones potential

$$\Phi(r) = 4\varepsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right] \ ,$$

where $r$ is the separation of two atoms and $\sigma$, $\varepsilon$ are adjustable constants. It was determined that the potential closely reproduces the behavior of argon if $\sigma=3.405\ {\rm\AA}$ and $\varepsilon =0.010323\ {\rm eV}$. Write a subroutine which determines the equilibrium separation of two argon atoms, $r_0$, for which the potential energy is at its absolute minimum. The subroutine should return both the equilibrium separation $r_0$ and the corresponding energy $\Phi(r_0)$.

Solution:

      program CALCMINE
      implicit none

      real sigma, eps, r0, Phi0
      parameter( sigma=3.405, eps=0.010323 )

      call LJ_MIN_E(sigma,eps,r0,Phi0)
      write(*,'("Equilibrium separation = ",F5.3," A")') r0
      write(*,'("Energy at equilibrium  = ",F8.6," eV")') Phi0

      end

!------------------------------------------------------------------

      subroutine LJ_MIN_E(sigma,eps,r0,Phi0)
      implicit none

      real sigma, eps, r0, Phi0

      r0 = sigma*2**(1.0/6.0)
      Phi0 = 4*eps*((sigma/r0)**12-(sigma/r0)**6)

      end

The equilibrium separation $r_0$ corresponds to an extremum of the L-J potential, ${\rm d}\Phi(r)/{\rm d}r=0$. This gives $r_0=\sigma\sqrt[6]{2}$. The corresponding energy is then obtained by substituting $r_0$ in the expression for $\Phi(r)$. The output should look as follows:

Equilibrium separation = 3.822 A
Energy at equilibrium  = -.010323 eV

Last command which is worth mentioning here is return that causes an immediate end of a running subroutine and return to the program which originally called it. This is very convenient once you want to leave a subroutine in the middle and forget the remaining commands. Without using return, you would need to use the goto command for an unconditional jump to the end identifier of your subroutine.