0

I'm a novice in Fortran and typically code in python. I think i'm asking a simple question, but I haven't been able to find the answer anywhere. I would like to assign a subroutine to a variable in Fortran. I've seen how to do this for functions, but subroutines are different from functions since they do not return values but instead have variables with intent(out).

Context:

We are writing code in fortran for the speed and then interfacing to python with f2py. I am writing a numerical integrator for gravitational systems.

Code summary:

There are three levels of subroutines. The first is integrator, which takes initial conditions and returns final positions. integrator has a do-loop that calls integration_step, which updates the positions and velocities after one step. Lastly, there's the potential which gives us the accelerations of the particles at a given time-step. potential could be many different subroutines, potential_a, potential_b, etc.

the first argument for integrator is potential. How can I make this a variable?

Current implementation

This is pseudo code and simplified.

MODULE integration
contains 
  SUBROUTINE integrator(potenname,...,dt,nstep,x,v,N,xf,vf)
    IMPLICIT NONE
    CHARACTER*200,INTENT(IN)::potenname
    INTEGER,INTENT(IN)::NSTEP
    REAL*8,INTENT(IN),DIMENSION(N) :: x,v
    REAL*8,INTENT(IN)::dt
    REAL*8,INTENT(OUT),DIMENSION(N):: xf,vf
    REAL*8,DIMENSION(N):: a ! acceleration
    INTEGER :: i

    DO i=1,nstep
      ! HOW DO I AVOID NESTING THIS IF STATEMENT?
      IF (potenname.EQ."potential_a") then
        CALL integration_step(potential_a,x,a) 
      ELSE IF (potenname.EQ."potential_b") then
        CALL integration_step(potential_b,x,a)
      END IF
      xf=x
      vf=v
    END DO
  END SUBROUTINE integrator

  SUBROUTINE integration_step(potential,dt,x,v,N,xf,vf)
    IMPLICIT NONE
    EXTERNAL :: POTEN 
    ! other declarations 
    CALL POTEN(x,a)
    xf = x+v*dt
    vf = v+a*dt
  END SUBROUTINE

  SUBROUTINE potential_a(x,N,a)
    IMPLICIT NONE
    REAL*8, DIMENSION(N), INTENT(IN) :: x
    REAL*8, DIMENSION(N), INTENT(OUT) :: a
    ! compute the acceleration given the position
  END SUBROUTINE potential_a

  SUBROUTINE potential_b(x,N,a)
    IMPLICIT NONE
    REAL*8, DIMENSION(N), INTENT(IN) :: x
    REAL*8, DIMENSION(N), INTENT(OUT) :: a
    ! compute the acceleration given the position
  END SUBROUTINE potential_b

END MODULE integration 
   

Desired solution

I would LOVE to just assign the subroutine to a variable instead of nesting if-statements in a do-loop.

perhaps in integration, I could just create a variable that stores the subroutine?

  SUBROUTINE integrator(potenname,...,dt,nstep,x,v,N,xf,vf)
    ! same declarations as above plus more parameters
    HowDoIDeclareThis :: poten
    if (potenname.EQ."potential_a") THEN
      poten=potential_a
    ELSE IF (potenname.EQ."potential_b") THEN
      poten=potential_b
    END IF
 
    DO i=1,nstep
      CALL integration_step(poten,...)
    END DO 
  END SUBROUTINE
  • You will be wanting _procedure pointers_. The linked question has a focus on functions, but subroutines follow in the same way. – francescalus May 17 '23 at 09:42

0 Answers0