求助修改程序

小飞侠 · 发表于 2006-5-20 16:29

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

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

x

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

小飞侠 · 发表于 2006-5-21 15:00

怎么没有人帮我呀!

Silence · 发表于 2006-5-22 13:55

你的 allocate 语句用得有问题 allocate应该在你明确了动态数组的大小(维数)后再开辟。 例如： integer,allocatalbe :: N(:) !定义动态数组N ...... !!得到了N维数N1 allocate(N(N1)) !开辟N的大小 ......

账号		自动登录	找回密码
密码			我要加入

[Fortran] 求助修改程序

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