Ошибка Фортрана; forrtl: суровый (408): fort(2)

Я могу нормально собрать свой проект (ОС IA-32, Windows 7. Intel Visual Fortran 11.1.048, интегрированный с Microsoft Visual Studio 2008). Когда я запускаю файл .exe, выдает следующую ошибку:

forrtl: severe (408): fort(2): Subscript #1 of the array dProductIDev has value 5 which is greater than the upper bound of 4

Он также говорит, что ошибка происходит в функции clcMatA при вычислении clcMatA в последних строках (я отметил см. ниже). Код функции:

        function clcMatA(cStress,D,I_dev,dtime,props,ndi,ntens)

        implicit none

        integer :: ndi,ntens
        real*8 :: Z(2,ntens), dProductIDev(ntens,ntens), &
        clcMatA(ntens,ntens),D(ntens,ntens),I_dev(ntens,ntens),&
        cStress(ntens),dProductSigmadev2(ntens,ntens),&
        sigmaDevDyadicProduct(ntens,ntens),identity(ntens,ntens),&
        sigmaDev(ntens),props(5),alpha, beta,dtime,coeff_1,coeff_2,coeff_3


        call identityMatrix(identity,ntens)

        if (normm(cStress,ntens)==0) then

            clcMatA = identity

        else
            alpha = expValVolume/(2*expValDStran)
            beta = (6*expValVolume/(pi*expValPStran))**(1/props(4))*props(3)
            sigmaDev = dev(cStress,I_dev,ntens)
            dProductIDev = matmul(D,I_dev)

            do i=1,ntens
                do j=1,ntens
                    sigmaDevDyadicProduct(i,j)= sigmaDev(j)*sigmaDev(i)
                end do
            end do

            dProductSigmadev2 = matmul(D,sigmaDevDyadicProduct)


            call zVals(Z,sigmaDev,props,ntens)

            do i=1,ntens
                do j=1,ntens
                    clcMatA(i,j) = identity(i,j) + dtime*( (alpha+beta* &
                    normm(sigmaDev,ntens)**(1./props(4)-1.))*dProductIDev(i,j) + & ! The line causing the error
                    beta*(1./props(4)-1.)*normm(sigmaDev,ntens)**(1./props(4)-3.)* &
                    dProductSigmadev2(i,j) )
                end do    
            end do

        end if

        end function

Переменные i, j, expValVolume, expValDStran и expValPStran определены в модуле, в котором содержатся функции clcMatA, dev, identityMatrix и normm.

ntens значение, верхняя граница индексов, передается в функцию со значением 4. Я также проверил его в режиме останова, и он проверяется. Я также заменил ntens на 4!!! при вычислении clcMatA, но я получил ту же ошибку.

Я напечатал параметры, переданные функции, и наткнулся на нечто странное. D должна быть симметричной матрицей и иметь следующие компоненты:

 D:
   174999994368.000        74999996416.0000        74999996416.0000       0.000000000000000E+000
   74999996416.0000        174999994368.000        74999996416.0000       0.000000000000000E+000
   74999996416.0000        74999996416.0000        174999994368.000       0.000000000000000E+000
  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000   49999998976.0000 

 D:
   174999994368.000        74999996416.0000        74999996416.0000       1.853911331209891E-307
   74999996416.0000        174999994368.000        74999996416.0000       2.228101728310706E-312
   74999996416.0000        74999996416.0000        174999994368.000       7.115376174740906E-307
  1.879376978297863E-307  0.000000000000000E+000  0.000000000000000E+000   49999998976.0000 

Поэтому я также изменил D на то, что должно быть, но снова получил ту же ошибку.

(Я должен использовать этот тип переменной, потому что эта функция вызывается коммерческим программным обеспечением конечных элементов (ABAQUS) косвенно, и тип обмениваемых переменных должен совпадать. Я также пробовал использовать real(8), но ничего не изменилось.)

Я не знаю, что там происходит. Есть идеи?