 # Numerical resolution of coupled partial differential equations governing natural flow of fluid in a rectangular cavity

Hello Dear all, I create this topic to have a help about Numerical Simulations of coupled partial differential equations governing natural convection in a rectangular cavity filled with fluid and heated from below. I need one exemple of Fortran code using finite volume method or Lattice Boltzmann method or Differential Quadrature method or any other method. In the following picture, I presented the Equations (1) to (9) that I want to solve.

Hi Justin, is this a homework problem? Please tag it as homework if so.

What does your current solution look like? What did you try and where did you get stuck? If you provide as much detail as possible, people here will be able to help with exactly the problem that you have.

This is a problem that I have been instructed to solve in order to familiarize myself with the numerical solution of similar equations. I did the discretization by the differential quadrature method but I have difficulties writing a Fortran code.

Did you write, compile, and run any Fortran code before? If not, Learn Fortran - Fortran Programming Language should get you going with the basics. Give it a try and please provide information about a specific thing that you get stuck with.

The code I wrote for the sample case when Ha = 0 is the following, but when I compile it by taking:

• N=21
• M=21
• NTT=1
• H=0.001
• EPS1=0.000001
• EPS2=0.000001
• Pr=0.71
• Ra=1000
• ALAY=0.66

the error is very large and I can’t get the exact results. It would be difficult to explain the discretization method and the various parameters used as well as the acronyms, which is why I am asking if anyone has any other sample code of any method to provide me with detail.

``````       Program pdqnc

implicit real*8(A-H,O-Z)

dimension AX(50,50),BX(50,50),XX(50),AY(50,50),BY(50,50),YY(50)
dimension U(50,50),V(50,50),PD(401,401),P(50,50)
dimension T(50,50),W(50,50),DT(50,50),DW(50,50)

! Imput number of mesh points, maximum iteration number, time step
! size, convergence criteria, prandtl number and rayleigh number

write(6,*) 'N,M,NTT,H,EPS1,EPS2,PR,RA'

! Mesh generation and set-up of initial flow field

call GRID(N,M,XX,YY,U,V,W,T,DT,DW)

! Computation of weighting coefficient of the first and second order
! derivative in the x and y direction.

call WCOEF(1,N,XX,AX,BX)
call WCOEF(1,M,YY,AY,BY)

! Form the matrix of equation system 6.21. The obtained matrix is
! then decomposed by LU decomposition (see equation 5.27).

call PDCO(N,M,AX,BX,AY,BY,PD)

NT=1

! Start iteration

88 continue

! Implement boundary conditions for vorticity and temperature

call BOCO(N,M,U,V,W,T,AX,AY)

!4-Stage Runge-Kutta method to solve resultant ordinary diffrencial
! equationqs for vortivity and temperature

call RUNGU(N,M,H,RA,PR,AX,BX,AY,BY,U,V,T,W,DT,DW)

! Compute stream function

call COREP(N,M,W,P,PD)

!Update stream function on the boundary

call BPCO(N,M,AX,AY,P)

!Compute velocity U,V

call CORUV(N,M,AX,AY,U,V,P)

!Find out maximum residuals for equations 6.22 and 6.23

ERR1=DT(1,1)**2
ERR2=DW(1,1)**2
do 91 I=1,N
do 91 J=1,M
ERR1=DT(I,J)**2+ERR1
ERR2=DW(I,J)**2+ERR2

91 continue
NM=N*M
ERR1=sqrt(ERR1/NM)
ERR2=sqrt(ERR2/NM)
write(6,113) NT,ERR1,ERR2
113 format(1X,I10,5X,'  ERR1=  ',E10.5,5X,'  ERR2=  ',E10.5)

! Check conergence criteria

if (ERR1.le.EPS1.and.ERR2.le.EPS2) go to 92
if (NT.ge.NTT) go to 92
NT=NT+1
go to 88
92 continue

! Compute nusselt numbers

call NUSEC(N,M,XX,YY,AX,U,T)

! Interpolate numerical results from coarse mesh to fine mesh for
! plotting and other purpose

call OUTPLT(N,M,XX,YY,U,V,W,T,P)

stop

end program pdqnc

! This subroutine is to compute the PDQ weighting coefficients of
! the first and second order derivatives when the coordinates of
! mesh points are input

subroutine WCOEF(N0,NN,XX,AX,BX)

implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),XX(50)
dimension AIX(50)

NC=NN-N0+1

! Compute the weighting coefficients of the first order derivatives

do 72 I=N0,NN
II=I-N0+1
AX1=1.

do 73 J=N0,NN
if (J.eq.I) go to 73
AX1=AX1*(XX(I)-XX(J))

73 continue
72 AIX(II)=AX1

do 74 I=1,NC
II=I+N0-1
do 74 J=1,NC
JJ=J+N0-1

if (J.eq.I) go to 74
AX(I,J)=AIX(I)/((XX(II)-XX(JJ))*AIX(J))
74 continue

do 77 I=1,NC
AXI=0.

do 78 J=1,NC
if (J.eq.I) go to 78
AXI=AXI+AX(I,J)
78 continue
AX(I,I)=-AXI
77 continue

! Compute the weighting of the second order derivate

do 300 I=1,NC
II=I+N0-1
do 300 J=1,NC
JJ=J+N0-1

if (J.eq.I) go to 300
BX(I,J)=2.*AX(I,J)*(AX(I,I)-1./(XX(II)-XX(JJ)))
300 continue

do 400 I=1,NC
BXI=0.

do 410 J=1,NC
if(J.eq.I) go to 410
BXI=BXI+BX(I,J)
410 continue
BX(I,I)=-BXI
400 continue

return
end

! This subroutine is to generate the coordinates of mesh points

subroutine GRID(N,M,XX,YY,U,V,W,T,DT,DW)

implicit real*8(A-H,O-Z)
dimension XX(50),YY(50),U(50,50),V(50,50)
dimension T(50,50),W(50,50),DT(50,50),DW(50,50)

write(6,*) 'ALAX,ALAY'

PI=4.*ATAN(1.)
!  Lx=1.
! Ly=1.
do 71 I=1,N
ALFA=PI*FLOAT(I-1)/FLOAT(N-1)
! ALFA1=(cos(PI/(2.*N))-cos((((2.*I)-1.)*PI)/(2.*N)))
! ALFA2=(cos(PI/(2.*N))-cos((((2.*N)-1.)*PI)/(2.*N)))
71 XX(I)=0.5*(1.-cos(ALFA))
! 71 XX(I)=(ALFA1/ALFA2)*Lx
X11=XX(1)

! Stretch the mesh towards the boundary, ALAX.lt.1

do 76 I=1,N
76 XX(I)=(1.-ALAX)*(3.*PI*XX(I)**2-2.*XX(I)**3)/PI**2+ALAX*XX(I)

do 81 J=1,M

ALFA=PI*FLOAT(J-1)/FLOAT(M-1)
! ALFA1=(cos(PI/(2.*M))-cos((((2.*J)-1.)*PI)/(2.*M)))
!  ALFA2=(cos(PI/(2.*M))-cos((((2.*M)-1.)*PI)/(2.*M)))
81 YY(J)=0.5*(1.-cos(ALFA))
! 81  YY(J)=(ALFA1/ALFA2)*Ly
! Stretch the mesh towards the boundary, ALAY.lt.1

do 86 J=1,M

86 YY(J)=(1.-ALAY)*(3.*PI*YY(J)**2-2.*YY(J)**3)/PI**2+ALAY*YY(J)
YY(MH)=0.5

! Set up initial flow field

do 18 J=1,M
do 10 I=1,N

U(I,J)=0.
V(I,J)=0.
T(I,J)=0.
W(I,J)=0.
DT(I,J)=0.
DW(I,J)=0.

10 continue
T(1,J)=1.0

18 continue

return
end

! Implement boudary conditions for vorticity and temperature

subroutine BOCO(N,M,U,V,W,T,AX,AY)
implicit real*8(A-H,O-Z)
dimension AX(50,50),AY(50,50),U(50,50),V(50,50)
dimension T(50,50),W(50,50)

! Update vorticity on the Bottom and top boundaries

do 72 I=1,N

SUM1=0.
SUM2=0.

do 10 J=1,M
SUM1=SUM1+AY(1,J)*U(I,J)
10 SUM2=SUM2+AY(M,J)*U(I,J)
W(I,1)=SUM1
72 W(I,M)=SUM2

! Update vorticity on the left and right boundaries

do 73 J=1,M
SUM3=0.
SUM4=0.

do 20 I=1,N
SUM3=SUM3+AX(1,I)*V(I,J)
20 SUM4=SUM4+AX(N,I)*V(I,J)
W(1,J)=-SUM3
W(N,J)=-SUM4

! Temperature condition at left and right boundaries

T(1,J)=1.
T(N,J)=0.

73 continue
return
end

! 4-Runge-Kutta method to solve ordinary differential equations
! for vorticity and temperature. The scheme is shown in equation 5.15

subroutine RUNGU(N,M,H,RA,PR,AX,BX,AY,BY,U,V,T,W,DT,DW)

implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),AY(50,50),BY(50,50)
dimension U(50,50),V(50,50),T(50,50),W(50,50)
dimension CT(50,50),CW(50,50),DT(50,50),DW(50,50)

! Stage 1

call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,T,W,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,T,DT)

do 3 I=1,N
do 3 J=1,M
CW(I,J)=W(I,J)+0.25*H*DW(I,J)
3 CT(I,J)=T(I,J)+0.25*H*DT(I,J)

! Stage 2

call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,CT,CW,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,CT,DT)

do 5 I=1,N
do 5 J=1,M
CW(I,J)=W(I,J)+H*DW(I,J)/3.
5 CT(I,J)=T(I,J)+H*DT(I,J)/3.

! Stage 3

call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,CT,CW,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,CT,DT)

do 7 I=1,N
do 7 J=1,M
CW(I,J)=W(I,J)+H*DW(I,J)/2.
7 CT(I,J)=T(I,J)+H*DT(I,J)/2.

! Stage 4

call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,CT,CW,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,CT,DT)

do 9 I=2,N-1
do 9 J=2,M-1
W(I,J)=W(I,J)+H*DW(I,J)
9 T(I,J)=T(I,J)+H*DT(I,J)

call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,T,W,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,T,DT)

return
end

! Compute the residual of vorticity equation 6.22

subroutine TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,CT,CW,DW)

implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),AY(50,50),BY(50,50)
dimension U(50,50),V(50,50),CT(50,50),CW(50,50),DW(50,50)

do 10 I=2,N-1
do 10 J=2,M-1
DWI=0.
DTI=0.
DDWI=0.

do 20 K=1,N
DDWI=DDWI+BX(I,K)*CW(K,J)
DTI=DTI+AX(I,K)*CT(K,J)
20 DWI=DWI+AX(I,K)*CW(K,J)
DWJ=0.
DDWJ=0.

do 30 K=1,M
DWJ=DWJ+AY(J,K)*CW(I,K)
30 DDWJ=DDWJ+BY(J,K)*CW(I,K)
DW(I,J)=(DDWI+DDWJ)*PR-U(I,J)*DWI-V(I,J)*DWJ-DTI*PR*RA

10 continue

return
end

! Compute the residual of temperature equation 6.23

subroutine TRFT(N,M,AX,BX,AY,BY,U,V,CT,DT)

implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),AY(50,50),BY(50,50)
dimension U(50,50),V(50,50),CT(50,50),DT(50,50)

! Update the temperature on the top and bottom boundaries using
! equations 6.42b and 6.42c

do 70 I=2,N-1
C2=0.
C3=0.

do 80 J=2,M-1
C2=C2+(AY(1,J)*AY(M,M)-AY(M,J)*AY(1,M))*CT(I,J)
C3=C3+(AY(1,J)*AY(M,1)-AY(M,J)*AY(1,1))*CT(I,J)
80 continue
CT(I,1)=-C2/(AY(1,1)*AY(M,M)-AY(M,1)*AY(1,M))

70 CT(I,M)=-C3/(AY(1,M)*AY(M,1)-AY(M,M)*AY(1,1))

! Compute residuals of temperature equation 6.23 for all interior
! points

do 10 I=2,N-1
do 10 J=2,M-1
DTI=0.
DDTI=0.

do 20 K=1,N
DDTI=DDTI+BX(I,K)*CT(K,J)
20 DTI=DTI+AX(I,K)*CT(K,J)
DTJ=0.
DDTJ=0.

do 30 K=1,M
DTJ= DTJ+AY(J,K)*CT(I,K)
30 DDTJ=DDTJ+BY(J,K)*CT(I,K)
DT(I,J)=DDTI+DDTJ-U(I,J)*DTI-V(I,J)*DTJ

10 continue
return
end

! Form the matrix of equation system 6.21. the two layer approach
! is used to implement the boundary condition for stream function.
! The obtained matrix is decomposed by LU decomposition

subroutine PDCO(N,M,AX,BX,AY,BY,PD)

implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),AY(50,50),BY(50,50)
dimension PD(401,401),PL(401,401),PU(401,401)
NN=N-4
MM=M-4
NM=NN*MM
N1=N-1
M1=M-1

! AXN and AYM are used in equation 6.30 and 6.31 for two-layer
! approach

AXN=AX(1,2)*AX(N,N1)-AX(N,2)*AX(1,N1)
AYM=AY(1,2)*AY(M,M1)-AY(M,2)*AY(1,M1)

do 10 I=3,N-2
do 10 J=3,M-2
IJ=MM*(I-3)+J-2

do 20 K=1,NM
20 PD(IJ,K)=0.

do 30 K1=3,N-2
KJ=MM*(K1-3)+J-2

! Equation 6.30 is directly substituted into equation 6.21

AXK1=AX(1,K1)*AX(N,N1)-AX(N,K1)*AX(1,N1)
AXKN=AX(N,K1)*AX(1,2)-AX(1,K1)*AX(N,2)

30 PD(IJ,KJ)=BX(I,K1)-(AXK1*BX(I,2)+AXKN*BX(I,N1))/AXN

do 40 K2=3,M-2
IK=MM*(I-3)+K2-2

! Equation 6.31 is directly substituted into equation 6.21

AYK1=AY(1,K2)*AY(M,M1)-AY(M,K2)*AY(1,M1)
AYKM=AY(M,K2)*AY(1,2)-AY(1,K2)*AY(M,2)

40 PD(IJ,IK)=PD(IJ,IK)+BY(J,K2)-(AYK1*BY(J,2)+AYKM*BY(J,M1))/AYM

10 continue

!LU decomposition using equation 5.27

do 50 I=1,NM
do 50 J=1,NM

if(J.ge.I) go to 51
J1=J-1
SUM1=0.

if(J1.eq.0) go to 57

do 55 JJ=1,J1
55 SUM1=SUM1+PL(I,JJ)*PU(JJ,J)

57 PL(I,J)=(PD(I,J)-SUM1)/PU(J,J)
go to 50

51 if(I.eq.1) go to 52
I1=I-1
SUM2=0.

do 53 II=1,I1
53 SUM2=SUM2+PL(I,II)*PU(II,J)
PU(I,J)=PD(I,J)-SUM2
go to 50

52 PU(I,J)=PD(I,J)
50 continue

! The tringular matrices [PL] and [PU] are stored in [PD]

do 60 I=1,NM
do 60 J=1,NM

if(J.ge.I) go to 61
PD(I,J)=PL(I,J)
go to 60

61 PD(I,J)=PU(I,J)
60 continue

return
end

! The algebric equation system for stream function is solved by LU
! decomposition

subroutine COREP(N,M,W,P,PD)

implicit real*8(A-H,O-Z)
dimension PD(401,401),P(50,50),W(50,50)
dimension B(401),BI(401),BJ(401)
NN=N-4
MM=M-4
NM=NN*MM

do 10 I=3,N-2
do 10 J=3,M-2
IJ=MM*(I-3)+J-2
B(IJ)=W(I,J)
10 continue

! Equation 5.26a is used to get the solution at the first step

do 40 I=1,NM
if(I.eq.1) go to 41
I1=I-1
SUM1=0.

do 42 J=1,I1
42 SUM1=SUM1+PD(I,J)*BI(J)
BI(I)=B(I)-SUM1
go to 40
41 BI(1)=B(1)
40 continue

! Equation 5.26b is used to get the solution at the second step

do 50 I=1,NM
IV=NM-I+1

if(I.eq.1) GO TO 51
I2=I-1
SUM2=0.

do 52 J=1,I2
IVJ=NM-J+1
52 SUM2=SUM2+PD(IV,IVJ)*BJ(IVJ)
BJ(IV)=(BI(IV)-SUM2)/PD(IV,IV)

go to 50

51 BJ(IV)=BI(IV)/PD(IV,IV)

50 continue

do 60 I=3,N-2
do 60 J=3,M-2
IJ=MM*(I-3)+J-2
60 P(I,J)=BJ(IJ)

return
end

! Implement boundary condition for stream function. The two-layer
! approach is used

subroutine BPCO(N,M,AX,AY,P)

implicit real*8(A-H,O-Z)
dimension AX(50,50),AY(50,50),P(50,50)
N1=N-1
M1=M-1

! AXN and AYM are used in equations 6.30 and 6.31 for two-layer
!approach

AXN=AX(1,2)*AX(N,N1)-AX(N,2)*AX(1,N1)
AYM=AY(1,2)*AY(M,M1)-AY(M,2)*AY(1,M1)

! Equation 6.31 is used to update the stream function

do 20 I=3,N-2
SUM3=0.
SUM4=0.

do 21 J=3,M-2
SUM3=SUM3-(AY(1,J)*AY(M,M1)-AY(M,J)*AY(1,M1))*P(I,J)
21 SUM4=SUM4-(AY(M,J)*AY(1,2)-AY(1,J)*AY(M,2))*P(I,J)
P(I,2)=SUM3/AYM
20 P(I,M1)=SUM4/AYM

! Equation 6.30 is used to update the stream function

do 10 J=2,M-1
SUM1=0.
SUM2=0.

do 11 I=3,N-2
SUM1=SUM1-(AX(1,I)*AX(N,N1)-AX(N,I)*AX(1,N1))*P(I,J)

11 SUM2=SUM2-(AX(N,I)*AX(1,2)-AX(1,I)*AX(N,2))*P(I,J)
P(2,J)=SUM1/AXN
10 P(N1,J)=SUM2/AXN

do 30 J=1,M
P(1,J)=0.
30 P(N,J)=0.

do 40 I=2,N1
P(I,1)=0.
40 P(I,M)=0.

return
end

! Compute velocity U,V from stream fonction

subroutine CORUV(N,M,AX,AY,U,V,P)

implicit real*8(A-H,O-Z)
dimension AX(50,50),AY(50,50),U(50,50),V(50,50),P(50,50)

do 10 I=2,N-1
do 10 J=2,M-1
DPI=0.

! Compute V velocity for interior points

DO 20 k=1,N
20 DPI=DPI+AX(I,K)*P(K,J)
V(I,J)=-DPI
DPJ=0.

! Compute U velocity for interior points

do 30 k=1,M
30 DPJ=DPJ+AY(J,K)*P(I,K)
10 U(I,J)=DPJ

do 40 I=1,N
U(I,1)=0.
V(I,1)=0.
U(I,M)=0.
40 V(I,M)=0.

do 50 J=2,M-1
U(1,J)=0.
V(1,J)=0.
U(N,J)=0.
50 V(N,J)=0.

return
end

! Compute the Nusselt numbers

subroutine NUSEC(N,M,XX,YY,AX,U,T)

implicit real*8(A-H,O-Z)
dimension XX(50),YY(50),X2(201),Y2(201),AX(50,50),U(50,50)
dimension T(50,50),QQ(50,50),QI(50),QT(50,201),CNU(50),CNUT(201)

! Compute the local heat flux in a horizontal direction Q(X,Y)

do 10 I=1,N
do 10 J=1,M
DTI=0.

do 11 K=1,N

11 DTI=DTI+AX(I,K)*T(K,J)
10 QQ(I,J)=U(I,J)*T(I,J)-DTI
! 10 QQ(I,J)=DTI
! Set up a fine mesh

NN=101
MM=101
TDX=(XX(N)-XX(1))/FLOAT(NN-1)
TDY=(YY(M)-YY(1))/FLOAT(MM-1)

do 20 I=1,NN
20 X2(I)=FLOAT(I-1)*TDX

do 30 J=1,MM
30 Y2(J)=FLOAT(J-1)*TDY

! Interpolate Q(X,Y) to the fine mesh in the Y direction

do 40 I=1,N
do 41 J=1,M
QI(J)=QQ(I,J)

41 continue

do 42 J=1,MM
YJ=Y2(J)

call CZ(M,YJ,YY,QI,QOUT)

QT(I,J)=QOUT

42 continue
40 continue

! Compute average Nusselt number on the vertical line

do 188 I=1,N
SUM=0.0

do 198 J=2,MM
SUM=SUM+0.5*(QT(I,J)+QT(I,J-1))*TDY

198 continue
CNU(I)=SUM

188 continue

do 52 I=1,NN
XI=X2(I)

call CZ(N,XI,XX,CNU,CNUOUT)

CNUT(I)=CNUOUT

52 continue

! Compute average Nusselt number throughout the cavity

CNUAVG=0.0

do 918 I=2,NN
CNUAVG=CNUAVG+0.5*(CNUT(I)+CNUT(I-1))*TDX
918 continue

! Find out the maximum and minimum local Nusselt numbers at X=0
! together with their positions

AMIN=QT(1,1)
AMAX=QT(1,1)

do 518 J=1,MM
if(QT(1,J).le.AMAX) go to 888
AMAX=QT(1,J)
YMAX=Y2(J)

888 if (QT(1,J).ge.AMIN) go to 518
AMIN=QT(1,J)
YMAX=Y2(J)

518 continue
NH=(N+1)/2

! OUTPUT THE Nusselt Numbers

write(6,70) CNU(1)
70 FORMAT(2X,'AVERAGE NUSSELT NUMBER AT X=0',5X,F15.4)
write(6,72) CNU(NH)
72 FORMAT(2X,'AVERAGE NUSSELT NUMBER AT X=0.5',3X,F15.4)
write(6,74) CNUAVG
74 FORMAT(2X,'AVERAGE NUSSELT NUMBER OVER CAVITY',F15.4)
write(6,76) AMAX,YMAX
76 FORMAT(2X,'NUMAX AND POSITION AT X=0',5X,F10.4,5X,F10.4)
write(6,78) AMIN,YMIN
78 FORMAT(2X,'NUMIN AND POSITION AT X=0',5X,F10.4,5X,F10.4)

return
end

! Lagrange interpolation

subroutine CZ(N,XX,X,F,FV)

implicit real*8(A-H,O-Z)
dimension X(50),F(50)

FV=0.

do 10 I=1,N

do 20 J=1,N
if(J.eq.I) go to 21

go to 22

20 continue

10 continue

return
end

! Interpolation the numerical results from a coarse mesh to a fine
!  mesh

subroutine OUTPLT(N,M,X1,Y1,U,V,W,G,P)

implicit real*8(A-H,O-Z)
dimension X1(50),Y1(50),X2(101),Y2(101),P(50,50),W(50,50)
dimension U(50,50),V(50,50),G(50,50),UT(101,101),VT(101,101)
dimension GT(101,101),PT(101,101),WT(101,101),UI(51),UJ(51)
dimension VI(51),VJ(51),WI(51),WJ(51),PI(51),PJ(51),GI(50),GJ(50)

X2(1)=X1(1)
Y2(1)=Y1(1)

write(6,*) 'Read the values of NN, MM for interpolation'

TDX=(X1(N)-X1(1))/FLOAT(NN-1)
TDY=(Y1(M)-Y1(1))/FLOAT(MM-1)

! Generate fine mesh

do 20 I=2,NN
20 X2(I)=X2(I-1)+TDX

do 30 J=2,MM
30 Y2(J)=Y2(J-1)+TDY

! Interpolate in the X direction

do 40 J=1,M
do 41 I=1,N

UI(I)=U(I,J)
VI(I)=V(I,J)
WI(I)=W(I,J)
GI(I)=G(I,J)
PI(I)=P(I,J)

41 continue

do 42 I=1,NN
XI=X2(I)

call CZ(N,XI,X1,UI,UOUT)
call CZ(N,XI,X1,VI,VOUT)
call CZ(N,XI,X1,WI,WOUT)
call CZ(N,XI,X1,GI,GOUT)
call CZ(N,XI,X1,PI,POUT)

UT(I,J)=UOUT
VT(I,J)=VOUT
WT(I,J)=WOUT
GT(I,J)=GOUT
PT(I,J)=POUT

42 continue
40 continue

! Interpolate in the Y direction

do 50 I=1,NN
do 51 J=1,M

UJ(J)=UT(I,J)
VJ(J)=VT(I,J)
WJ(J)=WT(I,J)
GJ(J)=GT(I,J)
PJ(J)=PT(I,J)

51 continue

do 52 J=1,MM
YJ=Y2(J)

call CZ(M,YJ,Y1,UJ,UOUT)
call CZ(M,YJ,Y1,VJ,VOUT)
call CZ(M,YJ,Y1,WJ,WOUT)
call CZ(M,YJ,Y1,GJ,GOUT)
call CZ(M,YJ,Y1,PJ,POUT)

UT(I,J)=UOUT
VT(I,J)=VOUT
WT(I,J)=WOUT
GT(I,J)=GOUT
PT(I,J)=POUT

52 continue
50 continue

! Output data for contour plotting by tecplot (at software)

open(7,file='NSTO.dat')
!  Write(7,*) 'TITLE= "  " '
!  write(7,*) 'VARIABLE=X,Y,U,V,T,P,W'
!  write(7,*) 'ZONE I=',NN,'J=',MM,',F=POINT'

do 850 J=1,MM
do 850 I=1,NN

write(7,950) X2(I),Y2(J),UT(I,J),VT(I,J),GT(I,J),PT(I,J),WT(I,J)

850 continue

950 FORMAT(F7.4,1X,F7.4,1X,E11.5,2X,E11.5,2X,E11.5,2X,E11.5,2X,E11.5)

close(7)
return
end
``````

My generic advice is to print out intermediate results of your program to figure out where things are going wrong.

It is better to declare variables explicitly and use IMPLICIT NONE instead of writing

`implicit real*8(A-H,O-Z)`

The code

``````dimension AX(50,50),BX(50,50),XX(50),AY(50,50),BY(50,50),YY(50)
dimension U(50,50),V(50,50),PD(401,401),P(50,50)
dimension T(50,50),W(50,50),DT(50,50),DW(50,50)
``````

should be rewritten by defining a PARAMETER, say `ngrid` equal to 50 and using it in the declarations. Is the 401 related to the 50 in some way? Does the code work if the discretization is made finer by increasing the value of ngrid?

``````write(6,*) 'N,M,NTT,H,EPS1,EPS2,PR,RA'
``````

can you check that the values are read properly by printing them out? Reading many parameters in a single statement is error-prone.

1 Like

Thanks for your comments. I will use it to review the code and come back. But I can’t check that the values are read properly by printing them out.
The code don’t work if the discretization is made finer by increasing the value of ngrid.

You MUST use explicit typing if you want assistance from me. There is an undefined variable reference in the code you posted. Luckily I have a tool that helps me find it. It is so easy to make a typo and use a new name (whose value is undefined) when you meant to use an existing name with a defined value. With implicit typing, how are you going to find it?

1 Like

I agree with @themos about adding IMPLICIT NONE. Compilng the code with gfortran -Wall, I get

``````  202 |        YY(MH)=0.5
|
Warning: 'mh' is used uninitialized [-Wuninitialized]
``````

If a user sees many compiler warnings that can be disregarded, she may miss the important ones. Before the message above, gfortran emits numerous messages such as

``````  797 |         dimension GT(101,101),PT(101,101),WT(101,101),UI(51),UJ(51)
|                    1
Warning: Array 'gt' at (1) is larger than limit set by '-fmax-stack-var-size=', moved from stack to static storage. This makes the procedure unsafe when called recursively, or concurrently from multiple threads. Consider using '-frecursive', or increase the '-fmax-stack-var-size=' limit, or change the code to use an ALLOCATABLE array. [-Wsurprising]
xcavity.f90:797:32:
``````

Compiling without the -Wall option, there is no warning about uninitialized mh and the only messages are ones such as

``````   62 |        do 91 J=1,M
|                  1
Warning: Fortran 2018 deleted feature: Shared DO termination label 91 at (1)
xcavity.f90:118:18:``````

And the first advice ("use `implicit none`") is free. For a Fortran beginner, I bet that’s tremendous help already.