컴파일이나 빌드시 오류가 없는데 결과에서 오류가 뜹니다..

[김현수]

프로그램으로 논문을 쓰고 있는데...계속해서 오류가 뜨고 어렵네요..
며칠째 오류잡는다고 모르는 책 뒤져가며..쉽지 않네요..

컴파일이나 빌드시에 오류는 없었는데요..
결과치가 맨 아래처럼 나오더군요..
main compute부분에서의 값이 계산되서 나와야 되는건 아닌가 싶은데요...
시간은 없고 해결은 안되고..해서 또다시 올립니다..
죄송하구요..보시고 조그마한 한마디라도 부탁드립니다...
아침에도 정말 많은 도움이 됬었는데요...부탁드려요...

program servo_1
cc..This solves system of equations describing dynamics of
cc..an electrohydrodynamic servo valve.
cc..All the physical variables are in SI(MKS) unit,
cc..except the final output values.
cc..If the program shows an overflow, increase imax,
cc..or decrease tmax.
implicit real*8 (a-h,o-z)
dimension f(1001),g(1001)
open(10,file='oservo')
open(12,file='ottt')
cc..Constants.
pi=2.d0*dasin(1.d0)
imax=1000001
tmax=0.2
dt=tmax/(imax-1.)
iw=(imax-1)/200
cc..Input data.
**..ps: supply pressure from a pump.
**..paa: pressure from load supplied through a port.
ps=100.e05
paa=ps/2.
cc..Fluid properties.
**..rmu: dynamic viscosity, rho: density.
rmu=25e-03
rho=850.
cc..Torque motor spec.
**..rk: bent-spring constant, rki: ampere constant,
**..rjfp: moment of rotational inertia, cfp: damping coeff,
**..r, b: distance of the nozzle and spring_end of spool
**..from the pivot point of the flapper plate's rotation,
**..xnm: distance between plate and nozzle at neutral.
rk=30000.
rki=1.6
rjfp=5.e-06
cfp=8.
c cfp=40.
r=19.e-03
b=45.e-03
rkcoma=0
rkfp=rkcoma+b**2*rk
xnm=0.04E-03
cc..Nozzle-orifice spec.
**..dor: orifice diameter, dn: nozzle diameter,
**..dd: pipe diameter(in fact, d_1 and d_2 in the text),
**..ao, an, aa: areas,
**..cqo, cqn: discharge coefficient of orifice and nozzle,
**..cc: contraction coefficient of the nozzle flow,
**..cf, z: constants(see text).
c dor=0.7e-03
dor=1.3e-03
ao=pi*dor**2/4.
dn=3.2e-03
an=pi*dn**2/4.
dd=7.e-03
aa=pi*dd**2/4.
cqo=0.67
cqn=0.60
cc=cqn
cf=2.*cc/(1.-(cc*an/aa)**2)
z=16.*(cqn/cqo)**2*(dn/dor)**2*(xnm/dor)**2
cc..Valve_spool properties.
**..rm: spool mass, ds: sppol diameter,
**..w: width of the outlet/inlet port,
**..rlv: length of the contact between spool and housing,
**..eps: clearance between spool and housing,
**..cs: damping coefficient of spool movement,
**..cq: discharge coefficient of spool port.
**..gamma: ratio of two flow rates through ports.
rm=0.03
ds=10e-03
as=pi*ds**2/4.
w=pi*ds
rlv=40.e-03
eps=1.3e-05
cs=pi*rmu*ds*rlv/eps
c cs=cs*5.
cq=0.60
gamma=0.7
cc..Flow constants.
**..cv: velocity coefficient, thf: jet angle.
cv=1.0
thf=69.
thf=thf*pi/180.
cc..Eval‎uate effective spring constant rkse and effective
cc..damping cse of the flapper plate motion.
rkse=rk+2.*cq*cv*w*(1.+gamma**2)*(ps-paa)*dcos(thf)
cse=cs
cc..prepare the function f and g for linear interpolations.
km=101
pstart=1.d-20
pend=ps-1.d-20
dpp=(pend-pstart)/(km-1.d0)
do k=1,km
pp=pstart+dpp*(k-1.d0)
f(k)=cqo*ao*dsqrt(2.d0*(ps-pp)/rho)
g(k)=cqn*pi*dn*dsqrt(2.d0*pp/rho)
yy=f(k)-xnm*g(k)
write(12,960) pp,yy
960 format(2e13.5)
end do
cc..Get the estimated steady current for the xs desired,
cc..and set the initial current.
xss=1.e-03
rkxth1=(rk+2.*cq*cv*w*(1.+gamma**2)*(ps-paa)*dcos(thf))
rkxth2=(4.*r*z*as*ps/(xnm*(1+z)**2)-b*rk)
rkxth=rkxth1/rkxth2
rkix=rki/(b*rk+rkxth*(4.*cf*r**2*z*an*ps/(xnm*(1+z)**2)+rkfp))
ri=xss/rkix
write(*,*) ' estimated initial current for xs=',xss,' is'
write(*,*) ri, ' ampere'
write(*,*) ' Set the initial current i; '
read(*,*) ri
cc..Initialize xs, u, th, om, pa, and pb.
xs=0.
u=0.
th=0.
om=0.
pa=ps/(1.+z)
pb=pa
cc..Now, start the main computation.
do i=1,imax
if(mod(i,10000).eq.0) write(*,*) i
t=dt*(i-1.)
xsn=xs+dt*u
un=u+dt*(-cse*u-rkse*xs+as*(pa-pb)-b*rk*th)/rm
thn=th+dt*om
omn=om+dt*(-cfp*om-rkfp*th+rki*ri-b*rk*xs-r*cf*an*(pa-pb))/rjfp
**..Get pan and pbn by linear interpolations.
asu=as*un
do k=1,km-1
ff1=f(k)-(xnm-r*thn)*g(k)
ff2=f(k+1)-(xnm-r*thn)*g(k+1)
if(ff1.ge.asu.and.ff2.le.asu) go to 10
910 format(i5,3e13.5)
end do
write(*,*) ' Failed in the first interpolation.'
write(*,*) ' ffstart=',f(1)-(xnm-r*thn)*g(1)
write(*,*) ' ffend=',f(km)-(xnm-r*thn)*g(km)
write(*,*) ' asu=',asu
stop
10 pan=pstart+dpp*(k-1.)+(ff1-asu)/(ff1-ff2)*dpp
asu=-as*un
do k=1,km-1
ff1=f(k)-(xnm+r*thn)*g(k)
ff2=f(k+1)-(xnm+r*thn)*g(k+1)
if(ff1.ge.asu.and.ff2.le.asu) go to 20
end do
write(*,*) ' Failed in the 2nd interpolation.'
write(*,*) ' ffstart=',f(1)-(xnm+r*thn)*g(1)
write(*,*) ' ffend=',f(km)-(xnm+r*thn)*g(km)
write(*,*) ' asu=',asu
stop
20 pbn=pstart+dpp*(k-1.)+(ff1-asu)/(ff1-ff2)*dpp
**..Now update the variables.
xs=xsn
u=un
th=thn
om=omn
pa=pan
pb=pbn
if(mod(i,iw).eq.1) write(10,900) t,xs*1000,pa*1e-07,pb*1e-07
900 format(8e13.5)
end do
stop
end












결과치 입니다..

estimatedinitialcurrentfor xs= 1.000000047497451E-003 is
1.09302748728335 ampire
Set the initial current i;

[최권일]

전기쪽인 것 같은데... 제가 보기엔 전류값 초기치(가정치일 듯 싶습니다)를 다시 입력하면 될것 같은데요 예를 들어 1.09302를 도스창에서 다시 입력하고 엔터를 치면 계속 돌아갈 겁니다.

[김현수]

감사합니다..계속 노력해볼께요..

[꿈꾸는가인]

free format 은 되도록 사용하시지 말고 엄격하게 프로그래밍 하셨으면 합니다. 우선 숫자를 나타낼 때 지수부분에 e 를 사용하셨는데 double precision 이므로 d 라고 사용하셔야 합니다. dasin -> dsin 이 되어야 할거 같습니다.

[꿈꾸는가인]

제가 free format을 사용하지 않아서 드리는 질문입니다만 comment 는 어떻게 나타내나요? f77 에서는 첫줄에 C를 사용하면 되는데 free format 에서는 cc를 사용하나요? 아니면 여기서도 c 만 쓰면 comment 를 의미하는 것인지요?

[꿈꾸는가인]

위의 프로그램을 보면 cc 또는 ** 가 comment 를 나타내는 거 같은데 맞습니까? 그렇다면 프로그램 중간부분에 cc = cqn 이라고 되어서 다음 식에 변수로 들어가는데 이 경우는 comment 로 간주 하나요? 아니면 cc = cqn 이라는 식으로 인식을 하나요? free format 을 사용하지 않으면 혼돈이 없을거 같습니다.

[꿈꾸는가인]

꼬리말을 달고 보니 제가 너무 나선거 같네요. 글자가 150자로 제한되어 있다 보니... 제가 틀린 곳이 있다면 고수님들께서 수정해 주시길...

[김현수]

아..감사해요..그리고 전 90에서 사용하고 있구요..c, cc, ** ,모두 앞에 !를 같이 붙여서 comment로 사용하고 있어요..^^