附件是主程序，下面是前处理程序，主要进行网格划分。最后输出应力时只能出现几个结果而且是错误的，剩下的都是*号和NaN,不知道问题出现在哪里？应该怎么解决呢？谢谢！
!输出节点编号，坐标，约束，节点荷载的程序
program zuobiao
character*64 fname1,fname2
dimension p(20000)
dimension x(1000),y(1000),jz(1000,2),mea(200,9),zhong(100,100,2)
write(*,'(a22)')'input data file:'
read(*,'(a)')fname1
write(*,'(a22)')'output data file:'
read(*,'(a)')fname2
open(1,file=fname1)
open(2,file=fname2,status='new')
read(1,11)a,b,c,tt,q
read(1,12)eo,vo
read(1,13)nf1,nf2,nf3
11 format(5f10.2)
12 format(2f10.2)
13 format(3i3)
   close(1)
   nm=1
   fangle=3.141592654/(4*nf1)
   fr=(b-a)/(4*nf2)
   do 20 i=1,2*nf1+1
   angle=(i-1)*fangle
   r=a
   nbei=1
   do 20 j=1,2*nf2+1
      zhong(j,i,1)=r*cos(angle)
      zhong(j,i,2)=r*sin(angle)
      if(j.gt.nf2*2*nbei/3)nbei=nbei+1
      r=r+nbei*fr
20 continue
   zhong(2*nf2+1,nf1+1,1)=b
   zhong(2*nf2+1,nf1+1,2)=b
   zhong(2*nf2,nf1+1,1)=(b+zhong(2*nf2-1,nf1+1,1))/2.
   zhong(2*nf2,nf1+1,2)=(b+zhong(2*nf2-1,nf1+1,2))/2.
   do 30 i=1,nf1-1
      zhong(2*nf2+1,i+1,1)=b
      zhong(2*nf2+1,i+1,2)=b*tan(i*fangle)
      zhong(2*nf2,i+1,1)=(b+zhong(2*nf2-1,i+1,1))/2.
      zhong(2*nf2,i+1,2)=(zhong(2*nf2+1,i+1,2)+zhong(2*nf2-1,i+1,2))/2.
    zhong(2*nf2+1,2*nf1+1-i,2)=b
    zhong(2*nf2+1,2*nf1+1-i,1)=zhong(2*nf2+1,i+1,2)
    zhong(2*nf2,2*nf1+1-i,2)=(b+zhong(2*nf2-1,2*nf1+1-i,2))/2.
    zhong(2*nf2,2*nf1+1-i,1)=(zhong(2*nf2+1,2*nf1+1-i,1)+zhong(2*nf2-1,2*nf1+1-i,1))/2.
30 continue
   nn=1
   do 40 i=1,2*nf2+1
      if(i/2.eq.i/2.)then
         k=2
       else
         k=1
      endif
      do 40 j=1,2*nf1+1,k
         x(nn)=zhong(i,j,1)
         y(nn)=zhong(i,j,2)
         nn=nn+1
40 continue
   ne=1
do 50 i=1,nf2
  do 50 j=1,nf1
     k1=(3*nf1+2)*(i-1)+j*2-1
     k2=k1+nf1*2-j+2
     k3=k1+3*nf1+2
     mea(ne,1)=k1
     mea(ne,2)=k3
     mea(ne,3)=k3+2
     mea(ne,4)=k1+2
     mea(ne,5)=k2
     mea(ne,6)=k3+1
     mea(ne,7)=k2+1
     mea(ne,8)=k1+1
     ne=ne+1
50 continue
   nc=1
   do 60 i=0,nf2
      jz(nc,1)=1+i*(3*nf1+2)
      jz(nc,2)=2
      jz(nc+1,1)=2*nf1+1+i*(3*nf1+2)
      jz(nc+1,2)=1
      nc=nc+2
60 continue
   do 70 i=0,nf2-1
      jz(nc,1)=2*nf1+2+i*(3*nf1+2)   
      jz(nc,2)=2
      jz(nc+1,1)=3*nf1+2+i*(3*nf1+2)
      jz(nc+1,2)=1
      nc=nc+2
70 continue
   do 80 i=1,nf1+1
      zhong(1,i,1)=b
      zhong(1,i,2)=zhong(2*nf2+1,i,2)
80 continue
  flong=(c-b)/(2*nf3)
  do 90 i=2,nf3*2+1
     do 90 j=1,nf1+1  
     zhong(i,j,1)=b+(i-1)*flong
     zhong(i,j,2)=zhong(1,j,2)
90 continue
   nn1=nn
   do 100 i=2,2*nf3+1
       if(i/2.eq.i/2.)then
         k=2
       else
         k=1
       endif
      do 100 j=1,nf1+1,k
         x(nn)=zhong(i,j,1)
         y(nn)=zhong(i,j,2)
         nn=nn+1
100 continue
  do 110 i=1,nf3
   do 110 j=1,nf1/2
     if(i.eq.1)then
       k4=nn1-2*nf1-1+(j-1)*2
       mea(ne,1)=k4
       mea(ne,4)=k4+2
       mea(ne,8)=k4+1
       mea(ne,5)=nn1+(j-1)
       mea(ne,7)=nn1+j
       mea(ne,2)=nn1+nf1/2+1+(j-1)*2
       mea(ne,6)=nn1+nf1/2+2+(j-1)*2
       mea(ne,3)=nn1+nf1/2+3+(j-1)*2
       ne=ne+1
     else
       k1=(3*nf1/2+2)*(i-1)+j*2-3-nf1+nn1
       k2=k1+nf1-j+2
       k3=k1+3*nf1/2+2
       mea(ne,1)=k1
       mea(ne,2)=k3
       mea(ne,3)=k3+2
       mea(ne,4)=k1+2
       mea(ne,5)=k2
       mea(ne,6)=k3+1
       mea(ne,7)=k2+1
       mea(ne,8)=k1+1
       ne=ne+1
     endif
110 continue
   do 120 i=1,nf3
      jz(nc,1)=nn1+(i-1)*(3*nf1/2+2)
      jz(nc,2)=2
      jz(nc+1,1)=nn1+nf1/2+1+(i-1)*(3*nf1/2+2)
      jz(nc+1,2)=2
      nc=nc+2
120 continue
   nc=nc-1
   ne=ne-1
   nn=nn-1
   nn1=nn
   do 500 i=1,20000
   p(i)=0
500 continue
    do 130 i=1,nf1
    p((nn1-1)*2+1)=p((nn1-1)*2+1)+(y(nn1)-y(nn1-1))*q*0.5
    p((nn1-2)*2+1)=p((nn1-2)*2+1)+(y(nn1)-y(nn1-1))*q*0.5
      nn1=nn1-1
130 continue
    do 140 i=1,ne
       mea(i,9)=1
140 continue
    write(2,460)nn,ne,nm,nc
    write(2,465)(i,x(i),y(i),p(i*2-1),p(i*2),i=1,nn)
    write(2,470)((jz(i,j),j=1,2),i=1,nc)
    write(2,475)(i,(mea(i,j),j=1,9),i=1,ne)
    write(2,480)eo,vo,tt
460 format(4i5)
465 format(i5,4f15.6)
470 format(5x,i5,5x,i5)
475 format(5x,i5,5x,9i5)
480 format(3f15.6)
    end