from Jerome Huck
Good morning.
I did use Gfortran GCC 9.0 on Windows 7 64 bits on an old Fortran77. I
attached both the source code and the GCC outputs.
I have warning with COMMON. The COMMONS are the same along the code...
When compiled I have warnings with different sizes !!!
I do not know if it is a bug. I will say it is an issue.
There were some bugs/issues in the past see :
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=45044
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=45045
If you can give me some clues/hints on the subject.
Best regards.
Jerome Huck.
C PROGRAM No. 12: RECTANGULAR LIFTING SURFACE (VLM)
C -------------------------------------------------
C 3D-VLM CODE FOR SIMPLE WING PLANFORMS WITH GROUND EFFECT (BY JOE KATZ, 1974).
DIMENSION QF(6,14,3),QC(4,13,3),DS(4,13,4)
DIMENSION GAMA(4,13),DL(4,13),DD(4,13),DP(4,13)
DIMENSION A(52,52),GAMA1(52),DW(52),IP(52)
DIMENSION A1(5,13),DLY(13),GAMA1J(5),X(4)
COMMON/NO1/ DS,X,B,C,S,AR,SN1,CS1
COMMON/NO2/ IB,JB,CH,SIGN
COMMON/NO3/ A1
COMMON/NO4/ QF,QC,DXW
C
C ==========
C INPUT DATA
C ==========
C
IB=4
JB=13
X(1)=0.
X(2)=0.
X(3)=4.
X(4)=4.
B=13.
VT=1.0
ALPHA1=5.0
1000 continue
print*,' '
WRITE(6,101)
print*,' '
print*,' Input the x coordinates of the wing''s 4 corner points.'
print*,' '
print*,'X(1) is the root L.E., X(2) is the tip L.E., '
print*,'X(3) is the tip T.E., and X(4) is the root T.E.'
print*,' '
print*,' 2------3 '
print*,' | | '
print*,' | | '
print*,' | | '
print*,' | | '
print*,' 1------4 '
print*,' '
read*,x(1),x(2),x(3),x(4)
print*, ' Input the Wing Semi-Span, B (0 to stop)'
read*,b
print*,' Input alpha in degrees.'
read*,alpha1
if(b.eq.0.)stop
CH=1000.
C X(1) TO X(4) ARE X-COORDINATES OF THE WING'S FOUR CORNERPOINTS.
C B - WING SPAN, VT - FREE STREAM SPEED, B - WING SPAN,
C CH - HEIGHT ABOVE GROUND
C
C CONSTANTS
DXW=100.0*B
DO 1 I=1,IB
DO 1 J=1,JB
C GAMA(I,J)=1.0 IS REQUIRED FOR INFLUENCE MATRIX CALCULATIONS.
1 GAMA(I,J)=1.0
RO=1.
PAY=3.141592654
ALPHA=ALPHA1*PAY/180.
SN1=SIN(ALPHA)
CS1=COS(ALPHA)
IB1=IB+1
IB2=IB+2
JB1=JB+1
C
C =============
C WING GEOMETRY
C =============
C
CALL GRID
WRITE(6,101)
WRITE(6,102) ALPHA1,B,C,S,AR,VT,IB,JB,CH
C
C ========================
C AERODYNAMIC CALCULATIONS
C ========================
C
C INFLUENCE COEFFICIENTS CALCULATION
C
K=0
DO 14 I=1,IB
DO 14 J=1,JB
SIGN=0.0
K=K+1
CALL WING(QC(I,J,1),QC(I,J,2),QC(I,J,3),GAMA,U,V,W,1.0,I,J)
L=0
DO 10 I1=1,IB
DO 10 J1=1,JB
L=L+1
C A(K,L) - IS THE NORMAL VELOCITY COMPONENT DUE TO A UNIT VORTEX
C LATTICE.
10 A(K,L)=A1(I1,J1)
C ADD INFLUENCE OF WING'S OTHER HALF
CALL WING(QC(I,J,1),-QC(I,J,2),QC(I,J,3),GAMA,U,V,W,1.0,I,J)
L=0
DO 11 I1=1,IB
DO 11 J1=1,JB
L=L+1
11 A(K,L)=A(K,L)+A1(I1,J1)
IF(CH.GT.100.0) GOTO 12
C ADD INFLUENCE OF MIRROR IMAGE (DUE TO GROUND)
SIGN=10.0
CALL WING(QC(I,J,1),QC(I,J,2),-QC(I,J,3),GAMA,U,V,W,1.0,I,J)
L=0
DO 8 I1=1,IB
DO 8 J1=1,JB
L=L+1
8 A(K,L)=A(K,L)+A1(I1,J1)
C ADD MIRROR IMAGE INFLUENCE OF WING'S OTHER HALF.
CALL WING(QC(I,J,1),-QC(I,J,2),-QC(I,J,3),GAMA,U,V,W,1.0,I,J)
L=0
DO 9 I1=1,IB
DO 9 J1=1,JB
L=L+1
9 A(K,L)=A(K,L)+A1(I1,J1)
SIGN=0.0
12 CONTINUE
C
13 CONTINUE
C
C CALCULATE WING GEOMETRICAL DOWNWASH
C
UINF=VT
VINF=0.0
WINF=0.0
C THIS IS THE GENERAL FORMULATION FOR RIGHT HAND SIDE.
DW(K)=-(UINF*DS(I,J,1)+VINF*DS(I,J,2)+WINF*DS(I,J,3))
14 CONTINUE
C
C SOLUTION OF THE PROBLEM: DW(I)=A(I,J)*GAMA(I)
C
K1=IB*JB
DO 15 K=1,K1
15 GAMA1(K)=DW(K)
CALL DECOMP(K1,52,A,IP)
16 CONTINUE
CALL SOLVER(K1,52,A,GAMA1,IP)
C HERE * THE SAME ARRAY SIZE IS REQUIRED,
C AS SPECIFIED IN THE BEGINNING OF THE CODE
C
C WING VORTEX LATTICE LISTING
C
K=0
DO 17 I=1,IB
DO 17 J=1,JB
K=K+1
17 GAMA(I,J)=GAMA1(K)
C
C ==================
C FORCES CALCULATION
C ==================
C
FL=0.
FD=0.
FM=0.
QUE=0.5*RO*VT*VT
DO 20 J=1,JB
DLY(J)=0.
DO 20 I=1,IB
IF(I.EQ.1) GAMAIJ=GAMA(I,J)
IF(I.GT.1) GAMAIJ=GAMA(I,J)-GAMA(I-1,J)
DYM=QF(I,J+1,2)-QF(I,J,2)
DL(I,J)=RO*VT*GAMAIJ*DYM
C INDUCED DRAG CALCULATION
CALL WING(QC(I,J,1),QC(I,J,2),QC(I,J,3),GAMA,U1,V1,W1,0.0,I,J)
CALL WING(QC(I,J,1),-QC(I,J,2),QC(I,J,3),GAMA,U2,V2,W2,0.0,I,J)
IF(CH.GT.100.0) GOTO 194
CALL WING(QC(I,J,1),QC(I,J,2),-QC(I,J,3),GAMA,U3,V3,W3,0.0,I,J)
CALL WING(QC(I,J,1),-QC(I,J,2),-QC(I,J,3),GAMA,U4,V4,W4,0.0,I,J)
GOTO 195
194 W3=0.
W4=0.
195 WIND=W1+W2-W3-W4
C ADD INFLUENCE OF MIRROR IMAGE (GROUND).
ALFI=-WIND/VT
DD(I,J)=RO*DYM*VT*GAMAIJ*ALFI
C
DP(I,J)=DL(I,J)/DS(I,J,4)/QUE
DLY(J)=DLY(J)+DL(I,J)
FL=FL+DL(I,J)
FD=FD+DD(I,J)
FM=FM+DL(I,J)*(QF(I,J,1)-X(1))
20 CONTINUE
CL=FL/(QUE*S)
CD=FD/(QUE*S)
CM=FM/(QUE*S*C)
C
C OUTPUT
C
WRITE(6,104) CL,FL,CM,CD
WRITE(6,110)
DO 21 J=1,JB
DO 211 I=2,IB
211 GAMA1J(I)=GAMA(I,J)-GAMA(I-1,J)
DLYJ=DLY(J)/B*JB
21 WRITE(6,103) J,DLYJ,DP(1,J),DP(2,J),DP(3,J),DP(4,J),GAMA(1,J),
1GAMA1J(2),GAMA1J(3),GAMA1J(4)
C
C END OF PROGRAM
100 CONTINUE
C
C FORMATS
C
101 FORMAT(1H ,/,20X,'WING LIFT DISTRIBUTION CALCULATION (WITH GROUND
1 EFFECT)',/,20X,56('_'))
102 FORMAT(1H ,/,10X,'ALFA:',F10.2,8X,'B :',
1F10.2,8X,'C : ',F13.2,/,10X,
2'S :',F10.2,8X,'AR :',F10.2,8X,'V(INF) :',F10.2,/,10X,
3'IB :',I10,8X,'JB :',I10,8X,'L.E. HEIGHT:', F8.2,/)
103 FORMAT(1H ,I3,' I ',F9.3,' II ',4(F9.3,' I '),' I ',4(F9.3,' I '))
104 FORMAT(/,1H ,'CL=',F10.4,2X,'L=',F10.4,4X,'CM=',F10.4,3X,
1'CD=',F13.7)
110 FORMAT(1H ,/,5X,'I DL',4X,'II',22X,'DCP',22X,'I I',25X,
1'GAMA',/,118('='),/,5X,'I',15X,'I= 1',11X,'2',11X,'3',11X,
2'4',5X,'I I',5X,'1',11X,'2',11X,'3',11X,'4',/,118('='))
112 FORMAT(1H ,'QF(I=',I2,',J,X.Y.Z)= ',15(F6.1))
113 FORMAT(1H ,110('='))
C
GOTO 1000
END
C
SUBROUTINE GRID
DIMENSION QF(6,14,3),QC(4,13,3),DS(4,13,4),X(4)
COMMON/NO1/ DS,X,B,C,S,AR,SN1,CS1
COMMON/NO2/ IB,JB,CH,SIGN
COMMON/NO4/ QF,QC,DXW
C
PAY=3.141592654
C X(1) - IS ROOT L.E., X(2) TIP L.E., X(3) TIP T.E., AND X(4) IS ROOT T.E.
C IB: NO. OF CHORDWISE BOXES, JB: NO. OF SPANWISE BOXES
IB1=IB+1
IB2=IB+2
JB1=JB+1
C
C WING FIXED VORTICES LOCATION ( QF(I,J,(X,Y,Z))...)
C
DY=B/JB
DO 3 J=1,JB1
YLE=DY*(J-1)
XLE=X(1)+(X(2)-X(1))*YLE/B
XTE=X(4)+(X(3)-X(4))*YLE/B
C XLE AND XTE ARE L.E. AND T.E. X-COORDINATES
DX=(XTE-XLE)/IB
DO 1 I=1,IB1
QF(I,J,1)=(XLE+DX*(I-0.75))*CS1
QF(I,J,2)=YLE
C Note this vvv is different from the original Katz version
QF(I,J,3)=-QF(I,J,1)*SN1+CH
1 CONTINUE
C WAKE FAR FIELD POINTS
QF(IB2,J,1)=XTE+DXW
QF(IB2,J,2)=QF(IB1,J,2)
3 QF(IB2,J,3)=QF(IB1,J,3)
C
C WING COLOCATION POINTS
C
DO 4 J=1,JB
DO 4 I=1,IB
QC(I,J,1)=(QF(I,J,1)+QF(I,J+1,1)+QF(I+1,J+1,1)+QF(I+1,J,1))/4
QC(I,J,2)=(QF(I,J,2)+QF(I,J+1,2)+QF(I+1,J+1,2)+QF(I+1,J,2))/4
QC(I,J,3)=(QF(I,J,3)+QF(I,J+1,3)+QF(I+1,J+1,3)+QF(I+1,J,3))/4
C
C COMPUTATION OF NORMAL VECTORS
C
CALL PANEL(QF(I,J,1),QF(I,J,2),QF(I,J,3),QF(I+1,J,1),QF(I+1,J,2),
1QF(I+1,J,3),QF(I,J+1,1),QF(I,J+1,2),QF(I,J+1,3),QF(I+1,J+1,1),
2QF(I+1,J+1,2),QF(I+1,J+1,3),DS(I,J,1),DS(I,J,2),DS(I,J,3),
3DS(I,J,4))
4 CONTINUE
C
C B -IS SEMI SPAN, C -AV. CHORD, S - AREA
S=0.5*(X(3)-X(2)+X(4)-X(1))*B
C=S/B
AR=2.*B*B/S
C
RETURN
END
C
SUBROUTINE PANEL(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4,C1,C2,C3,S)
C CALCULATION OF PANEL AREA AND NORMAL VECTOR.
A1=X2-X3
A2=Y2-Y3
A3=Z2-Z3
B1=X4-X1
B2=Y4-Y1
B3=Z4-Z1
C NORMAL VECTOR
X=A2*B3-A3*B2
Y=B1*A3-A1*B3
Z=A1*B2-A2*B1
A=SQRT(X**2+Y**2+Z**2)
C1=X/A
C2=Y/A
C3=Z/A
C CALCULATION OF PANEL AREA
E1=X3-X1
E2=Y3-Y1
E3=Z3-Z1
F1=X2-X1
F2=Y2-Y1
F3=Z2-Z1
C NORMAL AREAS (F*B+B*E)
S11=F2*B3-F3*B2
S12=B1*F3-F1*B3
S13=F1*B2-F2*B1
S21=B2*E3-B3*E2
S22=E1*B3-B1*E3
S23=B1*E2-B2*E1
S=0.5*(SQRT(S11**2+S12**2+S13**2)+SQRT(S21**2+S22**2+S23**2))
RETURN
END
C
SUBROUTINE VORTEX(X,Y,Z,X1,Y1,Z1,X2,Y2,Z2,GAMA,U,V,W)
C SUBROUTINE VORTEX CALCULATES THE INDUCED VELOCITY (U,V,W) AT A POI
C (X,Y,Z) DUE TO A VORTEX ELEMENT VITH STRENGTH GAMA PER UNIT LENGTH
C POINTING TO THE DIRECTION (X2,Y2,Z2)-(X1,Y1,Z1).
PAY=3.141592654
RCUT=1.0E-10
C CALCULATION OF R1 X R2
R1R2X=(Y-Y1)*(Z-Z2)-(Z-Z1)*(Y-Y2)
R1R2Y=-((X-X1)*(Z-Z2)-(Z-Z1)*(X-X2))
R1R2Z=(X-X1)*(Y-Y2)-(Y-Y1)*(X-X2)
C CALCULATION OF (R1 X R2 )**2
SQUARE=R1R2X*R1R2X+R1R2Y*R1R2Y+R1R2Z*R1R2Z
C CALCULATION OF R0(R1/R(R1)-R2/R(R2))
R1=SQRT((X-X1)*(X-X1)+(Y-Y1)*(Y-Y1)+(Z-Z1)*(Z-Z1))
R2=SQRT((X-X2)*(X-X2)+(Y-Y2)*(Y-Y2)+(Z-Z2)*(Z-Z2))
IF((R1.LT.RCUT).OR.(R2.LT.RCUT).OR.(SQUARE.LT.RCUT)) GOTO 1
R0R1=(X2-X1)*(X-X1)+(Y2-Y1)*(Y-Y1)+(Z2-Z1)*(Z-Z1)
R0R2=(X2-X1)*(X-X2)+(Y2-Y1)*(Y-Y2)+(Z2-Z1)*(Z-Z2)
COEF=GAMA/(4.0*PAY*SQUARE)*(R0R1/R1-R0R2/R2)
U=R1R2X*COEF
V=R1R2Y*COEF
W=R1R2Z*COEF
GO TO 2
C WHEN POINT (X,Y,Z) LIES ON VORTEX ELEMENT; ITS INDUCED VELOCITY IS
1 U=0.
V=0.
W=0.
2 CONTINUE
RETURN
END
C
SUBROUTINE WING(X,Y,Z,GAMA,U,V,W,ONOFF,I1,J1)
DIMENSION GAMA(4,13),QF(6,14,3),A1(5,13)
DIMENSION DS(4,13,4)
COMMON/NO1/ DS
COMMON/NO2/ IB,JB,CH,SIGN
COMMON/NO3/ A1
COMMON/NO4/ QF
C
C CALCULATES INDUCED VELOCITY AT A POINT (X,Y,Z), DUE TO VORTICITY
C DISTRIBUTION GAMA(I,J), OF SEMI-CONFIGURATION - IN A WING FIXED
C COORDINATE SYSTEM.
U=0
V=0
W=0
IB1=IB+1
DO 1 I=1,IB1
DO 1 J=1,JB
C I3 IS WAKE VORTEX COUNTER
I3=I
IF(I.EQ.IB1) I3=IB
VORTIC=GAMA(I3,J)
IF(ONOFF.LT.0.1) GOTO 2
CALL VORTEX(X,Y,Z,QF(I,J,1),QF(I,J,2),QF(I,J,3),QF(I,J+1,1),QF(I,J
1+1,2),QF(I,J+1,3),VORTIC,U1,V1,W1)
CALL VORTEX(X,Y,Z,QF(I+1,J+1,1),QF(I+1,J+1,2),QF(I+1,J+1,3),
3QF(I+1,J,1),QF(I+1,J,2),QF(I+1,J,3),VORTIC,U3,V3,W3)
2 CALL VORTEX(X,Y,Z,QF(I,J+1,1),QF(I,J+1,2),QF(I,J+1,3),QF(I+1,J+1,1
2),QF(I+1,J+1,2),QF(I+1,J+1,3),VORTIC,U2,V2,W2)
CALL VORTEX(X,Y,Z,QF(I+1,J,1),QF(I+1,J,2),QF(I+1,J,3),QF(I,J,1),
4QF(I,J,2),QF(I,J,3),VORTIC,U4,V4,W4)
C
U0=U2+U4+(U1+U3)*ONOFF
V0=V2+V4+(V1+V3)*ONOFF
W0=W2+W4+(W1+W3)*ONOFF
A1(I,J)=U0*DS(I1,J1,1)+V0*DS(I1,J1,2)+W0*DS(I1,J1,3)
IF(SIGN.GE.1.0)
* A1(I,J)=U0*DS(I1,J1,1)+V0*DS(I1,J1,2)-W0*DS(I1,J1,3)
IF(I.EQ.IB1) A1(IB,J)=A1(IB,J)+A1(IB1,J)
U=U+U0
V=V+V0
W=W+W0
C
1 CONTINUE
RETURN
END
C
SUBROUTINE DECOMP(N,NDIM,A,IP)
REAL A(NDIM,NDIM),T
INTEGER IP(NDIM)
C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
C N = ORDER OF MATRIX. NDIM = DECLARED DIMENSION OF ARRAY A.
C A = MATRIX TO BE TRIANGULARIZED.
C IP(K) , K .LT. N = INDEX OF K-TH PIVOT ROW.
C
IP(N) = 1
DO 6 K = 1, N
IF(K.EQ.N) GOTO 5
KP1 = K + 1
M = K
DO 1 I = KP1, N
IF( ABS(A(I,K)).GT.ABS(A(M,K))) M=I
1 CONTINUE
IP(K) = M
IF(M.NE.K) IP(N) = -IP(N)
T = A(M,K)
A(M,K) = A(K,K)
A(K,K) = T
IF(T.EQ.0.E0) GO TO 5
DO 2 I = KP1, N
2 A(I,K) = -A(I,K)/T
DO 4 J = KP1, N
T = A(M,J)
A(M,J) = A(K,J)
A(K,J) = T
IF(T .EQ. 0.E0) GO TO 4
DO 3 I = KP1, N
3 A(I,J) = A(I,J) + A(I,K)*T
4 CONTINUE
5 IF(A(K,K) .EQ. 0.E0) IP(N) = 0
6 CONTINUE
RETURN
END
C
SUBROUTINE SOLVER(N,NDIM,A,B,IP)
REAL A(NDIM,NDIM), B(NDIM), T
INTEGER IP(NDIM)
C SOLUTION OF LINEAR SYSTEM, A*X = B.
C N = ORDER OF MATRIX.
C NDIM = DECLARED DIMENSION OF THE ARRAY A.
C B = RIGHT HAND SIDE VECTOR.
C IP = PIVOT VECTOR OBTAINED FROM SUBROUTINE DECOMP.
C B = SOLUTION VECTOR, X.
C
IF(N.EQ.1) GOTO 9
NM1 = N - 1
DO 7 K = 1, NM1
KP1 = K + 1
M = IP(K)
T = B(M)
B(M) = B(K)
B(K) = T
DO 7 I = KP1, N
7 B(I) = B(I) + A(I,K)*T
DO 8 KB = 1, NM1
KM1 = N - KB
K = KM1 + 1
B(K) = B(K)/A(K,K)
T = -B(K)
DO 8 I = 1, KM1
8 B(I) = B(I) + A(I,K)*T
9 B(1) = B(1)/A(1,1)
RETURN
ENDC:\Users\Jerome Huck\Desktop\VLM\vlm.f:235.17:
COMMON/NO1/ DS,X,B,C,S,AR,SN1,CS1
1
Warning: Named COMMON block 'no1' at (1) shall be of the same size as elsewhere (872 vs 832 bytes)
C:\Users\Jerome Huck\Desktop\VLM\vlm.f:237.17:
COMMON/NO4/ QF,QC,DXW
1
Warning: Named COMMON block 'no4' at (1) shall be of the same size as elsewhere (1636 vs 1008 bytes)
