声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2506|回复: 2

[Fortran] 求助修改程序

[复制链接]
发表于 2006-5-20 16:29 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
<P>以下的程序是用for90编的,但是有错误,我怎么都改不对,还请高手指点,先谢谢了!<BR><BR>implicit none<BR>integer b,i,j,m<BR>integer s3,b3,b2,l0,m0,n0,ln,ss                 <BR>real gn,m1,ml,ln1,l,a1,a12,a22,a21,a11,m21,m22,m11,m12,a01,a02,b1,&amp;<BR>     c0,rs,rt,rxo,ia,cc1,aa1,aa2,ra,bz,bm,sp,kkk,e,f,a2,l1<BR>external rn,fn,fk,evs,optsol<BR>real,DIMENSION(:), ALLOCATABLE::r,z,n,u,c1,c2,c3,c4,c5,c6,c7,c8,c9,k,qx<BR>real,DIMENSION(:,:), ALLOCATABLE::p,kk  <BR>m0=5;n0=5;ln=5;ss=5<BR>allocate(r(5),z(5),n(5),u(10),c1(5),c2(5),c3(5),c4(5),&amp;<BR>       c5(10),c6(10),c7(10),c8(10),c9(10),k(5),p(5,3),kk(10,2),qx(5))<BR>      ! (r(m0),z(n0),n(ln),u(n0+1:2*n0),c1(l0),c2(l0),c3(l0),c4(l0),&amp;<BR>        ! c5(l0),c6(l0),c7(l0),c8(l0),c9(l0),k(ss),b(5) ,p(ln,3),kk(l0,2),qx(5))<BR>print*,'input s3'<BR>read*,s3<BR>!***********************read有问题***********************************************<BR>!gn分块求解的总块数减1,m1每一块处理的非公用网线的条数<BR>!l0电极系的套数,n0横向网线的最大编号数,ml最后一块处理的纵向网线总条数减1<BR>!read(4,gn,m1,n0,l0,ml)<BR>print*,'inputgn,m1,n0,l0,ml'<BR>read*,gn,m1,n0,l0,ml      <BR>     m0=gn*(m1+1)+ml+1<BR>  ln=(m1+2)*n0<BR>  ln1=(ml+1)*n0<BR>  ss=n0*(n0+1)*(m1+1.5)<BR>!read(4,b2,b3)<BR>print*,'input b2,b3'<BR>read*,b2,b3<BR>!read(4,r,z,c1,c2,c3,c4,c5,c6,c7,c8,c9,kk)<BR>print*,'input  r,z,c1,c2,c3,c4,c5,c6,c7,c8,c9,kk '<BR>read*,r,z,c1,c2,c3,c4,c5,c6,c7,c8,c9,kk<BR>     do l=0,m0<BR>        r(l)=r(l)/1000<BR>  enddo<BR>  do l=0,n0<BR>     z(l)=z(l)/1000<BR>  enddo<BR>print '(5x,5i5)', gn,m1,n0,l0,ml<BR>print '(5x,12i15.7)', r,z,c1,c2,c3,c4,c5,c6,c7,c8,c9,kk<BR>     if(s3==1)then                        !控制台变量(s3)若为1计算视电阻率<BR>     rs=1<BR>  rt=rs<BR>  rxo=rt<BR>  b1=5<BR>  c0=b1<BR>  else<BR>     ! read(4,rs,rt,rxo,c0,b1)<BR>  print*,'input rs,rt,rxo,c0,b1'<BR>  read*,rs,rt,rxo,c0,b1<BR>   print '(5x,7i15.7)',b1,b2,b3,c0,rs,rt,rxo<BR>   call fn(ln,0)<BR>   rt=1./rt<BR>   rs=1./rs<BR>   rxo=1./rxo<BR>   cc1=5<BR>   call fk(1,n0,1)<BR>   call optsol(k,p,n,3,2,n0,1)<BR>   do i=1,gn       !从第一块到第gn块刚度阵的逐块形成及分解,公共三角块的移置<BR>      a1=(i-1)*(ln-n0)+n0+1<BR>   a2=i*(ln-n0)+n0<BR>   call fk(a1,a2,n0+1)<BR>   call optsol(k,p,n,3,n0+1,ln,1)<BR>   print '(5x,i5)',i<BR>   do l=1,n0<BR>      do bm=1,l<BR>         a1=(m1+1)*n0+l<BR>      k(n(l)-l+bm)=k(n(a1)-l+bm)<BR>      enddo<BR>   enddo<BR>   enddo<BR>!   endif</P>
<P>  !******************************************************************<BR>  !计算电极边界节点在最后一块里相应的行号;计算最后一块节点相应的主对角元地址;形成系数阵<BR>         aa1=gn*(ln-n0)+n0+1<BR>   aa2=m0*n0<BR>   do l1=1,l0<BR>      a21=n0+c1(l1)<BR>   a22=n0+c2(l1)<BR>   a11=n0+c3(l1)<BR>   a12=n0+c4(l1)<BR>   m21=n0+c5(l1)<BR>   m22=n0+c6(l1)<BR>   m11=n0+c7(l1)<BR>   m12=n0+c8(l1)<BR>   a01=n0+c9(l1)<BR>   a02=2*n0<BR>   call fn(ln1,1)<BR>   cc1=c1(l1)<BR>   call fk(aa1,aa2,n0+1)<BR>   enddo<BR>   do i=1,ml-1      !对仪器棒内诸节点所对应的主对角元系数进行处理<BR>      j=ln1-(i-1)*n0<BR>   do m=j,j-n0+cc1+1 ,-1<BR>      k(n(m))=1<BR>   enddo<BR>   enddo<BR>   call evs(a21,a22)     !对电极进行等位面处理<BR>   call evs(m21,m22)<BR>   call evs(m11,m12)<BR>   call evs(a11,a12)<BR>   call evs(a01,a02)<BR>   do l=1,ln       !形成最后一块的右端向量<BR>      do i=1,3<BR>      p(l,i)=0<BR>      p(a22,3)=3.14159265/4<BR>      p(a12,2)=p(a22,3)<BR>      p(a02,1)=p(a12,2)<BR>   enddo<BR>   enddo<BR>   call optsol(k,p,n,3,n0+1,ln1,3)           !对最后一块进行三角分解和回代<BR>   do  l=1,3                                       !打印电极各分场的电位值<BR>       qx(1)=a22<BR>    qx(2)=a12<BR>    qx(3)=m22<BR>    qx(4)=m12<BR>    qx(5)=a02<BR>    do j=1,5<BR>       sp=p(qx(j),l)<BR>    print '(5x,f10.4)','sp=',sp<BR>    enddo<BR>   enddo<BR>   !进行浅双测向的计算并输出有关值<BR>   e=p(m12,1)-p(m12,3)-p(m22,1)+p(m22,3)   <BR>   f=p(m22,2)-p(m22,3)-p(m12,2)-p(m12,3)<BR>   ia=e/f<BR>   do l=n0+1,2*n0<BR>      u(l)=p(l,1)-p(l,3)+ia*(p(l,2)-p(l,3))<BR>   enddo<BR>   print*,u<BR>   if(s3==1)then      !这个地方问题不会********88<BR>      kk(l1,1)=1/u(m12)<BR>   kkk=kk(l1,1)<BR>   print '(5x,i5,2f15.7)',l1,kkk,ia<BR>   else<BR>      ra=kk(l1,1)*u(m12)<BR>   kkk=kk(l1,1)<BR>   bz=1/(rt*ra)<BR>   print '(5x,f15.7)',bz<BR>   print '(5x,i5,3f15.7)',l1,ra,kkk,ia<BR>   endif<BR> endif<BR> deallocate(r,z,n,u,c1,c2,c3,c4,c5,c6,c7,c8,c9,k,p,kk,qx)</P>
<P>end<BR>!*********************主程序结束**************************************</P>
<P>!*********************进行外部子程序**********************************</P>
<P>!*********************************************************************<BR>!**本程序:为等位面处理的需要在节点编号从lb到le的等位面上,对第le行的非零宽度加大<BR>subroutine rn(lb,le)<BR>integer lb,le<BR>real,DIMENSION(:), ALLOCATABLE::n<BR>allocate(n(:))<BR>  a=n(le)<BR>  do l=lb,le-1<BR>     a1=le-l+n(l)<BR>     if(a&lt;a1)a=a1<BR>  enddo<BR>  n(le)=a<BR>deallocate(n)<BR>end</P>
<P>!**形成一维紧缩排列下刚度阵的主对角元的地址****************************<BR>subroutine fn(ll,ind)<BR>integer ll,ind<BR>real,DIMENSION(:), ALLOCATABLE::n<BR>allocate(n(:))</P>
<P>do l=1,ll                       !计算刚度阵各行的非零宽度<BR>   if(l&lt;=n0)then<BR>      n(l)=l-1<BR>   else<BR>      n(l)=n0<BR> endif<BR>enddo<BR>if(ind==0)goto 3<BR>call rn(a21,a22)      !调用rn过程<BR>call rn(a11,a12)<BR>call rn(m21,m22)<BR>call rn(m11,m12)<BR>call rn(a01,a02)<BR>3  n(0)=0        !形成刚度阵的主对角元地址<BR>   do l=1,ll<BR>      n(l)=n(l)+n(l-1)+1<BR> enddo<BR>deallocate(n)<BR>end</P>
<P>!**形成并存储刚度阵的过程*************************************************<BR>SUBROUTINE fk(lb,le,lsb)<BR>integer lb,le,lsb<BR>real aa,bb,cc,dd,oo,f1,f2,f3,f4,f5,f6,f7,f8,he,hs,hw,hn,s,t<BR>real,DIMENSION(:), ALLOCATABLE::  k<BR>allocate(k(:))<BR>!将刚度阵从lsb行开始充零<BR>do l=lb,le<BR>   a1=lsb-lb+l<BR>   do a2=n(a1-1)+1,n(a1)<BR>      k(a2)=0<BR>   enddo<BR>enddo<BR>!计算节点所在纵横向网线的序号<BR>do l=lb,le<BR>   s=(l-0.1)/n0+1<BR>   t=l-(s-1)*n0</P>
<P>!计算节点四周四个象限的电阻率<BR>if(s&lt;=b2.and.t&lt;=c0)then<BR>     ss1=rs<BR>elseif(s&lt;=b1.and.c0&lt;t)then<BR>     ss1=rt<BR>elseif(b1&lt;s.and.s&lt;=b2.and.c0&lt;t)then<BR>     ss1=rxo<BR>elseif(b2&lt;s.and.s&lt;=b3.or.b3&lt;s.and.t&lt;cc1)then<BR>     ss1=1<BR>else<BR>     ss1=0<BR>endif<BR>if(s&lt;=b2.and.t&lt;=c0)then<BR>     ss2=rs<BR>elseif(s&lt;=b1.and.c0&lt;t)then<BR>     ss2=rt<BR>elseif(b1&lt;=s.and.s&lt;b2.and.c0&lt;t)then<BR>     ss2=rxo<BR>elseif(b2&lt;=s.and.s&lt;b3.or.b3&lt;=s.and.s&lt;=m0-1.and.t&lt;=cc1)then<BR>     ss2=1<BR>else<BR>     ss2=0<BR>endif<BR>if(s&lt;b2.and.t&lt;c0)then<BR>     ss3=rs<BR>elseif(s&lt;b1.and.c0&lt;=t.and.t&lt;n0)then<BR>     ss3=rt<BR>elseif(b1&lt;=s.and.s&lt;b2.and.c0&lt;=t.and.t&lt;n0)then<BR>     ss3=rxo<BR>elseif(b2&lt;=s.and.s&lt;b3.and.t&lt;n0.or.b3&lt;=s.and.s&lt;=m0-1.and.t&lt;cc1)then<BR>     ss3=1<BR>else<BR>     ss3=0<BR>endif<BR>if(s&lt;=b2.and.t&lt;c0)then<BR>     ss4=rs<BR>elseif(s&lt;=b1.and.c0&lt;=t.and.t&lt;n0)then<BR>     ss4=rt<BR>elseif(b1&lt;s.and.s&lt;=b2.and.c0&lt;=t.and.t&lt;n0)then<BR>     ss4=rxo<BR>elseif(b2&lt;s.and.s&lt;=b3.and.t&lt;n0.or.b3&lt;s.and.t&lt;cc1)then<BR>     ss4=1<BR>else<BR>     ss4=0<BR>endif<BR>!计算节点的r坐标、步长he,hs,hn,hw以及有关的因子<BR>rrs=r(s)<BR>he=r(s-1)-r(s)<BR>if(t==n0)then<BR>    hs=1<BR>else<BR>    hs=z(t)-z(t+1)<BR>endif<BR>if(s==m0)then<BR>    hw=1<BR>else<BR>    hw=r(s)-r(s+1)<BR>endif<BR>hn=z(t-1)-z(t)<BR>f1=(3*rrs+2*he)*ss4/(he*6)<BR>f2=(3*rrs+he)*ss1/(he*6)<BR>f3=(3*rrs+he)*ss4/(hs*6)<BR>f4=(3*rrs-hw)*ss3/(hs*6)<BR>f5=(3*rrs-hw)*ss3/(hw*6)<BR>f6=(3*rrs-2*hw)*ss2/(hw*6)<BR>f7=(3*rrs+he)*ss1/(hn*6)<BR>f8=(3*rrs-hw)*ss2/(hn*6)<BR>!计算节点及四邻点有关的系数<BR>aa=f1*hs+f2*hn<BR>bb=f3*he+f4*hw<BR>cc=f5*hs+f6*hn<BR>dd=f7*he+f8*hw<BR>oo=aa+bb+cc+dd<BR>!存储所得的系数阵<BR>a1=lsb-lb+l<BR>a2=n(a1)<BR>k(a2)=oo<BR>if(1&lt;t)then<BR>      a3=a2-1<BR>   k(a3)=-dd<BR>endif<BR>if(1&lt;s)then<BR>      a3=a2-n0<BR>   k(a3)=-aa<BR>endif<BR>enddo<BR>deallocate(k)<BR>end<BR>!****等位面处理过程******************************************<BR>subroutine evs(lb,le)<BR>integer lb,le<BR>real,DIMENSION(:), ALLOCATABLE::k<BR>allocate(k(:))<BR>!对等位面上从lb到le节点对应的主对角元进行处理<BR>aaa=k(n(lb))<BR>k(n(lb))=1<BR>do l=lb+1,le<BR>   a1=n(l)<BR>   a2=a1-1<BR>   aaa=aaa+k(a1)+2*k(a2)<BR>   k(a1)=1<BR>   k(a2)=0<BR>enddo<BR>k(n(le))=aaa<BR>!对等位面上从lb到le节点对应的行的元素进行处理<BR>do l=lb,le-1<BR>   a1=n(l-1)+1<BR>   a2=n(le)-le+l-n0<BR>   k(a2)=k(a2)+k(a1)<BR>   k(a1)=0<BR>enddo<BR>a1=n(lb)-1<BR>a2=n(le)-le+lb-1<BR>k(a2)=k(a2)+k(a1)<BR>k(a1)=0<BR>!对等位面上从lb到le节点对应的列的元素进行处理<BR>do l=lb,le-1<BR>   a1=n(l+n0-1)+1<BR>   a2=n(l+n0)-l+le-n0<BR>   k(a2)=k(a2)+k(a1)<BR>   k(a1)=0<BR>enddo<BR>deallocate(k)<BR>end<BR>!****Crout三角分解的标准过程**********************<BR>subroutine optsol(a,b,na,m,nb,neq,kex)<BR>!real,DIMENSION(:,:), ALLOCATABLE::  b<BR>!real,DIMENSION(:), ALLOCATABLE::a,na<BR>real najp,j,il,jk,naj,fi,fi1,jia,naip,jka,kl,i,nai,ik,ii,kf,jj,aa,k,cc,bb,p<BR>integer nb,neq,m,kex<BR>real a(:),b(:,:),na(:)<BR>allocate(a(:),b(:,:),na(:))<BR>if(kex==0) goto 110<BR>10 najp=na(nb-1)<BR>do j=nb,neq<BR>   il=j-1<BR>   naj=na(j)<BR>   if(kex==2.or.kex==5)goto 100<BR>   jk=naj-j<BR>   fi=1-jk+najp<BR>   if(j&lt;=fi)goto 80<BR>   fi1=fi+1<BR>   if(il&lt;fi1)goto 60<BR>   jia=2+najp<BR>   naip=na(fi)<BR>   kl=fi+jk<BR>   do i=fi1,il<BR>      nai=na(i)<BR>   ik=nai-i<BR>   ii=i+1-nai+naip<BR>   if(i&lt;=ii)goto 40<BR>   if(ii&lt;=fi)then <BR>      kf=fi+jk<BR>   else<BR>      kf=ii+jk<BR>   endif<BR>   jj=ik-jk<BR>   aa=0<BR>   do k=kf,kl<BR>      aa=aa+a(k)*a(jj+k)<BR>   enddo<BR>   a(jia)=a(jia)-aa<BR>40    jia=jia+1<BR>   kl=kl+1<BR>   naip=nai<BR> enddo<BR>60 kf=jk+fi<BR>    kl=naj-1<BR>100 if(kex==1)goto 70<BR>    ! if(kex==2.or.kex==5)then<BR>      do k=fi,j-1<BR>      do p=1,m<BR>      b(j,p)=b(j,p)-a(k+jk)*b(k,p)*a(na(k))<BR>    enddo<BR>   enddo<BR>  if(kex==2) goto 80<BR>70     aa=0 <BR>       do k=kf,kl<BR>      nai=na(fi)<BR>   cc=a(k)/a(nai)<BR>   aa=aa+a(k)*cc<BR>   a(k)=cc<BR>   fi=fi+1<BR>  enddo<BR>       a(naj)=a(naj)-aa<BR>80  najp=naj<BR>    if(kex==1)goto 90<BR>     do p=1,m<BR>     b(j,p)=b(j,p)/a(naj)<BR>  enddo<BR>90 enddo !j循环结束<BR>     if(kex==1.or.kex==2.or.kex==4)goto 200<BR>110  do j=neq,nb,-1        !分块解法的回代<BR>        ii=j-na(j)+na(j-1)+1<BR>  jka=na(j-1)+1<BR>  if(j&lt;=ii)goto 180<BR>  kl=j-1<BR>  do k=ii,kl<BR>     do p=1,m<BR>        bb=b(j,p)<BR>     b(k,p)=b(k,p)-a(jka)*bb<BR>     enddo<BR>     jka=jka+1  <BR>  enddo<BR>180  enddo<BR>  !   deallocate(a,b,na)<BR>200 end</P>



<P><BR> </P>
回复
分享到:

使用道具 举报

 楼主| 发表于 2006-5-21 15:00 | 显示全部楼层
怎么没有人帮我呀!
发表于 2006-5-22 13:55 | 显示全部楼层
<P>你的<BR>allocate<BR>语句用得有问题<BR>allocate应该在你明确了动态数组的大小(维数)后再开辟。<BR>例如:<BR><BR>integer,allocatalbe :: N(:) !定义动态数组N<BR>......<BR>!!得到了N维数N1<BR>allocate(N(N1))  !开辟N的大小<BR>......<BR><BR></P>
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-25 10:50 , Processed in 0.073136 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表