[子程序] 求助，子程序中一个自定义函数的问题

subroutine Front
    dimension fixed(160),equat(60),vecrv(160),gload(60),gstif(1830),estif(16,16),&
               iffix(160),nacva(60),locel(16),ndest(16)
    common/Contro/npoin,nelem,nnode,ndofn,ndime,nstre,ntype,ngaus,nprop,nmats,nvfix,nevab,&
                   icase,ncase,itemp,iprob,nprob
    common/Lgdata/coord(80,2),props(10,5),presc(40,2),asdis(160),eload(25,16),strin(3,225),&
                   nofix(40),ifpre(40,2),lnods(25,8),matno(25)
    common/Work/elcod(2,8),shape(8),deriv(2,8),dmatx(3,3),cartd(2,8),dbmat(3,16),bmatx(3,16),&
                 smatx(3,16,9),posgp(3),weigp(3),gpcod(2,9),neror(24)            
    !integer nfunc

    call Openfile


    nfunc(i,j)=(j*j-j)/2+i            !*********************************问题在此行


    mfron=60
    mstif=1830
!
!  interpret fixity data in vector form
!
    ntotv=npoin*ndofn
    do 100 itotv=1,ntotv
    iffix(itotv)=0
100 fixed(itotv)=0.0
    do 110 idofn=1,ndofn
    ngash=nloca+idofn
    iffix(ngash)=ifpre(ivfix,idofn)
110 fixed(ngash)=presc(ivfix,jdofn)
!
!  change the sign of the last appearance
!  of each node
!
    do 140 ipoin=1,npoin
    klast=0
    do 130 ielem=1,nelem
    do 120 inode=1,nnode
    if(lnods(ielem,inode).ne.ipoin) goto 120
    klast=ielem
    nlast=inode
120 continue
130 continue
    if(klast.ne.0) lnods(klast,nlast)=-ipoin
140 continue
!
!  start by initializing everything that
!  matters to zero
!
    do 150 istif=1,mstif
150 gstif(istif)=0.0
    do 160 ifron=1,mfron
    gload(ifron)=0.0
    equat(ifron)=0.0
    vecrv(ifron)=0.0
160 nacva(ifron)=0
!
!  and prepare for disc reading and writing
!  operations
!
    rewind 2
    rewind 4
    rewind 6
    rewind 8
!
!  enter main element assembly-reduction loop
!
    nfron=0
    kelva=0
    do 380 ielem=1,nelem
    kevab=0
    read(2,*) estif
    do 170 inode=1,nnode
    do 170 idofn=1,ndofn
    nposi=(inode-1)*ndofn+idofn
    locno=lnods(ielem,inode)
    if(locno.gt.0) locel(nposi)=(locno-1)*ndofn-idofn
    if(locno.lt.0) locel(nposi)=(locno+1)*ndofn-idofn
170 continue
!
!  start by looking for existing destinations
!
    do 210 ievab=1,nevab
    nikno=iabs(locel(ievab))
    kexis=0
    do 180 ifron=1,nfron
    if(nikno.ne.nacva(ifron)) goto 180
    kevab=kevab+1
    kexis=1
    ndest(kevab)=ifron
180 continue
    if(kexis.ne.0) goto 210
!
!  we now seek new empty places for
!  destination vector
!
    do 190 ifron=1,mfron
    if(nacva(ifron).ne.0) goto 190
    nacva(ifron)=nikno
    kevab=kevab+1
    ndest(kevab)=ifron
190 continue
!
!  the new places may demand an increase
!  in current frontwidth
!
200 if(nedst(kevab).gt.nfron) nfron=ndest(kevab)
210 continue
!
!  assemble element loads
!
    do 240 ievab=1,nevab
    idest=ndest(ievab)
    gload(idest)=gload(idest)+eload(ielem,ievab)
!
!  assemble the element stifnesses
!  but not in resolution
!
    if(icase.gt.1) goto 230
    do 220 jevab=1,ievab
    jdest=ndest(jevab)
    ngash=nfunc(idest,jdest)
    ngish=nfunc(jdest,idest)
    if(jdest.ge.idest) gstif(ngash)=gstif(ngash)+estif(ievab,jevab)
    if(jdest.lt.idest) gstif(ngish)=gstif(ngish)+estif(ievab,jevab)
220 continue
230 continue
240 continue
!
!  re-examine each element node, to
!  enquire which can be eliminated
!
    do 370 ievab=1,nevab
    nikno=-locel(ievab)
    if(nikno.le.0) goto 370
!
!  find positions of variables ready
!  for elimination
!
    do 350 ifron=1,nfron
    if(nacva(ifron).ne.nikno) goto 350
!
!  extract the cofficients of the 
!  new equation for elimination
!
    if(icase.gt.1) goto 260
    do 250 jfron=1,mfron
    if(ifron.lt.jfron) nloca=nfunc(ifron,jfron)
    if(ifron.ge.jfron) nloca=nfunc(jfron,ifron)
    equat(jfron)=gstif(nloca)
250 gstif(nloca)=0.0
260 continue
!
!  and fxtract the corresponding right
!  hand sides
!
    eqrhs=gload(ifron)
    gload(ifron)=0.0
    kelva=kelva+1
!
!  write equations to disc or to tape
!
    if(icase.gt.1) goto 270
    write(4,*) equat,eqrhs,ifron,nikno
    goto 280
270 write(8,*) eqrhs
    read(4,*) wquat,dummy,idumm,nikno
280 continue
!
!  deal with pivot
!
    pivot=equat(ifron)
    equat(ifron)=0.0
!
!  enquire whether present varizable is
!  free or prescribed
!
    if(iffix(nikno).eq.0) goto 300
!
!  deal with a prescribed deflection
!
    do 290 jfron=1,nfron
290 gload(jfron)=gload(jfron)-fixed(nikno)*equat(jfron)
    goto 340
!
!  eliminate a free variable-deal with 
!  the right hand side first
!
300 do 330 jfron=1,nfron
    gload(jfron)=gload(jfron)-equat(jfron)*eqrhs/pivot
!
!  now deal with the coefficients in core
!
    if(icase.gt.1) goto 320
    if(equat(jfron).eq.0.0) goto 330
    nloca=nfunc(0,jfron)
    do 310 lfron=1,jfron
    ngash=lfron+nloca 
310 gstif(ngash)=gstif(ngash)-equat(jfron)*equat(lfron)/pivot
320 continue
330 continue
340 equat(ifron)=pivot
!
!  record the new vacant space, and reduce
!  frontwidth if possible
!
    nacva(ifron)=0
    goto 360
!
!  complete the element loop in the forward
!  elimination
!
350 continue
360 if(nacva(nfron).ne.0) goto 370
    nfron=nfron-1
    if(nfron.gt.0) goto 360
370 continue
380 continue
!
!  enter back-substitution phase, loop 
!  backwards through variables
!
    do 410 ielva=1,kelva
!
!  read a new equation
!
    backspace 4
    read(4,*) equat,eqrhs,ifron,nikno
    backspace 4
    if(icase.eq.1) goto 390
    backspace 8
    read(8,*) eqrhs
    backspace 8
390 continue
!
!  prepare to back-substitute from the
!  current equation
!
    pivot=equat(ifron)
    if(iffix(nikno).eq.1) vecrv(ifron)=fixed(nikno)
    if(iffix(nikno).eq.0) equat(ifron)=0.0
!
!  back-substitute in the current equation
!
    do 400 jfron=1,mfron 
400 eqrhs=eqrhs-vecrv(jfron)*equat(jfron)
!
!  put the final values where they belong
!
    if(iffix(nikno).eq.0) vecrv(ifron)=eqrhs/pivot
    if(iffix(nikno).eq.1) fixed(nikno)=-eqrhs
    asdis(nikno)=vecrv(ifron)
410 continue
    write(12,900)
 900 format(1ho,5x,13hDisplacements)
    if(ndofn.ne.2) goto 430
    if(ndime.ne.1) goto 420
    write(12,905)
 905 format(1ho,5x,4hnode,6x,5hDisp.,7x,8hRotation)
    goto 440
420 write(12,910)
 910 format(1ho,5x,4hNode,5x,7hX-disp.,7x,7hY-disp.)
430 write(12,915)
 915 format(1ho,5x,4hNode,5x,7hX-disp.,7x,7hY-disp.)
440 continue
    do 450 ipoin=1,npoin
    ngash=ipoin*ndofn
    ngish=ngash-ndofn+1
450 write(12,920) ipoin,(asdis(igash),igash=ngish,ngash)
 920 format(i10,3e14.6)
    write(12,925)
 925 format(1ho,5x,9hReactions)
    if(ndofn.ne.2) goto 470
    if(ndime.ne.1) goto 460
    write(12,930)
 930 format(1ho,5x,4hNode,6x,5hForce,8x,6hMomnet)
    goto 480
460 write(12,935)
 935 format(1ho,5x,4hNode,5x,7hX-force,7x,7hY-froce)
    goto 480
470 write(12,940)
 940 format(1ho ,5x,4hNode,6x,5hForce,6x,9hXz-moment,5x,9hYz-moment)
480 continue
    do 510 ipoin=1,npoin
    nloca=(ipoin-1)*ndofn
    do 490 idofn=1,ndofn
    ngush=nloca+idofn
    if(iffix(ngush).gt.0) goto 500
490 continue
    goto 510
500 ngash=nloca+ndofn
    ngish=nloca+1
    write(12,945) ipoin,(fixed(igash),igash=ngish,ngash)
510 continue
 945 format(i10,3e14.6)
!
!  post front-reset all element connection numbers to positive
!  values for subsequent use in stress calculation
!
    do 520 ielem=1,nelem
    do 520 inode=1,nnode
520 lnods(ielem,inode)=iabs(lnods(ielem,inode))
    return

  end