材料学中,运用fortran计算疲劳裂纹扩展速率(毕设中用)。希望哪位高手帮忙修改下。

c  fatigue crack growth rates programme
c  this programme file name-GB.FOR
cc m-total test points
cc Pmax-maximun load
cc Pmin- minimun load
cc b-specimen thickness
cc w-specimen widthness
cc a0-initial crack length
cc ys-yield strength
cc Sratio-Stress ratio
cc Ni(I)- load cycle number
cc Ai(I)-a set of crack length
1 dimension ni(100),ai(100),dadn(100),
2 &delk(100),aa(10),nn(10),bb(3),name(20)
3  integer qq
4  parameter (PI=3.1415926)
5  Character*15 finname,foutname,key*5,
6  &ancurve*15,dkdadn*15,name*50
7  write(*,*)'input Title:'
8  read(*,'(A)')name(8)
9  write(*,*)'input data file name:'
10 read(*,'(A)')finname
11 write(*,*)'output data file name:'
12 read(*,'(A)')foutname
13 write(*,*)'enter A-N curve name:'
14 read(*,'(A)')ancurve
15 write(*,*)'enter DELT K-da/dN curve name:'
16 read(*,'(A)')dkdadn
17 open(15,file=finname,status='old')
18 open(16,file=foutname,status='new')
19 open(17,file=ancurve,status='new')
20 open(18,file=dkdadn,statu='new')
21 write(*,*)'Material Name:'
22 read(*,'(A)')name(1)
23 write(*,*)'No.of Specimen'
24 read(*,'(A)')name(2)
25 write(*,*)'environment'
26 read(*,'(A)')name(3)
27 write(*,*)'test frequency'
28 read(*,'(A)')name(9)
29 write(*,*)'Temperature'
30 read(*,'(A)')name(4)
31 write(*,*)'Moisture'
32 read(*,'(A)')name(5)
33 write(*,*)'Test Date'
34 read(*,'(A)')name(6)
35 write(*,*)'Testers'
36 read(*,'(A)')name(7)
37 write(16, 300)name(8)
38 300 format(1x,'******TEST REPORT:',A35,'******')
39 write(16,299)
40 299 format(/)
41 write(16,400)name(1)
42 format(1X,'Material Name & Orientation:',2x,A40)
43 write (16,410)name(2)
44 410 format(1x,'No.of Specimen:',2x,A20)
45 write (16,420)name(3)
46 420 format(1x,'environments:',2x,A20)
47 write (16,430)name(9)
48 430 format(1x,'Test Frequency:',3x,A3,'Hz')
49 write (l6,440)name(4)
50 440 format(1x,'Temperature:',2x,A3,'c')
51 write(16,450)name(5)
52 450 format(1x,'Moisture:',2x,A4,'%','RH')
53 write(16,460)name(6)
54 460 format(1x,'Test Date:',2x,A15)
55 write(16,470)name(7)
56 470 format(1x,'Testers:',2x,A20)
57 write(*,*)'asking for specimen kind:SEB,CT OR MT'
58 read(*,'(A)')key
59 c begin to operate seven points poly
60 read(15,*)m,pmax,pmin,b,w,ys,a0,sratio
61 write(16,20)
62 20 format(/)
63 write(16,21)
64 21 format(5x,'========',1x,'RESULTS OF TEST DATA PROCESSING',
65 &1x,'========')
66 write(16,22)
67 22 format(/)
68 write(16,25)m,pmax,pmin, b, w ,a0 ,ys ,Sratic,
69 25 format (1x,'the test points:',i4,2x,'pmax=',
70 &f8.1,'N',2x,'pmin=',f8.1,'N',/,1x,'b=',f5.2,'mm',2x,'w=',
71 &f6.2,'mm',2x,'A0=',f6.3,'mm',/1x,'yield limit ys=',f8.2 ,
72 &'MPa',5x,'Stress Ratio R=',f6.3)
73 read(15,*)(Ni(I),Ai(I),I=1,m )
74 DO 30 I=1,m
75 AI(I)=A0+A(I)
76 30 Continue
77 write(16,105)
78 105 format(1x,'cbs.no',1x,'cycle',2x,'a(means)',2x,'a(areg.)'
79 &,3x,'m.c.c.',4x,'delk',6x,'dada')
80 k=o
81 pp=pmax-pmin
82 r=pmin/pmax
83 write (16,95)(i,Ni(i),Ai(I),I=1,3)
84 m=m-6
85 do 100 i=1,m
86 1=0
87 k=k+1
88 k1=k+6
89 do 64 j=k,k1
90 1=1+1
91 aa(1)- ai(j)
92 nn(1)一ni(j)
93 64 continue
94 c1=0.5*(nn(1)+nn(7))
95 c2=0.5*(nn(7)-nn(1))
96 sx=0.0
97 sx2=0.0
98 sx3=0.0
99 sx4=0.0
100 sy= 0.0
101 syx=0.0
102 syx2=0.0
103 do 70 j=1,7
104 x=(NN(j)- c1)/c2
105 yy=aa(j)
106 sx=sx+x
107 sx2=sx2+x**2
108 sx3=sx3+x**3
109 sx4=sx4+x**4
110 sy=sy+yy
111 syx=syx+x*yy
112 sYx2=syx2十yy* (x**2)
113 70 continue
114 den=7.0*(sx2*sx4-sx3**2)-sx*(sx*sx4-sx2*sx3)+sx2*(sx*sx3
115 &-sx2**2)
116 t2=sy*(sx2*sx4-sx3**2)-syx*(sx*sx4-sx2*sx3)+syx2*(sx*sx3
117 &-sx2**2)
118 bb(1)=t2/dem
119 t3=7.0*(syx*sx4一syx2*sx3)-sx*(sy*sx4-syx2*sx2)+
120 &sx2*(sy*sx3-syx*sx2)
121 bb(2)=t3/den
122 t4=7.0*(sx2*syx2-sx3*syx)-sx*(sx*syx2-sx3*sy)+sx2*(sx*syx
123 &-sx2*sy)
124 bb(3)=t4/den
125 yb=sy/7.0
126 rss=0.0
127 tss=0.0
128 do 75j=1,7
129 x=(nn(j)-cl)/c2
130 yhat=bb(1)+bb(2)*x+bb(3)*(x**2)
131 rss=rss+(aa(j)- yhat)**2
132 tss=tss+(aa(j)一yb)**2
133 75 continue
134 r2=1.0-rss/tSS
135 dadn(i)=(bb(2)/c2+2.0*bb(3)*(nn(4)-c1)/(c2**2))
136 x=(nn(4)-cl)/c2
137 ar=bb(1)+bb(2)*x+bb(3)*(x**2)
138 s=1el0
139 snet=0.0
140 qq=i+3
141 IF (key.se.'SEB'.or .key.eq .'seb')THEN
142 t=ar/w
143 ft=6*t**0.5/(1.0 +2.0 *t)/(1.0-t)**1.5*
144 &(1.99-t*(1-t)*(2.15-3.93*t 12.7*t**2))
145 delk(i)=((ft*pp(/(b*sqrt(w)))/(10**1.5)
146 s=ys*(w*(1一t))**2
147 ax=pmax*6.0*w/h
148 snet=pmax/(b*w*(l.-t)
149 else if(key.eq.'ct'.or.key.eq.'CT')then
150 t=ar/w
151 ft=(2+t)*(0.8866+4.64*t-13.32*t*t+14.72*t**3  5.6*t**4
152 &/((1-t**1.5)
153 delk(i)-(ft*pp)/(b*sqrt(w))/10**1.5
154 s=ys*sqrt(pi*w*(1-t))/2
155 sent=pmax/(b*w*(1.-t))
156 ax=delk(I)/(1.-r)
157 else
158 if(key.EQ.'MT'.or.key.eq.'mt')then
159 t=2.*ar/w
160 sec=1./(cos(pi*t/2.))
161 ft=sqrt((pi*t*sec)/2.)
162 snet=pmax/(b*w*(1.-t)
163 delk(i)=(ft*pp)/(B*sqrt(w))/10**1.5
164 ax=delk(I)/(1.-t)/1.25
165 end if
166 end if
167 if(ax.ge.s)go to 398
168 if(snet.ge.ys)go to 398
169 392 write(16,92)qq,Ni(qq),Ai(qq),ar,r2,delk(i),dadn(i)
170 go to 100
171 398 write (16,98)qq,Ni(qq),Ai(qq),ar,r2,delk(i),dadn(i)
172 100 continue
173 92 format(i4,lx,i8,2x,f7.3,2x,f7.4,2x,f7.5,2x,f8.4,2x,e10.4)
174 98 format(i4,lx,i8,2x,f7.3,2x,f7.4,2x,f7.5,2x,f7.4,2x,e10.4,
175 2x,'***')
176 j=m+4
177 k=m+6
178 write(16,95)(i,Ni(i),Aai(i),i=j,k)
179 format(14,1X,18,2X,F7.3)
180 write(17,200(Ni(i),Ai(i),i=l,m+6)
181 200 format (1x,18,5X,F6.3 )
182 write(18,350)(delk(i),dadn(i),i=l,qq-3)
183 350 format(1x,F8.5,5X,E12.5)
184 close(17)
185 close(18)
186 close(16)
187 close(15)
188 end