Combining two Fortran subroutines for use with ABAQUS

Question

I have two small user-defined subroutines that I have implemented in the finite element software ABAQUS (they define two special types (JTYPE=1,2) of multi-point constraint (MPC) in my finite element model). These subroutines are written in FORTRAN 66/77 as required and are compiled at runtime by ABAQUS using the Intel FORTRAN Compiler. I have verified that they compile and work correctly.

However, I can only use one of these subroutines at a time in my model. This is because they must both use the following name and arguments (in order for ABAQUS to know when to call them and how to use them):

      SUBROUTINE MPC(UE,A,JDOF,MDOF,N,JTYPE,X,U,UINIT,MAXDOF,
     * LMPC,KSTEP,KINC,TIME,NT,NF,TEMP,FIELD,LTRAN,TRAN)

It turns out, however, that I need to be able to simultaneously use both types in a given analysis. Since JTYPE, which is specified for each of the many MPC instances used in my finite element model, is either 1 or 2, the ideal solution is to combine my currently separate subroutines, and switch between the two types inside this new subroutine. I figured that I could do this using an IF(JTYPE .EQ. 1) THEN type syntax.

The problem is that the required blocks of code differ for each JTYPE, even though they are both SUBROUTINE MPC. For my first subroutine (JTYPE=1), the manual requires the following interface:

      SUBROUTINE MPC(UE,A,JDOF,MDOF,N,JTYPE,X,U,UINIT,MAXDOF,
     * LMPC,KSTEP,KINC,TIME,NT,NF,TEMP,FIELD,LTRAN,TRAN)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION A(N),JDOF(N),X(6,N),U(MAXDOF,N),UINIT(MAXDOF,N),
     * TIME(2),TEMP(NT,N),FIELD(NF,NT,N),LTRAN(N),TRAN(3,3,N)


      user coding to define A, JDOF, and, optionally, LMPC


      RETURN
      END

For the second subroutine (JTYPE=2), the manual requires the following interface:

      SUBROUTINE MPC(UE,A,JDOF,MDOF,N,JTYPE,X,U,UINIT,MAXDOF,
     * LMPC,KSTEP,KINC,TIME,NT,NF,TEMP,FIELD,LTRAN,TRAN)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION UE(MDOF),A(MDOF,MDOF,N),JDOF(MDOF,N),X(6,N),
     * U(MAXDOF,N),UINIT(MAXDOF,N),TIME(2),TEMP(NT,N),
     * FIELD(NF,NT,N),LTRAN(N),TRAN(3,3,N)


      user coding to define JDOF, A and, optionally, LMPC


      RETURN
      END

The important difference is that the dimensions of the A and JDOF arrays are different for each type.

I had hoped to get around this by declaring these variables inside the IF loop, as follows:

      SUBROUTINE MPC(UE,A,JDOF,MDOF,N,JTYPE,X,U,UINIT,MAXDOF,
     * LMPC,KSTEP,KINC,TIME,NT,NF,TEMP,FIELD,LTRAN,TRAN)
C
C     CUSTOM MULTI POINT CONSTRAINT
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION A(N),JDOF(N),X(6,N),U(MAXDOF,N),UINIT(MAXDOF,N),
     * TIME(2),TEMP(NT,N),FIELD(NF,NT,N),LTRAN(N),TRAN(3,3,N),
     * A_2(MDOF,MDOF,N),JDOF_2(MDOF,N),UE(MDOF)
C
      IF (JTYPE .EQ. 1) THEN
      DIMENSION A(N),JDOF(N),X(6,N),U(MAXDOF,N),UINIT(MAXDOF,N),
     * TIME(2),TEMP(NT,N),FIELD(NF,NT,N),LTRAN(N),TRAN(3,3,N),
     * A_2(MDOF,MDOF,N),JDOF_2(MDOF,N),UE(MDOF)
        A(1) = 1.
        A(2) = -1.
        JDOF(1) = 3
        JDOF(2) = 3
      END IF
C
      IF (JTYPE .EQ. 2) THEN
      DIMENSION UE(MDOF),A(MDOF,MDOF,N),JDOF(MDOF,N),X(6,N),
     * U(MAXDOF,N),UINIT(MAXDOF,N),TIME(2),TEMP(NT,N),
     * FIELD(NF,NT,N),LTRAN(N),TRAN(3,3,N)
        DATA H /2.000000e-03/
        A(1,1,1) = 1.
        A(2,2,1) = 1.
        A(3,3,1) = 1.
        A(1,1,2) = -1.
        A(1,5,2) = H
        A(2,2,2) = -1.
        A(2,4,2) = H
        A(3,3,2) = -1.
        JDOF(1,1) = 1
        JDOF(2,1) = 2
        JDOF(3,1) = 3
        JDOF(1,2) = 1
        JDOF(2,2) = 2
        JDOF(3,2) = 3
        JDOF(4,2) = 4
        JDOF(5,2) = 5
      END IF
C
      RETURN
      END

However, this causes compilation to fail at runtime, with ABAQUS spitting out the following error:

Error in job Job-1: Problem during compilation - C:\Users\kbodjo\Documents\Abaqus Analyses\custom_MPC_3.for Job Job-1 aborted due to errors.

How can I achieve my goal of combining these two MPC subroutines, given the challenges described above?

EDIT This is how the code was fixed, and it now works:

      SUBROUTINE MPC(UE,A,JDOF,MDOF,N,JTYPE,X,U,UINIT,MAXDOF,
     * LMPC,KSTEP,KINC,TIME,NT,NF,TEMP,FIELD,LTRAN,TRAN)
C
C     CUSTOM MULTI POINT CONSTRAINT
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION UE(MDOF),A(MDOF*MDOF*N),JDOF(MDOF*N),X(6,N),U(MAXDOF,N),
     * UINIT(MAXDOF,N),
     * TIME(2),TEMP(NT,N),FIELD(NF,NT,N),LTRAN(N),TRAN(3,3,N)
C
      IF (JTYPE .EQ. 1) THEN
        A(1) = 1.
        A(2) = -1.
        JDOF(1) = 3
        JDOF(2) = 3
      END IF
C
      IF (JTYPE .EQ. 2) THEN
C
        DATA H /2.000000e-03/
        A((1-1)*MDOF*MDOF + (1-1)*MDOF + 1) = 1.
        A((1-1)*MDOF*MDOF + (2-1)*MDOF + 2) = 1.
        A((1-1)*MDOF*MDOF + (3-1)*MDOF + 3) = 1.
        A((1-1)*MDOF*MDOF + (1-1)*MDOF + 2) = -1.
        A((2-1)*MDOF*MDOF + (5-1)*MDOF + 1) = H
        A((2-1)*MDOF*MDOF + (2-1)*MDOF + 2) = -1.
        A((2-1)*MDOF*MDOF + (4-1)*MDOF + 2) = H
        A((2-1)*MDOF*MDOF + (3-1)*MDOF + 3) = -1.
C
        JDOF((1-1)*MDOF + 1) = 1
        JDOF((1-1)*MDOF + 2) = 2
        JDOF((1-1)*MDOF + 3) = 3
        JDOF((2-1)*MDOF + 1) = 1
        JDOF((2-1)*MDOF + 2) = 2
        JDOF((2-1)*MDOF + 3) = 3
        JDOF((2-1)*MDOF + 4) = 4
        JDOF((2-1)*MDOF + 5) = 5
      END IF
C
      RETURN
      END

Answer 1

The compiling error seems normal to me. As far as I know, you first declare variables before the executable instructions. Now, since you are using fortran 66/77 that doesn't check type of arguments at compilation, one solution is to consider A and JDOF as one 1D array in both cases. Using the fact that fortran stores arrays columnwise, you are all set to get the position of each element by index translation. And if your code is really that simple, it will be very simple too. For an array argument, fortran actually pass the address of the first element.

for example, declare A and JDOF as

DIMENSION A(MDOF*MDOF*N),JDOF(MDOF*MDOF*N)

in both cases the case 1 is simple, and for the case 2, you have these indexes translations:

A(1,1,1) becomes A(1)
A(2,1,1) becomes A(2)

and so on. As a formula for larger cases, for your array that is initially

A(MDOF,MDOF,N)
A(i,j,k) becomes A( (k-1)*MDOF*MDOF + (j-1)*MDOF + i)