fortran 90 で w=curl(u) と書いてみたかった

速度場 u=(u_x(x,y,z),u_y(x,y,z),u_z(x,y,z)) の curl 演算 ω=∇×u の計算を fortran で

      w= curl(u)

と書いてみたかったので試したコード。手許の WinXP + MinGW g95 (デフォルト値コンパイル) で動作した。手許のノートPCで試した配列サイズは倍精度, 128³×3枚である。

      MODULE typedef
        TYPE VEC3Dxyz ! REAL(8), 3次元配列, 3成分へのポインタ
          REAL(8),POINTER,DIMENSION(:,:,:) :: x,y,z
        END TYPE
        TYPE pVEC3Dxyz ! VEC3Dxyz 構造体へのポインタ
          TYPE(VEC3Dxyz),POINTER :: P
        END TYPE

        INTEGER,PARAMETER::NN=128
        REAL(8),DIMENSION(0:NN-1,0:NN-1,0:NN-1),TARGET :: wk1, wk2, wk3
          ! ↑メモリ上の配列の実体をアロケートしているのはこれだけ
      END MODULE

      PROGRAM main
        USE typedef
        TYPE(VEC3Dxyz),TARGET :: u_
        TYPE(pVEC3Dxyz) :: u,w,curl
        u_%x=>wk1; u_%y=>wk2; u_%z=>wk3; u%P=>u_;

        w=curl(u) ! 実は u, w の指し示す実体は同じものなので危険かも

        DO k=0,NN-1; DO j=NN/2,NN/2; DO i=NN/2,NN/2;
          PRINT*,i,j,k,w%P%x(i,j,k),w%P%y(i,j,k),w%P%z(i,j,k)
        END DO; END DO; END DO;
      END PROGRAM

      TYPE(pVEC3Dxyz) FUNCTION curl(u)
        ! 配列にデータを入れているだけ(curl演算ではない)
        ! このコードを実行すると u のデータは書き換えられる。
        USE typedef
        TYPE(pVEC3Dxyz) :: u
        DO k=0,NN-1;
          PRINT*,k
          DO j=0,NN-1;
            DO i=0,NN-1; ! ↓ そのうちに curl 演算に書き換えたい
              u%P%x(i,j,k)=1000000000.+i*1000000+j*1000+k
              u%P%y(i,j,k)=2000000000.+i*1000000+j*1000+k
              u%P%z(i,j,k)=3000000000.+i*1000000+j*1000+k
            END DO;
          END DO;
        END DO;
        curl=u
      END FUNCTION

戻り値を TYPE(VEC3Dxyz) にしたコードも試作してみたが、手元の WinXP + MinGW g95 (デフォルト値コンパイル) で実行すると "Exception: Stack Overflow" を吐き出してコケた。