1988 Fortran code

HelderC · ‎03-22-2023

Hi, everybody!

Learned Fortran language a long time ago, back in 1992, or so.

Now, back in the business, i got an 1988 fortran code that i need to put to work.

Manage to change some lines and put it to work under Visual Studio Code debug application.

Strangly, i can't compile to an executable, using any of the compilers i tried till now, resulting in some really strange errors (at least for me).

What can be the problem, and how can i sort it out?

Really desperate with this...

Arjen_Markus · ‎03-22-2023

Without the code or at the very least the error messages it will be impossible for us to help you. Can you show the code? There may be any number of things causing such messages and such a failure to build the program.

HelderC · ‎03-22-2023

Hi, Arjen!

I can send you the code, no problem with that, but it have 1135 lines... please check if it was properly attach to this message.

The error codes were different from compiler to compiler and all very strange, like the one regarding the start sentence "program", which was expected to have a digit on it!

Thank you in advance for your kind attention.

Helder

andrew_4619 · ‎03-22-2023

We I renamed your file console3.f90 and compiled it. There are some errors with FORMAT statements with lower case "L" where "1" seems to be the intention and a couple of spurious spaces e.g. "F6 . 1" rather than "F6.1" having corrected those it compiles and runs up to the point it asks me for some input.... The Modified version below. It does have a number of warnings and I had to turn standards checking off also.

program DOFRUB

implicit none

! *********************************************************
!
!   NON-LINEAR THREE DEGREE OF FREEDOM SYSTEM RESPONSE
!          USING FOURTH ORDER RUNGE-KUTTA METHOD
!       & BILINEAR VERTICAL & HORIZONTAL STIFFNESSES
!       WITH HORIZONTAL/VERTICAL ACELERATION INPUT AND
!                 DISPLACEMENT OUTPUT FILES
!        (INCLUDES WALE SHORE EFFECTS & HIGH BUILDUPS
!                   AND THE USE OF RUBBER CAPS)
! *********************************************************

integer NN,l,mm,n,hull,nsys,flag10,ll
integer flag1, flag2, flag3, flag4, flag5, flag6, flag7, flag8
integer KY1, KY2, KY3, KY4, KY5, WWW1, YYY1, UUU1, WWW2, YYY2, UUU2, WWW3, YYY3
integer UUU3, WWW4, YYY4, UUU4, UUU5, WWW5, YYY5, decrr
integer*2 i, j
real*8 beta, weight, h,Ik, gravity, AAA, ks, sidearea, keelarea, plside
real ac(2002), acv(2002), xx(2002), yy(2002), tt(2002), rrr(2002)
real t1,t2, t3, t4, t5, t6, t7, t8, x1, x2, x3, x4, x5, x6, x7, x8
real y1, y2, y3, y4, y5, y6, y7, y8, zeta, ff
real*8 m(4,4), cx(4,4), k(4,4), ko(4,4), crit2, crit3
real*8 baseside, basekeel, htside, htkeel
real*8 dtau, maxx, maxt, maxy, timex, timet
real*8 rf1, rf2, rf3, hf1, hf2, hf3, ampacc, mass, ampacmax
real*8 kvs, kvk, kvkp, khs, khk, kshp, kkhp, kvsp, base, counter, time
real*8 time1, time2, time3, time4, time5, time6, time7, time8
real*8 x, t, y, xold, told, yold, XSCL(6)
real*8 bbb, ccc, w12, w1, w22, w2, w32, w3, mode1, mode3
real*8 mmx1, mmang1, mmx3, mmang3, crit4, alpha, LLL
real*8 timey, mmmmm1, mmmmm2, mmmmm3, mmmmm4
real*8 R, S, TAU, A(6), B(6), C(6), D(6), E(6), F(6), G(6), HH(6)
real*8 br, amp, plkeel, u1, u2, XPRIM, VEL
real*8 KU1, KD1, khkb, QD1, XEL1, XMAX1, XMIN1, RR1, ZZ1, WZ1, VEL1
real*8 KU2, KD2, khsb, QD2, XEL2, XMAX2, XMIN2, RR2, ZZ2, WZ2, YPRIM1
real*8 KU3, KD3, kvsb1, QD3, YEL1, YMAX1, YMIN1, RR3, ZZ3, WZ3, DELTA
real*8 KU4, KD4, kvsb2, YEL2, YMAX2, YMIN2, RR4, ZZ4, WZ4, YPRIM2, VEL2
real*8 KU5, KD5, kvkb, QD4, YEL3, YMAX3, YMIN3, RR5, ZZ5, WZ5, YPRIM3
CHARACTER*40 DEC, DECV, quakname, hname, vname
character*40 sbfname, aclfname, outfname, vfname

!   COLLECT FOLLOWING VESSEL AND DOCK INFO: VESSEL WEIGHT, KG, I (about keel),
!   TIME INCREMENT OF DATA POINTS, VERTICAL STIFNESS OF SIDE AND
!   KEEL PIERS, HORIZONTAL STIFNESS OF SIDE AND KEEL PIERS,
!   GRAVITATIONAL CONSTANT, SIDE BLOCK BASE AND HEIGHT,
!   KEEL BLOCK BASE AND HEIGHT,
!   BLOCK-BLOCK AND BLOCK-HULL FRICTION COEFICIENTS,
!   SIDE AND KEEL BLOCK S PROPORTIONAL LIMIT,
!   SIDE PIER-VESSEL CONTACT AREA, KEEL PIER-VESSEL CONTACT AREA,
!   CAP BLOCK INCLINATION ANGLE

! OPEN INPUT FILES AND READ DATA

write (*,'(a)') ' ENTER SHIP/BUILDUP FILE NAME . . . '
read (*,'(a)') sbfname

open (4, file= sbfname, status='old', form='formatted')
read(4,*) weight, h, Ik, kvs, kvsp, kvk, AAA, Ks
read (4,*) khs, khk, kshp, kkhp, QD1, QD2, QD3, gravity
read(4,*) baseside, basekeel, htside, htkeel, u1, u2
read(4,*) br, plside, plkeel, sidearea, keelarea, zeta
read(4,*) hull, nsys, beta, QD4, kvkp
CLOSE (4)

write (*,*) 'DO YOU WANT RESPONSE OUTPUT FILES? (Y OR N)'
read (*, '(a)') dec
if (dec.eq.'Y'.or.dec.eq.'y') then
write (*,*) 'INPUT DESIRED RESISTANCE OUTPUT: (1,2,3,4,5)'
write (*,*) 'KEEL HORIZONTAL FORCE = 1'
write (*,*) 'SIDE BLOCK HORIZONTAL FORCE = 2'
write (*,*) 'LEFT SIDE BLOCK VERT FORCE = 3'
write (*,*) 'RIGHT SIDE BLOCK VERT FORCE = 4'
write (*,*) 'KEEL BLOCK VERTICAL FORCE = 5'
write (*,*) decrr
end if

do 12, i=1,3
do 13, j=1,3
m(i, j)=0.0
k(i, j)=0.0
cx(i, j)=0.0
ko(i, j)=0.0
13 continue
12 continue

!   CALCULATE SYSTEM PARAMETERS

mass=weight/gravity
LLL=sqrt((htside-htkeel)**2D0+(br/2D0)**2D0)
alpha=asin((htside-htkeel)/LLL)

M(1, 1) =mass
m(1, 3) =h*mass
m(2, 2) =mass
m(3, 1) =mass*h
m(3, 3)=Ik

k(1, 1)=(2D0*Ks+2D0*khs+khk)
k(1, 3)=(2D0*Ks*AAA+2D0*khs*LLL*sin(alpha))
k(3, 1)=k( 1,3)
k(2, 2)=(2D0*kvs-t-kvk)
k(3, 3)=(2D0*Ks*AAA**2D0+2D0*khs*((LLL*sin( alpha))**2D0)+(2D0*kvs*((LLL*cos(alpha))**2D0)-(weight*h)))
ko(1, 1)=k(1, 1)
ko(1, 3)=k(1, 3)
ko(3, 1)=k(3, 1)
ko(2, 2)=k(2, 2)
ko(3, 3)=k(3, 3)

!   DETERMINE NATURAL FREQUENCIES OF SYSTEM

bbb=-(m(1, 1)*k(3, 3)+m(3,3)*k(1, 1)-m(1,3)*k(3, 1)-m(3, 1)*k(1,3))/(m(1, 1)*m(3, 3)-m(l, 3)*m(3, 1))
ccc = (k(1, 1)*k(3,3)-k(1,3)*k(3,1))/(m(1, 1)*m(3,3)-m(1,3)*m(3, 1))

! NATURAL FREQ. MODE #1

w12=(-bbb-sqrt(bbb**2-4D0*ccc))/2D0
w1=sqrt(w12)

! NATURAL FREQ. MODE #2

w22=k(2, 2)/m(2, 2)
w2=sqrt(w22)
! NATURAL FREQ. MODE #3

w32=(-bbb+sqrt(bbb**2-4D0*ccc) )/2D0
w3=sqrt(w32)

! MODE SHAPE #1 & #3

mode1=(m(1 ,3)*w12-k(1, 3))/(-m(1, 1)*w12+k(1, 1))
mode3=(m(1, 3)*w32-k(1, 3))/(-m(1, 1)*w32+k( 1,1))
! DETERMINE C11 , C13 , C31 , C33
mmx1=m(1, 1)+m(1,3) /mode1
mmang1=mode1*m(3, 1)+m(3, 3)
mmx3=m(1, 1)+m( 1, 3)/mode3
mmang3=mode3*m(3, 1)+m(3, 3)
mmmmm1=2D0*zeta*mmx1*w1
mmmmm2=2D0*zeta*mmx3*w3
mmmmm3=2D0*zeta*mmang1 *w1
mmmmm4=2D0*zeta*mmang3*w3

cx(1 , 3)=(mmmmm1-mmmmm2)/(1/mode1-1/mode3)
cx(1 , 1)=mmmmm1-(cx(1 , 3)/mode1)
cx(2, 2) =2D0*zeta*m(2, 2)*w2
CX(3, 1) =(mmmmm3-mmmmm4)/( mode1-mode3)
cx(3, 3) =mmmmm3-(cx(3, 1)*mode1)

! READ IN ACCELERATION DATA

CALL ACCLINPT(amp, ac, acv, dtau, quakname, hname, vname)

! ESTABLISH FAILURE CRITERIA AND FLAGS

crit2=min (u1,u2)
crit3= (6.6D-1*baseside-1.2D1)/htside
crit4=basekeel/(6D0*htkeel)
ampacc=1D0
counter=0.0
ampacmax=0.0
10000 continue
write(*,*) ampacc
flag1=0
flag2=0
flag3=0
flag4=0
flag5=0
flag6=0
flag7=0
flag8=0
flag10=0
maxx=0.0
maxt=0.0
maxy=0.0
mm=0
x=0.0
y=0.0
t=0.0
xold=0.0
yold=0.0
told=0.0
R=0.0
S=0.0
TAU=0.0

! INITIALIZING BILINEAR VARIABLES

! INITIALIZING DELTA

if (kvs.eq.kvsp) then
YEL1=0.0
elseif (kvs.ne.kvsp) then
YEL1=(QD3/(kvs-kvsp))
endif
if (kvk.eq.kvkp) then
YEL3=0.0
elseif (kvk.ne.kvkp) then
YEL3=QD4/(kvk-kvkp)
endif
DELTA=weight/(2D0*kvs+kvk)
if (QD3.ge.0.0.or.QD4.ge.0.0) then
kvsb1=kvs
kvkb=kvk
goto 100
endif
if (DELTA.lt.YEL3.and.DELTA.lt.YEL1) then
kvsb1=kvs
kvkb=kvk
elseif (DELTA.ge.YEL3.or.DELTA.ge.YEL1) then
kvsb1=kvsp
kvkb=kvkp
DELTA=YEL3+(weight-(YEL3*(2D0*kvs+kvk)))/(2D0*kvsp+kvkp)
endif

100 continue

! INITIALIZING KEEL HORIZONTAL STIFFNESS

KU1=khk
KD1=kkhp
khkb=KU1
if (QD1.eq.0.0) goto 101
KY1=0
XEL1=QD1/(KU1-KD1)
XMAX1=0.0
XMIN1=0.0
RR1=0.0
ZZ1=0.0
WZ1=0.0
WWW1=0.0
YYY1=0.0
UUU1=0.0

101 continue

! INITIALIZING SIDE BLOCK HORIZONTAL STIFFNESS

KU2=khs
KD2=kshp
khsb=KU2
if (QD2.eq.0.0) goto 102
KY2=0
XEL2=QD2/(KU2-KD2)
XMAX2=0.0
XMIN2=0.0
RR2=0.0
ZZ2=0.0
WZ2=0.0
WWW2=0.0
YYY2=0.0
UUU2=0.0

102 continue

! INITIALIZING LEFT SIDE BLOCK VERTICAL STIFFNESS

KU3=kvs
KD3=kvsp
if (QD3 .eq. 0.0) goto 103
KY3=0
YMAX1=0.0
YMIN1=0.0
RR3=kvsb1*DELTA
ZZ3=0.0
WZ3=0.0
WWW3=0.0
YYY3=0.0
UUU3=0.0

103 continue

! INITIALIZING RIGHT SIDE BLOCK VERTICAL STIFFNESS

KU4=kvs
KD4=kvsp
kvsb2=kvsb1
if (QD3.eq.0.0) goto 104
KY4=0
YEL2=YEL1
YMAX2=0.0
YMIN2=0.0
RR4=kvsb2*DELTA
ZZ4=0.0
WZ4=0.0
WWW4=0.0
YYY4=0.0
UUU4=0.0

104 continue

! INITIALIZING KEEL VERTICAL STIFFNESS

KU5=kvk
KD5=kvkp
if (QD4.eq.0.0) goto 105
KY5=0
YMAX3=0.0
YMIN3=0.0
RR5=kvkb*DELTA
ZZ5=0.0
WZ5=0.0
WWW5=0.0
YYY5=0.0
UUU5=0.0

105 continue

! IMPLEMENTATION OF EQUATIONS OF MOTION INTO THE
! RUNGE-KUTTA FORMULAS

do 301, l=1, 2000

! CALCULATE BILINEAR STIFFNESS AND RESISTANCE

! CALCULATE KEEL HORIZONTAL BILINEAR STIFFNESS

if (QD1 .eq. 0.0) goto 106

CALL BILINALL(x, S, khkb, RR1, KD1, QD1, KU1, XEL1, XMAX1, XMIN1, KY1, ZZ1, WZ1, WWW1, YYY1, UUU1)

106 continue

! CALCULATE SIDE BLOCK HORIZONTAL BILINEAR STIFFNESS

XPRIM=+x+LLL*t*sin(alpha)

if (QD2 .eq. 0.0) goto 107

VEL=+S+LLL*TAU*sin( alpha)

CALL BILINALL(XPRIM,VEL,khsb,RR2,KD2,QD2,KU2,XEL2,XMAX2,XMIN2,KY2,ZZ2, WZ2, WWW2, YYY2,UUU2)

107 continue

! CALCULATE LEFT SIDE BLOCK VERTICAL BILINEAR STIFFNESS

YPRIM1=-y-t*LLL*cos(alpha) +DELTA

if (QD3 .eq. 0.0) goto 108
if (QD3 .gt. 0.0) then
VEL1=-R-TAU*LLL*cos(alpha)

CALL BILINALL(YPRIM1,VEL1,kvsb1,RR3,KD3,QD3,KU3,YEL1, YMAX1,YMIN1,KY3,ZZ3,WZ3,WWW3,YYY3,UUU3)

elseif (QD3 .lt. 0.0) then

CALL RUBBER(YPRIM1,kvsb1,RR3,KD3,QD3,KU3,YEL1)

endif

108 continue

! CALCULATE RIGHT SIDE BLOCK VERTICAL BILINEAR STIFFNESS

YPRIM2=-y+t*LLL*cos(alpha)+DELTA

if (QD3 .eq. 0.0) goto 109
if (QD3 .gt. 0.0) then

VEL2=-R+TAU*LLL*cos(alpha)

CALL BILINALL(YPRIM2, VEL2,kvsb2,RR4,KD4,QD3,KU4,YEL2,YMAX2,YMIN2,KY4,ZZ4,WZ4, WWW4, YYY4,UUU4)

elseif (QD3 .lt. 0.0) then

CALL RUBBER (YPRIM2,kvsb2,RR4,KD4,QD3,KU4,YEL2)

endif

109 continue

! CALCULATE KEEL VERTICAL STIFFNESS

YPRIM3=-y+DELTA

if (QD4 .eq. 0.0) goto 110
if (QD4 .gt. 0.0) then

CALL BILINALL(YPRIM3,-R,kvkb,RR5,KD5,QD4,KU5,YEL3,YMAX3,YMIN3,KY5,ZZ5,WZ5,WWW5,YYY5,UUU5)

elseif (QD4 .lt. 0.0) then

CALL RUBBER(YPRIM3,kvkb,RR5,KD5,QD4,KU5,YEL3)

endif

110 continue

! RECALCULATION OF DELTA

if (QD3.ge.0.0.or.QD4.ge.0.0) then
DELTA=weight/(2D0*kvs+kvk)
goto 120
endif
if (kvkb.eq.kvk) then
DELTA=weight/(2D0*kvs+kvk)
elseif (kvkb.gt.kvk) then
DELTA=YEL3 + (weight-(YEL3*(2D0*kvs+kvk)))/(2D0*kvsp-t-kvkp)
end if

120 continue

if (QD1.eq.0.0.and.QD2.eq.0.0.and.QD3.eq.0.0.and.QD4.eq.0.0) goto 111

! RECALCULATION OF STIFFNESS MATRIX VALUES

k(1, 1)=(2D0*Ks+2D0*khsb+khkb)
k(1,3)=(2D0*Ks*AAA+2D0*khsb*LLL*sin(alpha))
k(3, 1)=k(1, 3)
k(2, 2)=(kvsb1+kvsb2+kvkb)
k(3, 3) = (2D0*Ks*AAA**2D0+2D0*khsb*((LLL*sin(alpha))**2D0)+((kvsb1+kvsb2)*((LLL*cos(alpha))**2D0)-(weight*h)))

111 do 3000, ll=0,5
A(ll)=0.0
B(ll)=0.0
C(ll)=0.0
D(ll)=0.0
E(ll)=0.0
F(ll)=0.0
G(ll)=0.0
HH(ll)=0.0
3000 CONTINUE
mm=mm+l
do 302, NN=1,4
IF(NN.EQ.1) THEN
FF=0.0
ELSE IF (NN.EQ.2 .OR. NN.EQ.3) THEN
FF=5D-1
ELSE IF (NN.EQ.4) THEN
FF=1D0
ENDIF
A(NN)=dtau*(R+FF*D(NN-1))
B(NN)=dtau*(S+FF*E(NN-1))
C(NN)=dtau*(TAU+FF*F(NN-1))
D(NN)=dtau*((-cx(2,2)/m(2,2))*(R+FF*D(NN-1))-(k(2,2)/m(2,2))*(y+FF*A(NN-1))-amp*ampacc*acv(l)/2.54D0)
G(NN)=dtau*((-cx(1,1)/m(1,1))*(S+FF*E(NN-1))-(cx(1,3)/m(1,1))*(TAU+FF*F(NN-1))-(k(1,1)&
/m(1,1))*(x+FF*B(NN-1))-(k(1,3)/m(1,1))*(t+FF*C(NN-1))-ampacc*ac(1)/2.54D0)
HH(NN)=dtau*((-cx(3,3)/m(3,3))*(TAU+FF*F(NN-1))-(cx(3,1)/m(3,3))&
*(S+FF*E(NN-1))-(k(3,3)/m(3,3))*(t+FF*C(NN-1))+(m(3,1)/m(3,3))&
*((-cx(2,2)/m(2,2))*(R+FF*D(NN-1))-(k(2,2)/m(2,2))*(y+FF*A(NN-1)))*(t+FF*C(NN-1))&
-(k(3,1)/m(3,3))*(x+FF*B(NN-1))-(m(3,1)/m(3,3))*ampacc*ac(1)/2.54D0)

E(NN)=(m(1,1)*m(3,3)*G(NN)-m(1,3)*m(3,3)*HH(NN))/(m(3,3)*m(1,1)-m(1,3)*m(3,1))
F(NN)=(HH(NN)-(m(3,1)/m(3,3))*E(NN))
302 continue

! DETERMINING SYSTEM RESPONSE

yold=y
y=yold+(A(1)+2D0*A(2)+2D0*A(3)+A(4))/6D0

xold=x
x=xold+(B(1)+2D0*B(2)+2D0*B(3)+B(4))/6D0

told=t
t=told+(C(1)+2D0*C(2)+2D0*C(3)+C(4))/6D0

R=R+(D(1)+2D0*D(2)+2D0*D(3)+D(4) )/6D0

S=S+(E(1)+2D0*E(2)+2D0*E(3)+E(4))/6D0

TAU=TAU+(F(1)+2D0*F(2)+2D0*F(3)+F(4))/6D0

! MAXIMUM VALUES FOR TRANSLATIONS AND ROTATION

if (abs(xold) .gt. abs(maxx) ) then
timex=dtau*(l-1)
maxx=xold
end if
if (abs(told) .gt. abs(maxt) ) then
timet=dtau*(l-1)
maxt=told
endif
if (abs(yold).gt.abs(maxy)) then
timey=dtau*(l-1)
maxy=yold
endif

! CALCULATE VERTICAL AND HORIZONTAL FORCES CAUSED BY VESSEL,
! TEST FOR FAILURE

! CALCULATE FORCES ON SIDE/KEEL BLOCKS
if (QD3.eq.0.0) then
rf1=kvs*((weight/k(2,2))-yold-(LLL*cos(alpha) )*told)
rf2=kvs*( (weight/k(2, 2))-yold+( LLL*cos( alpha) )*told)
elseif (QD3.ne.0.0) then
rf1=RR3
rf2=RR4
endif

if (QD4.eq.0.0) then
rf3=kvk*((weight/k(2, 2))-yold)
elseif (QD4.ne.0.0) then
rf3=RR5
endif

if (QD2.eq.0.0) then
hf1=khs*(xold+LLL*told*sin(alpha))
hf2=khs*(xold+LLL*told*sin(alpha))
elseif (QD2.gt.0.0) then
hf1=RR2
hf2=RR2
endif

if (QD1.eq.0.0) then
hf3=khk*(xold)
elseif (QD1.gt.0.0) then
hf3=RR1
endif

! TEST FOR SIDE BLOCK SLIDING

if (flag1.eq. 1) then
go to 400
else if (hf1.lt.0.0.and.rf1.gt.0.0.and. u1*rf1+hf1+u2*rf1*cos(beta)*sin(beta)-rf1*cos(beta)*sin(beta).lt. 0.0) then
time1= dtau*(l-1)
flag1=1
else if (hf2.gt.0.0.and.rf2.gt.0.0.and. -u1*rf2+hf2-u2*rf2*(cos(beta) *sin(beta))&
+rf2*cos(beta)*sin(beta) .gt. 0.0) then
time1=dtau*(l-1)
flag1=1
endif
x1=xold
y1=yold
t1=told
400 continue

! TEST FOR KEEL BLOCK SLIDING

if (flag2.eq.1) then
go to 410
else if (rf3.gt.0.0.and.abs(hf3/rf3).gt.crit2) then
time2=dtau*(l-1)
flag2=1
endif
x2=xold
y2=yold
t2=told
410 continue
! TEST FOR SIDE BLOCK OVERTURNING

if (flag3.eq.1) then
go to 420
else if (hf1.lt.0.0.and. rf1.gt.0.0.and.abs(hf1/rf1).gt.crit3) then
time3= dtau*(l-1)
flag3=1
else if (hf2.gt.0.0.and.rf2.gt.0.0.and.abs(hf2/rf2 ).gt.crit3) then
time3=dtau*(l-1)
flag3=1
endif
x3=xold
y3=yold
t3=told
420 continue

! TEST FOR KEEL BLOCK OVERTURNING

if (flag4.eq.1) then
go to 430
else if (rf3.gt.0.0.and.abs(hf3/rf3).gt.crit4) then
time4=dtau*(l-1)
flag4=1
end if
x4=xold
y4=yold
t4=told
430 continue

! TEST FOR SIDE BLOCK LIFTOFF

if (flag5.eq.1) then
go to 440
else if (rf1.lt. 0.0 .or. rf2.lt. 0.0) then
time5=dtau*(t-1)
flag5=1
end if
x5=xold
y5=yold
t5=told
440 continue

! TEST FOR KEEL BLOCK LIFTOFF

if (flag6.eq.1) then
go to 450
else if (rf3.lt. 0.0) then
time6=dtau*(l-1)
flag6=1
end if
x6=xold
y6=yold
t6=told
450 continue

! TEST FOR SIDE BLOCK CRUSHING

if (flag7.eq.1) then
go to 460
else if (rf1.gt.0.0.and. (rf1/sidearea).gt.plside) then
flag7=1
time7=dtau*(l-1)

else if (rf2.gt.0.0.and.(rf2/sidearea).gt.plside) then
flag7=1
time7=dtau*(l-1)
end if
x7=xold
y7=yold
t7=told
460 continue

! TEST FOR KEEL BLOCK CRUSHING

if (flag8.eq.1) then
go to 470
else if (rf3.gt.0.0.and.(rf3/keelarea).gt.plkeel ) then
flag8=1
time8=dtau*(l-1)
end if
x8=xold
y8=yold
t8=told
470 continue

! CAPTURE OF DISPLACEMENT, ROTATION & RESISTANCE OUTPUT:

if (dec.ne.' Y'.and.dec.ne.'y') goto 301
xx(mm)=xold
tt(mm)=told
goto (501, 502, 503, 504, 505), decrr
501 if (QD1.eq.0.0) then
rrr(mm)=hf3
elseif (QD1.gt.0.0) then
rrr(mm)=RR1
end if
yy(mm)=yold
goto 506
502 if (QD2.eq.0.0) then
rrr(mm)=hf1
elseif (QD2.gt.0.0) then
rrr(mm)=RR2
xx(mm)=XPRIM
endif
yy(mm)=yold
goto 506
503 if (QD3.eq.0.0) then
rrr(mm)=rf1
elseif (QD3.ne.0.0) then
rrr(mm)=RR3
endif
yy(mm)=YPRIM1
goto 506
504 if (QD3.eq.0.0) then
683 rrr(mm)=rf2
684 elseif (QD3.ne.0.0) then
rrr(mm)=RR4
endif
yy(mm)=YPRIM2
goto 506
505 if (QD4.eq.0.0) then
rrr(mm)=rf3
elseif (QD4.ne.0.0) then
rrr(mm)=RR5
endif
yy( mm) =YPRIM3

506 continue

301 continue

go to 999

60000 continue
if (dec.ne.'Y'.and.dec.ne.'y') then
write(*,'(A)') ' I AM FINISHING. '
goto 20000
endif

! CREATION OF DISPLACEMENT, ROTATION, & RESISTANCE OUTPUT FILES:
CALL RESPALL(xx, yy, tt, rrr, dtau)

998 go to 20000

999 CONTINUE

if (ampacc.eq.1D0) then

write(*,'(a)') ' ENTER OUTPUT FILENAME ...'
read( *,'(a)' ) outfname
open (46, file=outfname, status='new', form='formatted')


write(46,4000) nsys
4000 format(1x, /,28x, '**** System ', I2, 1x, '****')
write(46, 4050) hull
4050 format( 1x,/, 30x, '** Hull ' , I3, 1x, '**')
write(46, 4100)
4100 format( 1x, //, 28x, ' * Ship Parameters *')
write(46, 4150)
4150 format( 1x,/,5x, 'Weight',8x, 'Moment of Inertia',9x, ' K.G.')
write(46, 4200) weight, Ik,h
4200 format( 1x, f9. 1, 1x, 'kips' , 1x, f11. 1, 1x, 'kips-in-sec2',3x, f6. 1, 1X, ' ins')
write(46, 4250)
4250 format( 1x, //, 26x, '* Drydock Parameters *')
write(46, 4300)
4300 format( 1x, /, 1x, 'Side Block Height', 3x, 'Side Block Width',3x, 'Keel Block Height', 3x, ' Keel Block Width')
write( 46, 4350) htside, baseside, htkeel , basekeel
4350 format (2x, f6.1 , 1x, ' ins' , 11x, f6.1, 1x, ' ins' , 11x, f6 . 1,1x, ' ins',9x, f6. 1, 1x, ' ins' )
write(46, 4400)
4400 format(1x, /, 1x, 'Side-to-Side Pier Distance' , 3x, 'Wale Shore Ht. ',3x,'Wale Shore Stiffness' , 2x, ' Cap Angle')
write(46, 4450) br, AAA, Ks, beta
4450 format (1x, t7, f6. 1, 1x, ' ins' , 17x, f6.1, 1x, ' ins' , 8x, f8. 1, 1x,'kips/in' , 1x, f5.3, 1x, ' rad')
write(46, 4470)
4470 format(1x, /,' ISide Side Pier Contact Area',3x, 'Total Keel Pier Contact Area' , 6X, ' kkhp')
write( 46, 4475 ) sidearea, keelarea, kkhp
4475 format ( 1x, 8x, f11. 1, 1x, ' in2' , 14x, f11. 1 , 1x, 'in2 ' , 10x, f7. 1 , 1x,'kips/in')
write(46, 45000)
45000 format( 1x,/, 1x, 'B/B Friction Coeff',3x,'H/B Friction Coeff' , 5x, 'kshp', 10x, 'kvsp')
write(46, 4550) u1 , u2, kshp, kvsp
4550 format (6x, f7. 3, 13x, f7.3,7x, f7.1 , 1x, ' kips/in' , 1x, f7.1, 1x,'kips/in')
write(46, 4600)
4600 format(1x, /, 1x, 'Side Pier Fail Stress Limit ',4x, 'Keel Pier',' Fail Stress Limit', 6x, 'kvkp')
write(46, 4650) plside, plkeel , kvkp
4650 format (1x, 10x, f7.3 , 1x, 'kips/in2' 15x, f7.3, 1x, 'kips/in2',6x, f7. 1, 1x, 'kips/in')
write(46, 4700)
4700 format (1x,/, 1x,'Side Pier Vertical Stiffness ', 3x, 'Side Pier',' Horizontal Stiffness')
write (46,4750) kvs, khs
4750 format (1x, 3x, f11.1,1x, 'kips/in', 11x, f11.1, 1x, 'kips/in')
write(46,4775)
4775 format (1x,/,1x, 'Keel Pier Vertical Stiffness', 3x,'Keel Pier Horizontal Stiffness')
write (46,4780) kvk, khk

4780 format (1x, 3x, f11.1, 1x, 'kips/in', 11x, f11.1, 1x, 'kips/in')
write(46,4782)
4782 format (1x,/,6x, 'QD1', 17x, 'QD2', 18x,'QD3', 17x, 'QD4')
write (46,4785) QD1, QD2, QD3, QD4
4785 format (2x, f8.1, 1x, 'kips', 7x, f8.1, 1x, 'kips', 8x, f8.1, 1x, 'kips',7x, f8.1, 1x, 'kips')
write (46,4800)
4800 format(1x,//, 20x,'* System Parameters and Inputs *')
write (46,4850) quakname
4850 format (1x,/,1x, 'Vertical acceleration input is ', A40)
write(46,4852) hname
4852 format (1x,/,1x, 'Horizontal acceleration input is ', A40)
write (46,4854)
4854 format (1x,/,1x, 'Vertical acceleration input is ', A40)
write (46,4875)
4875 format (1x, 20x, 'Earthquake Acceleration Time History.')

write(46,4995)
4995 format (1x,/,1x, 'Vertical/horizontal Ground Acceleration Ratio',3x, 'Data Time Increment')
write (46,4990) amp, dtau
4990 format (1x, 10x, f6.3, t55, f6.3, 1X, 'sec')
write (46, 4900)
4900 format (1x,/,1x, 'Gravitational constant', 3x,'% System Damping')
write (46,4950) gravity, zeta*100.
4950 format (1x,7x, f6.2,1x, 'in/sec2', 10x, f6.2, 1x, '%')
write (46, 5000)
5000 format(1x,/, 25x, 'Mass Matrix',/)
do 5100 i=1,3
write(46,5050) m(i,1), m(i,2), m(i,3)
5050 format (1x, f15.4, 5x, f15.4, 5x, f15.4)
5100 continue
write (46,5200)
5200 format (1x,/,25x, 'Damping Matrix',/)
do 5300 i=1,3
write (46,5250) cx(i,1), cx(i,2), cx(i,3)
5250 format (1x, f15.4, 5x, f15.4, 5x, f15.4)
5300 continue
write (46, 5400)
5400 format (1x,/, 25x, 'Stifness Matrix',/)
do 5500 i=1,3
write (46,5450) ko(i,1), ko(i,2), ko(i,3)
5450 format (1x, f15.4,5x, f15.4,5x, f15.4)
5500 continue
write (46, 5700)
5700 format (1x,//)
write (46,6000)
6000 format (1x, 'Undamped Natural Frequencies', t35, 'Mode #1', t50,'Mode #2', t65, 'Mode #3')
write (46, 6001) w1, w3, w2
6001 format (1x, t31, f7.3, 1x, 'rad/sec', t46, f7.3, 1x, 'rad/sec', t62, f7.3,' rad/sec')
WRITE(46, 6002)
6002 FORMAT (1X, 'Damped Natural Frequencies ', t35, ' Mode #l',t50,'Mode #2',t65,' Mode #3')
WRITE(46, 6500) w1*sqrt(1-zeta**2 ) , w3*sqrt(1-zeta**2 ),w2*sqrt(1-zeta**2)
6500 format (1x, t31, f7.3, 1x, 'rad/sec', t46, f7.3,1x, 'rad/sec', t62, f7.3,' rad/sec')
endif

write(46, 10500) ampacc*100, quakname
10500 format(1x, ///, 1x, 'For Earthquake Acceleration of ',f6.2,' % ', 'of the ', A40, /)

write(46, 25000)
25000 format(1x, 'Maximums/Failures' ,t26, 'X ( ins )', t36 , 'Y (ins)',t51,'Theta (rads)' ,t65 ,' Time (sec)')
write(46, 25001)
25001 format (1x, '----------------', t25, '--------', t35, '--------', t50,&
'-------------', t64, '-----------')
write (46,310) maxx,timex
310 format (1x, ' Maximum X' , t25, f9.6, t65, f5.2)
write (46,311) maxy,timey
311 format (1x, ' Maximum Y' , t35, f9.6, t65, f5.2)
write (46,312) maxt, timet
312 format (1x, ' Maximum Rotation', t50, f9.6, t65, f5.2)

if (flag1.eq.1) then
flag10=flag10+1
write (46,313) x1 ,y1 ,t1 , time1
313 format (1x, 'Side block sliding', t25, f9.6, t35, f9.6, t50, f9.6,t65, f5.2)

endif

if (flag2.eq.1) then
flag10=flag10+1
write (46,314) x2 , y2 , t2 , time2
314 format (1x, 'Keel block sliding' , t25, f9.6, t35 , f9.6, t50, f9.6,t65,f5.2)
endif

if (flag3.eq.1) then
flag10=flag10+1
write (46,315) x3,y3,t3,time3
315 format (1x, 'Side block overturning' , t25 , f 9.6, t35 , f 9.6, t50, f9.6,t65,f5.2)
endif

if (flag4.eq.1) then
flag10=flag10+1
write (46,316) x4, y4, t4, time4
316 format (1x, 'Keel block overturning' ,t25, f9.6, t35, f9.6, t50, f9.6,t65,f5.2)
endif

if (flag5.eq.1) then
flag10=flag10+1
write (46,317) x5 ,y5 , t5,time5
317 format (1x,'Side block liftoff ',t25 ,f9.6 , t35, f9.6, t50, f9.6,t65, f5.2)
end if

if (flag6.eq.1) then
flag10=flag10+1
write (46,318) x6, y6, t6, time6
318 format (1x, 'Keel block liftoff ', t25, f9.6, t35, f9.6, t50, f9.6,t65, f5.2)
end if
if (flag7.eq. 1) then
flag10=flag10+1
write (46,319) x7, y7, t7, time7
319 format (1x, 'Side block crushing' , t25,f9.6, t35, f9.6, t50, f9.6,t65, f5.2)
end if
if (flag8.eq.1) then
flag10=flag10+1
write (46,320) x8, y8, t8 ,time8
320 format (1x, 'Keel block crushing' , t25, f9.6, t35, f9.6, t50, f9.6,t65, f5.2)
end if

if (flag10.eq.0) then
write(46, 11000)
11000 format(1x,/, 1x, 'No failures occurred. ')
if (counter.eq.1.0 .and.flag10.eq.0) then
go to 60000
end if
if (counter.eq.0.0) then
ampacmax=ampacc
ampacc=ampacc+1D-1
counter=1.0
write(*,' ( A)')' In secondary looping stage. '
end if
end if
if (ampacc.le.ampacmax) go to 20000
if (counter.eq.1.0) then
ampacc = ampacc-1D-2
else if (counter.eq.0.0) then
ampacc=ampacc-1D-1
end if
go to 10000
20000 continue
stop
END PROGRAM dofrub


! ***************************************************
!SUBROUTINE WHICH PROMPTS FOR AND READS IN HORIZONTAL
!AND VERTICAL ACCELERATION TIME HISTORY FILES
!AND THE TIME STEP AND EARTHQUAKE NAME
!****************************************************

SUBROUTINE ACCLINPT (amp, ac, acv, dtau, quakname, hname, vname)
integer N
real ac(2002), acv(2002)
real*8 amp, dtau, dtauh, dtauv
character*40 aclfname, vfname, decv, quakname, hname, vname
character*40 hquaknam, vquaknam

! READ IN ACCELERATION DATA

! HORIZONTAL ACCELERATION

700 write(*,'(a)') ' ENTER HORIZONTAL ACCELERATION FILE NAME...'
read (*, '(a)') aclfname
open(44, file=aclfname, status='old', form='formatted')
write(*,'(a)') 'READING HORIZONTAL ACCELERATION FILE...'
read(44, '(a)') hquaknam
read(44, '(a)') hname
read(44,'(f9.4)') dtauh
do 300, n=1,2000
read(44,*) ac(n)
300 continue

! VERTICAL ACCELERATION

307 write(*,'(a)') 'WILL YOU USE A VERTICAL ACCELERATION FILE? '
write(*, '(a)') ' (Y/N) '
read(*, '(a)') decv
if (decv.eq.'Y') then
    write(*, '(a)') ' ENTER VERTICAL ACCELERATION FILE NAME...'
    read(*, '(a)')vfname
    open(45,file=vfname,status='old', form='formatted')
    write(*,'(a)') 'READING VERTICAL ACCELERATION FILE... '
    amp=1.0
read(45, '(a)') vquaknam
read(45, '(a)') vname
read(45, '(f9.4)') dtauv
if(dtauv.ne.dtauv.or.vquaknam.ne.hquaknam) then
    write(*,'(a)') ' INCOMPATIBLE ACCELERATION FILES!!! '
    write(*,'(a)') 'REINPUT COMPATIBLE FILES '
    goto 700
endif
do 305,n=1,2000
read (45,*) acv(n)
305 continue
endif

if (decv.eq.'N') then
    do 306,n=1,2000
        acv(n)=ac(n)
    306 continue
    write(*, '(a)') 'INPUT DESIRED VERTICAL/HORIZONTAL ACCELERATION RATIO '
    read(*,*) amp
endif

if(decv.ne. 'Y'.and.decv.ne. 'N') then
    write(*, '(a)') 'TRY AGAIN '
    goto 307
endif

quakname=hquaknam
dtau=dtauh
 CLOSE (44)
CLOSE(45)

RETURN
END

! **********************************************************
! SUBTOUTINE WHICH CALCULATES THE BILINEAR HORIZONTAL OR
! VERTICAL STIFFNESS AND RESISTANCE
! **********************************************************

SUBROUTINE BILINALL (U, V, PK, RR,KD,QD,KU,UEL, UMAX, UMIN, KY, ZZ, WZ, WWW, YYY, UUU)

real*8 U,V, RR, KD,QD, KU, UEL, PK
real*8 UMAX, UMIN, ZZ, WZ
integer WWW, YYY,UUU, KY

! BEGGINING OF BILINEAR LOGIC

!CHECK IF RESPONSE STILL ON INITIAL ELASTIC LINE

if (KY.lt.0) goto 4040
if (KY.gt.0) goto 3480
     RR=KU*U
    PK=KU

!CHECK IF RESPONSE HAS GONE PLASTIC

if (U.gt.-UEL.and.U.lt.UEL) goto 4720

!RESPONSE IS NOW PLASTIC

if(U.lt.-UEL) goto 4040

!RESPONSE IS ON THE TOP PLASTIC LINE

3220 KY=1
    PK=KD
    RR=KD*U+QD
    WWW=0
    YYY=0
    ZZ=0.0
    goto 4720

! CHECKS IF VELOCITY SHIFTS FROM POSITIVE TO NEGATIVE

3480 if(V.gt.0) goto 3720

!CHECKS IF ON THE RIGHT ELASTIC LINE

if(YYY.gt.0) goto 3630

! CALCULATE VALUE OF UMAX

    ZZ=U2
3630 YYY=1
    UMAX=ZZ

! CHECK IF RESPONSE SHIFTS TO LOWER PLASTIC LINE

3720 if(U.lt.(UMAX-2*UEL)) goto 4040

! CHECK IF RESPONSE SHIFTS TO TOP PLASTIC LINE

if(U.gt.UMAX) goto 3220

!CHECKS IF RESPONSE RETURNS TO TOP PLASTIC LINE

if(YYY.eq.0) goto 3220

! RESPONSE IS ON THE RIGHT ELASTIC LINE

KY=1
PK=KU
RR=KU*U+(KD-KU)*UMAX+QD
goto 4720

! CHECKS IF VELOCITY SHIFTS TO POSITIVE

4040 if(V.gt.0) goto 4350

! CHECKS IF RESPONSE REMAINS ELASTIC

if(WWW.eq.1) goto 4350

! RESPONSE IS ON THE BOTTOM PLASTIC LINE

4150 KY=-1
    PK=KD
    RR=KD*U-QD
    UUU=0
    WZ=0.0
    goto 4720

! CHECK IF RESPONSE IS ON THE LEFT ELASTIC LINE

4350 if(UUU.gt.0) goto 4370
    WZ=U
4370 UUU=1
    UMIN=WZ

! CHECK IF RESPONSE RETURNING TO TOP PLASTIC LINE

if(U.gt.(UMIN+2*UEL)) goto 3220

! CHECK IF RESPONSE RETURNS TO BOTTOM PLASTIC LINE

if (U.lt.UMIN) goto 4150

! RESPONSE IS ON THE LEFT ELASTIC LINE

WWW=1
RR=KU*U+(KD-KU)+UMIN-QD
PK=KU

4720 continue
RETURN
END

! ***************************************************
!SUBROUTINE WHICH CALCULATRES THE RUBBER CAP VERTICAL
! STIFNESS AND RESISTANCE
! ***************************************************

SUBROUTINE RUBBER (U,PK, RR, KD, QD, KU, UEL)

real*8 U, RR, KD, QD, KU, UEL, PK

! BEGGINING OF RUBBER LOGIC

! CHECKS IF RESPONSE STILL ON INITIAL ELASTIC LINE

if(U.gt.UEL) goto 3220
    RR=KU*U
    PK=KU
    goto 4720

! RESPONSE IS ON THE SECONF ELASTIC LINE

3220 continue
    PK=KD
    RR=KD*U+QD

4720 continue
RETURN
END

! *********************************************************
! SUBROUTINE WHICH CREATES VERTICAL, ROTATIONAL, HORIZONTAL
! DISPLACEMENT AND DESIGNATED RESISTANCE OUTPUT FILES
! *********************************************************

SUBROUTINE RESPALL (xx, yy, tt, rrr, dtau)

real xx(2002), tt(2002), yy(2002), rrr(2002)
real*8 dtau, time
character*40 xname, yname, tname, rrname
integer N

! CREATION OF DISPLACEMENT & ROTATION OUTPUT FILES

write(*, '(a)') 'ENTER X OUTPUT FILE NAME...'
read(*, '(a)') xname
open(47, file=xname, status='new', form='formatted')

write(*, '(a)') 'ENTER Y OUTPUT FILE NAME...'
read(*, '(a)') yname
open(48, file=yname, status='new', form='formatted')

write(*, '(a)') 'ENTER THETA OUTPUT FILE NAME...'
read(*, '(a)') tname
open(49, file=tname, status='new', form='formatted')

write(*, '(a)') 'ENTER RESISTANCE OUTPUT FILE NAME...'
read(*, '(a)') rrname
open(41, file=rrname, status='new', form='formatted')

do 308,n=1,2000
time=dtau*(n-1)
write(47,7000) time,xx(n)
7000 format(f7.3,10x,e13.6)

write(48,7010) time,yy(n)
7010 format(f7.3, 10x,e13.6)

write(49,7020) time,tt(n)
7020 format(f7.3,10x,e13.6)
write(41,7030) time, rrr(n)
7030 format(f7.3,10x,e13.6)

308 continue
RETURN
END

HelderC · ‎03-22-2023

Thank you so much Andrew!

Had some issues with those before, too, but thought they were solved when i first saw the programm asking for the first imputs.

Thank you so much!

I'll work from here! Let's see what i'll get with the compiler.

Do you advise some compiler in particular?

Helder

Arjen_Markus · ‎03-22-2023

These extra spaces in formats might be due to the conversion from fixed form to free form. In fixed form spaces are only significant inside string literals.

andrew_4619 · ‎03-22-2023

I just used the latest OneApi compiler, it is free to download and use.