在這里附上子程序的內容，麻煩相關人士指正錯誤：

 subroutine vusdfld(

c Read only -

     *   nblock, nstatev, nfieldv, nprops, ndir, nshr, 

     *   jElemUid, kIntPt, kLayer, kSecPt, 

     *   stepTime, totalTime, dt, cmname, 

     *   coordMp, direct, T, charLength, props, 

     *   stateOld, 

c Write only -

     *   stateNew, field )

c

      include 'vaba_param.inc'

c

      dimension props(nprops),

     *          jElemUid(nblock), coordMp(nblock, *), 

     *          direct(nblock, 3, 3), T(nblock,3,3), 

     *          charLength(nblock),

     *          stateOld(nblock, nstatev), 

     *          stateNew(nblock, nstatev),

     *          field(nblock, nfieldv)

      character*80 cmname

      parameter( nrData=6 )

      character*3 cData(maxblk*6)

      dimension jData(maxblk*nrData)

      dimension temp(maxblk*6), eqps(maxblk*6)

   

 

      jStatus = 1

      call vgetvrm( 'TEMP', temp, jData, cData, jStatus )       

      

      jStatus = 1

      call vgetvrm( 'PEEQ', eqps, jData, cData, jStatus )  

       

      

      do k = 1, nblock

    

       if ( stepTime .eq. 0 ) then

           

           

      stateOld(k,1) = 0

      stateOld(k,2) = 0

      stateOld(k,3) = 0

      stateOld(k,4) = 0

      stateOld(k,5) = 20

      stateOld(k,6) = 0

      

       end if

       

        Te = temp(k)

        

        stateNew(k,1) = Te

    

        str = eqps(k)

        

        stateNew(k,2) =  str

        

        str1 = stateOld(k,2)

        

        strra = (str-str1)/dt

        

        stateNew(k,3) =  strra

              

          zpara = strra*exp(300/(Te+273))   

          

          stateNew(k,6) =  zpara 

          

 !       end if

             

           eqc= 0.00158*zpara**0.10662

                              

           eq05=0.018063*zpara**0.08518  

           

           dfull=1.27e6*zpara**(-0.351)

           

      

          if (str.ge. eqc) then

              

              de= (str-eqc)/(eq05-eqc)

              

              Xdrx=1-exp(-0.693*(de)**1.92)                   

              

              if (Xdrx .ge. 0.95) then

                  

              stateNew(k,4) = 1

              

              stateNew(k,5) = dfull

              

              else

              

              stateNew(k,4) =  Xdrx 

 

              stateNew(k,5) = dfull*Xdrx+(1-Xdrx)*20 

              

              endif                

          end if

          

      end do 

      

      return

      end