!##################################################################################################
program psspap
!plane stress/strain problem analysis program
!##################################################################################################
dimension lnd(500,3),x(200),y(200),jz(50,3),pj(50,3),p(500),tl(20),ak(500,100),ake(6,6)

open (6,file='psspapin.txt')
open (8,file='psspapout.txt')
read (6,10) tl
10 format(20a4)
read(6,*) nj,ne,nz,npj,ips
write(8,10) tl

if (ips==1) write(8,20)
if (ips==2) write(8,30)


20 format(/1x,'plane stress problem')
30 format(/1x,'plane strain problem')

npj0=npj
if (npj==0) npj=1
call read_in(nj,ne,nz,nd,npj,npj0,ips,e,pr,t,v,lnd,x,y,jz,pj)

n=2*nj

do i=1,n
 do j=1,nd
  ak(i,j)=0.0
 end do
end do

do ie=1,ne
 call mke(ie,nj,ne,lnd,x,y,e,pr,t,ake)
 do i=1,3
  do ii=1,2
   ih=2*(i-1)+ii
   idh=2*(lnd(ie,i)-1)+ii
   do j=1,3
    do jj=1,2
     l=2*(j-1)+jj
     il=2*(lnd(ie,j)-1)+jj
     idl=il-idh+1
     if (idl>0) then
      ak(idh,idl)=ak(idh,idl)+ake(ih,l)
     end if
    end do
   end do
  end do
 end do
end do

call mf(nj,ne,npj,npj0,n,t,v,lnd,x,y,pj,p)
call rkr(nz,nd,n,jz,ak,p)
call slov(nj,n,nd,ak,p)
call made(ne,nj,n,e,pr,lnd,x,y,p)

close(6)
close(8)

stop
end program psspap

!##################################################################################################
subroutine read_in(nj,ne,nz,nd,npj,npj0,ips,e,pr,t,v,lnd,x,y,jz,pj)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj),jz(nz,3),pj(npj,3)
 read(6,*) e,pr,t,v
 read(6,*) ((lnd(i,j),j=1,3),i=1,ne)
 read(6,*) (x(i),y(i),i=1,nj)
 read(6,*) ((jz(i,j),j=1,3),i=1,nz)
 
 if (npj0/=0) then
  read(6,*) ((pj(i,j),j=1,3),i=1,npj)
 end if
 
 nd=0
 
 do ie=1,ne
  do i=1,3
   do j=1,3
    iw=iabs(lnd(ie,i)-lnd(ie,j))
    if (nd<iw) then
     nd=iw
    end if
   end do
  end do
 end do

 nd=(nd+1)*2
 
 if (ips/=1) then
  e=e/(1.0-pr*pr)
  pr=pr/(1.0-pr)
 end if

end

!##################################################################################################
subroutine mke(ie,nj,ne,lnd,x,y,e,pr,t,ake)
!##################################################################################################
dimension lnd(ne,3),x(nj),y(nj),b(3,6),d(3,3),s(3,6),ake(6,6)
call ma(ie,nj,ne,lnd,x,y,ae)
call md(e,pr,d)
call mb(ie,nj,ne,lnd,x,y,ae,b)

do i=1,3
  do j=1,6
   s(i,j)=0.0
   do k=1,3
    s(i,j)=s(i,j)+d(i,k)*b(k,j)
   end do
  end do
 end do

 do i=1,6
  do j=1,6
   ake(i,j)=0.0
   do k=1,3
    ake(i,j)=ake(i,j)+s(k,i)*b(k,j)*ae*t
   end do
  end do
 end do

end

!##################################################################################################
subroutine ma(ie,nj,ne,lnd,x,y,ae)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj)

 i=lnd(ie,1)
 j=lnd(ie,2)
 k=lnd(ie,3)
 xij=x(j)-x(i)
 yij=y(j)-y(i)
 xik=x(k)-x(i)
 yik=y(k)-y(i)
 ae=0.5*(xij*yik-xik*yij)

return
end

!##################################################################################################
subroutine md(e,pr,d)
!##################################################################################################
dimension d(3,3)

 do i=1,3
  do j=1,3
   d(i,j)=0.0
  end do 
 end do

 d(1,1)=e/(1.0-pr*pr)
 d(1,2)=e*pr/(1.0-pr*pr)
 d(2,1)=d(1,2)
 d(2,2)=d(1,1)
 d(3,3)=0.5*e/(1.0+pr)

return
end

!##################################################################################################
subroutine mb(ie,nj,ne,lnd,x,y,ae,b)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj),b(3,6)

 i=lnd(ie,1)
 j=lnd(ie,2)
 k=lnd(ie,3)

 do ii=1,3
  do jj=1,6
   b(ii,jj)=0.0
  end do
 end do

 b(1,1)=y(j)-y(k)
 b(1,3)=y(k)-y(i)
 b(1,5)=y(i)-y(j)
 b(2,2)=x(k)-x(j)
 b(2,4)=x(i)-x(k)
 b(2,6)=x(j)-x(i)
 b(3,1)=b(2,2)
 b(3,2)=b(1,1)
 b(3,3)=b(2,4)
 b(3,4)=b(1,3)
 b(3,5)=b(2,6)
 b(3,6)=b(1,5)

 do i1=1,3
  do j1=1,6
   b(i1,j1)=0.5/ae*b(i1,j1)
  end do
 end do

return
end

!##################################################################################################
subroutine mf(nj,ne,npj,npj0,n,t,v,lnd,x,y,pj,p)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj),pj(npj,3),p(n)

 do i=1,n
  p(i)=0.0
 end do
 
 if (npj0/=0) then
  do i=1,npj
   ii=pj(i,1)
   p(2*ii-1)=pj(i,2)
   p(2*ii)=pj(i,3)
  end do
 end if

 if (v/=0) then
  do ie=1,ne
   call ma(ie,nj,ne,lnd,x,y,ae)
   pe=-v*ae*t/3.0
   do i=1,3
    ii=lnd(ie,i)
    p(2*ii)=p(2*ii)+pe
   end do
  end do
 end if

return
end

!##################################################################################################
subroutine rkr(nz,nd,n,jz,ak,p)
!##################################################################################################
dimension p(n),jz(nz,3),ak(500,100)

 do i=1,nz
  ir=jz(i,1)
  do j=2,3
   if (jz(i,j)/=0) then
    ii=2*ir+j-3
    ak(ii,1)=1.0
    do jj=2,nd
     ak(ii,jj)=0.0
    end do
    if (ii>nd) jo=nd
    if (ii<=nd) jo=ii
    do jj=2,jo
     ak(ii-jj+1,jj)=0.0
    end do
    p(ii)=0.0
   end if
  end do
 end do

return
end

!##################################################################################################
subroutine slov(nj,n,nd,ak,p)
!##################################################################################################
dimension p(n),ak(500,100)

 nj1=n-1
 
 do k=1,nj1
  if (n>k+nd-1) im=k+nd-1
  if (n<=k+nd-1) im=n
  
  k1=k+1
  do i=k1,im
   l=i-k+1
   c=ak(k,l)/ak(k,1)
   iw=nd-l+1
   do j=1,iw
    m=j+i-k
    ak(i,j)=ak(i,j)-c*ak(k,m)
   end do
   p(i)=p(i)-c*p(k)
  end do
 end do
 
 p(n)=p(n)/ak(n,1)

 do i1=1,nj1
  i=n-i1
  if (nd>n-i+1) jo=n-i+1
  if (nd>n-i+1) jo=nd
  
  do j=2,jo
   k=j+i-1
   p(i)=p(i)-ak(i,j)*p(k)
  end do
  p(i)=p(i)/ak(i,1)
 end do
 
 write(8,50)
50 format(/5x,5('**'),'result of calculation',5('**')//1x,'nodal displacements'//3x,'node',5x,'x-disp.',8x,'y-disp.')

 do i=1,nj
  write(8,70) i,p(2*i-1),p(2*i)
70 format(1x,i5,2e15.6)
 end do

return
end 

!##################################################################################################
subroutine made(ne,nj,n,e,pr,lnd,x,y,p)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj),d(3,3),b(3,6),s(3,6),st(3),p(n),de(6)

 write(8,10)
10 format(/1x,'elemente stresses'/)
 call md(e,pr,d)
  do ie=1,ne
   call ma(ie,nj,ne,lnd,x,y,ae)
   call mb(ie,nj,ne,lnd,x,y,ae,b)
   do i=1,3
    do j=1,6
     s(i,j)=0.0
     do k=1,3
      s(i,j)=s(i,j)+d(i,k)*b(k,j)
     end do
    end do
   end do
   
   do i=1,3
    do j=1,2
     ih=2*(i-1)+j
     iw=2*(lnd(ie,i)-1)+j
     de(ih)=p(iw)
    end do
   end do

   do i=1,3
    st(i)=0.0
    do j=1,6
     st(i)=st(i)+s(i,j)*de(j)
    end do
   end do

   stx=st(1)
   sty=st(2)
   txy=st(3)
   ast=(stx+sty)/2
   rst=sqrt(0.25*(stx-sty)**2+txy*txy)
   stma=ast+rst
   stmi=ast-rst

   if (sty==stmi) ceta=0.0
   if (sty/=stmi) ceta=90.0-57.29578*atan(txy/(sty-stmi))

   write(8,60) ie,stx,sty,txy,stma,stmi,ceta
60 format(1x,'elemetn no.=',i5/3x,'stx=',e15.6,2x,'sty=',e15.6,2x,'txy=',e15.6/3x,'stma=',e15.6,2x,'stmi=',e15.6,2x,'ceta=',e15.6)

 end do

return
end



提示我数组越界，然后我用debug查，错误出现在   p(i)=p(i)-ak(i,j)*p(k)这一行，可是我不知道如何去改正。
请哪位大哥帮小弟看看啊，不胜感激