Abaqus の UMAT ユーザーサブルーチンを implicit none を使ってすっきり書こう

## はじめに
ひょんなとこから Abaqus の UMAT をコーディングをすることになりました。言語は Fortran です。
Fortran なんて触ったこともないのでとりあえずググってみます。
すると Fortran の入門記事がいろいろと見つかるわけですが、どれもこれも最初に

> implicit none を使え

と書いてあるわけです。先人がみな言うのですから、これはぜひとも従いたいところです。

ところが、Abaqus のユーザーサブルーチンのサンプルコードを見ても `implicit none` を使用しているものはほぼありません。いったいこれはどうしたことでしょう？

Abaqus のユーザーサブルーチンでは `implicit none` は使えないのでしょうか？

## 要点
結論から言うと、Abaqus のユーザーサブルーチンで `implicit none` は使えます。使えますが、ひと工夫必要です。
Abaqus のインターフェースやドキュメントは Fortran77 を引きずっていますが、Fortran90 以降のモダンな機能は使えます。
というわけで `implicit none` を使いつつ、module や function や type を使って読みやすいコードを書くとよいです。

## 方法

Abaqus のユーザーサブルーチンで `implicit none` を使うのが難しい理由は2つあります。
1つはユーザーサブルーチンの引数の型が暗黙に決まっていて分からないこと。
1つはインクルードファイル中で宣言されているパラメータで方が宣言されていないこと、です。
`implicit none` を宣言しない場合 Fortran では変数名で型が決まります。変数名が `IJKLMN` で始まっていれば `INTEGER` それ以外は `REAL` です。
が、ここにも落とし穴があります。REAL は単精度浮動小数点数を表す型ですが、Abaqus のような FEM の計算には単精度では足りないような……？

こうした謎を解くカギは、インクルードファイルにあります。ABA_PARAM.INC の中身を見てみましょう。これは Abaqus2018 の ABA_PARAM.INC です。
``` fortran:ABA_PARAM.INC
      implicit real*8(a-h,o-z)
      parameter (nprecd=2)
```
1行目に `implicit real*8(a-h,o-z)` とあります。デフォルトでは `REAL` 型となるところを `REAL*8` 型にせよということです。
 `REAL*8` は倍精度です。Abaqus は倍精度浮動小数点数で計算をしているのがこれでわかります。
というわけで、上記の1つ目の理由はこれで解決です。変数名が `IJKLMN` で始まっていれば `INTEGER` それ以外は `REAL*8` ということになります。[^1]

[^1]: インクルードファイルは ABA_PARAM.INC 以外にもいくつかあります、Abaqus のバージョンによっても内容が変わるかもしれないのでそれぞれの環境で確認しましょう。

2行目の `parameter (nprecd=2)` ですが、これは定数の定義です。 n で始まっているのでこれは `INTEGER` です。
何に使っているのかわからないですが、おそらく Abaqus の組み込みサブルーチンなどで使うのでしょう。
`implicit none` を使うとこの定数の型が指定されていないためにエラーになります。これについては、インクルードファイルをあきらめて、コード中で都度定義することにします。


## 実際のコーディング

ここでは実例として弾性体の UMAT の例を見せます。
Fortran では自由記法と固定記法の二つの書き方がありますが、今回は固定記法を使います。理由はドキュメントやサンプルコードがほぼ固定記法になっているため、そこの変換の手間を避けるためです。
自由記法で不便がないなら自由記法を使う方がいいとは思います。

### UMAT サブルーチン内に実装を書かかない
UMAT サブルーチン自体は Abaqus のドキュメントでインターフェースが提供されています。
これをなるべく触らないようにするため、 UMAT サブルーチン内では実装せず umat_impl サブルーチンを呼び出すだけとします。
umat_impl サブルーチンは関連する function やその他のサブルーチンとまとめるために implements モジュールの中に置くことにします。
そのように実装した UMAT が以下です。
UMAT は `implicit none` にしません。

``` fortran:
      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
      use implements
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
     4 JSTEP(4)
      
      call umat_impl(STRESS,DDSDDE,DSTRAN,NTENS,PROPS,NPROPS)
      RETURN
      END
```

### サブルーチン内でインクルードファイルの定数を定義する

呼び出すモジュール umat_impl はこのように定義します。
UMAT にわたってきた引数のうち、必要なものだけを引数として渡しています。弾性体だと必要なものは少ないです。
それらについては INTEGER だったり REAL*8 だったり REAL*8 の配列だったりするので UMAT のドキュメントやインターフェースと整合するように定義します。

``` fortran:
      module implements
      implicit none
      contains
      subroutine umat_impl(STRESS,DDSDDE,DSTRAN,NTENS,PROPS,NPROPS)
      ! replace aba_param.inc for use implicit none
      integer,parameter :: nprecd=2
      integer :: NTENS,NPROPS
      real*8 :: STRESS(NTENS),DDSDDE(NTENS,NTENS),DSTRAN(NTENS),
     1 PROPS(NPROPS)
!     
!     Process
!     
      end subroutine umat_impl
      end module implements
```

さらに `integer,parameter :: nprecd=2` で ABA_PARAM.INC 内にあった定数を定義します。これについては umat_impl からさらにサブルーチンを呼ぶ場合はそのサブルーチン内でも定義する必要があるようです。

### 完成したコード

以上から、さらに内部の処理を umat_impl に書いたものがこれになります。type や function も使っていますが、これは使えるということを示す例として入れたもので正直使う必要はないです。
このコードは [Tutorial: Write a simple UMAT in ABAQUS](https://simplifiedfem.wordpress.com/about/tutorial-write-a-simple-umat-in-abaqus/) にあるコードをベースとして改変したものです。比較してみると何が変わったのかより分かりやすいです。

``` fortran:elasticUMAT_moduletest.f
      module implements
      implicit none
      type LAMBDA_struct
        real*8 :: A,B,C
      end type LAMBDA_struct
      contains
      function calc_LAMBDA(E,ANU)
      type(LAMBDA_struct)::calc_LAMBDA
      real*8::E,ANU
      integer,parameter :: ONE=1.0D0, TWO=2.0D0
      calc_LAMBDA%A=E/(ONE+ANU)/(ONE-TWO*ANU)
      calc_LAMBDA%B=(ONE-ANU)
      calc_LAMBDA%C=(ONE-TWO*ANU)
      end function calc_LAMBDA
      subroutine umat_impl(STRESS,DDSDDE,DSTRAN,NTENS,PROPS,NPROPS)
      ! replace aba_param.inc for use implicit none
      integer,parameter :: nprecd=2
      integer :: NTENS,NPROPS
      real*8 :: STRESS(NTENS),DDSDDE(NTENS,NTENS),DSTRAN(NTENS),
     1 PROPS(NPROPS)
      integer :: I,J
      type(LAMBDA_struct) :: LAMBDA
      real*8 :: E,ANU
C   ****************************************
C   *** ELASTIC ELEMENT
C   ****************************************
      E=PROPS(1)
      ANU=PROPS(2)
      LAMBDA=calc_LAMBDA(E,ANU)
C   ************************
      DDSDDE=0.0D0
C   ************************
      DDSDDE(1,1)=(LAMBDA%A*LAMBDA%B)
      DDSDDE(2,2)=(LAMBDA%A*LAMBDA%B)
      DDSDDE(3,3)=(LAMBDA%A*LAMBDA%B)
      DDSDDE(4,4)=(LAMBDA%A*LAMBDA%B)
      DDSDDE(5,5)=(LAMBDA%A*LAMBDA%B)
      DDSDDE(6,6)=(LAMBDA%A*LAMBDA%B)
      DDSDDE(1,2)=(LAMBDA%A*ANU)
      DDSDDE(1,3)=(LAMBDA%A*ANU)
      DDSDDE(2,1)=(LAMBDA%A*ANU)
      DDSDDE(2,3)=(LAMBDA%A*ANU)
      DDSDDE(3,1)=(LAMBDA%A*ANU)
      DDSDDE(3,2)=(LAMBDA%A*ANU)
C   ************************
      DO I=1,NTENS
        DO J=1,NTENS
            STRESS(I)=STRESS(I)+DDSDDE(I,J)*DSTRAN(J)
        ENDDO
      ENDDO
        
      end subroutine umat_impl
      end module implements

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
      use implements
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
     4 JSTEP(4)
      
      call umat_impl(STRESS,DDSDDE,DSTRAN,NTENS,PROPS,NPROPS)
      RETURN
      END
```

## 最後に

Abaqus のユーザーサブルーチン界隈は、いまだに固定形式が当たり前だったりとプログラミングの手法がアップデートされていない印象があります。
固定形式はまだ我慢できますが、`implicit none` を使わないのは少しでもやりたい処理が複雑になると変数名がわかりにくくなる上に、変数名にかかわるバグで悩まされることになりいいことはありません。
Abaqus は implicit な変数定義を強要してくる非道なプロダクトですが、この記事を生かして `implicit none` を使いましょう。
ついでに module や function や subroutine や type を使ってコードにより意味を持たせ、処理を分割し、読みやすいコードを書くことで、ユーザーサブルーチンの実装を楽してもらえればと思います。
コードが複雑になればなるほど、コーディングの時間に占めるバグ取りの時間、つまりはコードを読み返す時間は増えます。多少書くのが面倒でも、読みやすさを優先してコーディングしたいものです。

## 参考文献
[Tutorial: Write a simple UMAT in ABAQUS](https://simplifiedfem.wordpress.com/about/tutorial-write-a-simple-umat-in-abaqus/)
[编写ABAQUS子程序时启用implicit none](https://zhuanlan.zhihu.com/p/557929822)