[siddiquesil]

!===============================================
This is a general purpose transient 2-d CFD code (FORTRAN 90)based
on FVM.
!************************************************* ************
Note to the users:
Please go through this brief instruction before attempting anything with this code
!************************************************* ************
1. Most of the variables used in the code are listed in the file "list_of_var". However, some new variables are also used, typically for the USER portion. Those variables are explained just before opening the input file in subroutine GRID

2. The coordinate system used is right-handed. The systems are as:
mode=1 : cartesian; mode=2: cylindrical(axisym); mode=3: polar

3. The user is reqd. to use only six subroutines almost. They are:

a) Subroutine GRID: For specification of input data and initialization of arrays. It calls subroutine UGRID if one selects uniform gridding;
and calls subroutine MESH if one selects non-uniform gridding. The
selection is made by giving a value to a variable called NGRID introduced
in subroutine GRID. Do not forget to specify following in subroutine grid:
# mode: to specify coordinate system. For a system involving radial
direction, also specify r(1).

# lsolve: to select which variables you want to solve; eg. lsolve(4)=.true.
will make eqn corresponding to nf=4(ie., energy eqn) to be solved.
Note that nf=1,2,3,4,np,6 correspond to solving for u,v,pc,T,p,conc. in the code.

!# ngrid: to specify choice of grid. For non-uniform grid, subroutine MESH is
called. It takes input from a file "grid.dat", where proportional lengths of
non-uniform grids are specified. However, for class-problems; generally
uniform grid will suffice, for which specify NGRID =1
# no. of grids in coordinate dirns, domain size, time step(dt) etc.
# eru,erv,ert: For convergence criteria
# property values as reqd.
# initialise all arrays, logical variables

b) Subroutine START: For giving initial condns for variables (say u,v,T)
for the problem

c) Subroutine DENSE : For specifying density of the working medium.

d) Subroutine BOUND: For giving boundary condns. for the problem

e) Subroutine GAMSOR: For giving diffusion coeff and source term for
all nfs (i.e., all eqns). As per symbols used in Patankar's book;
the symbol 'ap' used in the code corresponds to 'Sp' and the symbol
'con' used in the code corresponds to 'Sc'. 'gam(:,: )' corresponds to
the diffusion coefficient 'gamma' as per general 'phi' formulation

f) Subroutine Printout: This subroutine will give the plot files which contain isotherms and velocity vectors plot compatible with matlab or tecplot (as you choose). There is a logical switch called ltecplot declared in the programme. If you declare ltecplot as true, then you will get matlab files only. You can run the ".m" files generated in MATLAB to see the outputs.However, you are encouraged to develop output format according to your own choice and needs.

4. For steady state soln; use a typically large value of time step(dt)
as per instruction of the instructor

5. The code deals with a sample problem of natural convection in a rectangular cavity.The cavity is heated from two vertical sides by employing a uniform heat
flux (symbolised as FLUX in the code). The top and bottom are insulated. Due to heating,
there is a natural convection driven flow inside the liquid that develops with time. We are interested in studying the evolution of temperature field and velocity field inside the cavity.

6.Once you get the code, first compile it as it is. For compilation of the code in UNIX Systems, either compile in IBM m/c using the compilation command: xlf90 code.f OR complile in a COMPAQ machine using the compilation command f90 -fast code.90

Note: for compiling and running in IBM, you have to rename the file from code.f90 to code.f
Once it is compiled, execute the file a.out by typing a.out in the command prompt of your terminal. You will soon see that in your screen certain figures are appearing. We are basically printing the values of u,v,T and their relative errors in every iteration to monitor convergence.
Once it converges within a timestep, it goes to the next timestep. In this programme, I have kept timestep, i.e., dt=10s. The time domain I have kept as 510s, and in subroutine printout
I have stated that the MATLAB output should be written every 50 timestep. Accordingly, after 500s, this programme will open a file out500.m You can run the file in MATLAB to see various plots. Note that the streamfunction is symbolised as pc in the output. Do not confuse it with the pressure correction that we are solving to get the fluid flow.The figure 1 will give you temperature contours, figure 2 will give you velocity vectors and figure 3 will give you streamfunctions.Once you start solving problems on your own, initially you will be solving steady state problems only. For that, use dt=1.0e30 (a typically large no.) and tlast=1.0e30. Change the name
!of the .m file you want to have by altering the filename in subroutine printout. Note that
your tlast (i.e. end of time domain) should be same as dt for a steady state problem

[siddiquesil]

!################################################# #################
!================================================= =================
!GLOBAL DECLARATIONS
!================================================= =================
MODULE GLOBAL
IMPLICIT NONE

INTEGER,PARAMETER :: ni=40
INTEGER,PARAMETER :: nj=40
INTEGER,PARAMETER :: nij=40
INTEGER,PARAMETER :: nfmax=10
INTEGER,PARAMETER :: np=11
INTEGER,PARAMETER :: nfx3=nfmax+1

REAL,SAVE,DIMENSION(NI,NJ) :: u,v,pc,t,p,conc,rho,gam,con,ur,vr
REAL,SAVE,DIMENSION(NI,NJ) :: aip,aim,ajp,ajm,ap,ap0
REAL,SAVE,DIMENSION(NI,NJ) :: du,dv
REAL,SAVE,DIMENSION(NI,NJ) :: tO,uO,vO,told,uold,vold,clold,uhat,vhat
REAL,SAVE,DIMENSION(NI,NJ) :: fle,flw,fln,fls

REAL,SAVE,DIMENSION(NI,NJ,6) :: cofu,cofv,cofp,cof

REAL,SAVE,DIMENSION(NI) :: x,xu,xdif,xcv,xcvs,fv,fvp,fx,fxm,xcvi,xcvip,ara
REAL,SAVE,DIMENSION(NJ) :: y,yv,ydif,ycv,ycvs,ycvr,ycvrs,arx,arxj,arxjp
REAL,SAVE,DIMENSION(NJ) :: r,rmn,sx,sxmn,fy,fym

REAL,SAVE,DIMENSION(NIJ) :: pt,qt

REAL,SAVE,DIMENSION(NFX3) :: relax
REAL,SAVE,DIMENSION(NFX3) :: ntimes

REAL,SAVE :: rhocon,ert,eru,erv,time,flux
REAL,SAVE :: xl,yl,dt,tlast,dt2,time2,dt3,time3,dt4,time4,dt5,t ime5
REAL,SAVE :: smax,ssum,flow,diff,acof
REAL,SAVE :: ts,tl,g,beta,tref,rhoref,amiu,constant1,constant2, constant3,tk,cp
REAL,SAVE :: ax1,ay1,rhoinitial, tinitial

REAL,SAVE :: deltmx,delumx,delvmx

INTEGER,SAVE :: itll,NF,NRHO,NGAM,L1,L2,L3,M1,M2,M3,ngrid,ny1
INTEGER,SAVE :: IST,JST,ITER,ISTP,IPREF,JPREF,MODE,ITERL,NTIME

CHARACTER(LEN=4),DIMENSION(NFX3) :: title

LOGICAL,SAVE,DIMENSION(NFX3) :: LSOLVE,LPRINT,LBLK
LOGICAL,SAVE :: LSTOP,lconv,ltecplot
END MODULE GLOBAL
!================================================= =================

[siddiquesil]

!================================================= =====================
!***************** PROGRAM MAIN ***********************************
!================================================= =====================
PROGRAM MAIN
USE GLOBAL
IMPLICIT NONE

CALL GRID
CALL geomet
CALL START
mainloop: DO WHILE(.NOT.lstop)
CALL DENSE
CALL BOUND
CALL OLDVAL
CALL COEFF
END DO mainloop
WRITE(*,*)" PROGRAM FINISHED. "
END PROGRAM MAIN

!================================================= =====================

[siddiquesil]

!===============================================
!this subroutine calculates discretization eqn coeffs as per power law
!************************************************* ************
subroutine diflow
use global
implicit none
real :: Temp

acof=diff
if(flow.eq.0.) return
temp=diff-abs(flow)*0.1
acof=0.
If(temp.le.0.) return
temp=temp/diff
acof=diff*temp**5
end subroutine diflow

!===============================================

[siddiquesil]

!================================================= =====================
!THIS SUBROUTINE SOLVES DISCRETISATION EQUATIONS BY 'LINE BY LINE TDMA'
!************************************************* *********************
SUBROUTINE SOLVE
USE GLOBAL
IMPLICIT NONE
INTEGER :: I,J,II,JJ,N,NT,ISTF,JSTF,IT1,IT2,JT1,JT2
REAL :: BL,BLP,BLM,BLC,TEMP,DENOM
REAL,DIMENSION(NI,NJ) ::f

SELECT CASE(nf)
CASE(1)
f = u
CASE(2)
f = v
CASE(3)
f = pc
CASE(4)
f = t
CASE(np)
f = p
CASE(6)
f =conc
END SELECT

!----------------------------------------------------------------------
ISTF=IST-1
JSTF=JST-1
IT1=L2+IST
IT2=L3+IST
JT1=M2+JST
JT2=M3+JST
!----------------------------------------------------------------------
DO 999 NT=1,NTIMES(NF)
N=NF
!----------------------------------------------------------------------
IF(.NOT.LBLK(NF)) GO TO 10
PT(ISTF)=0.
QT(ISTF)=0.
DO 11 I=IST,L2
BL=0.
BLP=0.
BLM=0.
BLC=0.
DO 12 J=JST,M2
BL=BL+AP(I,J)
IF(J.NE.M2) BL=BL-AJP(I,J)
IF(J.NE.JST) BL=BL-AJM(I,J)
BLP=BLP+AIP(I,J)
BLM=BLM+AIM(I,J)
BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J)+AIM(I,J)*F(I-1,J) &
& +AJP(I,J)*F(I,J+1)+AJM(I,J)*F(I,J-1)-AP(I,J)*F(I,J)
12 END DO
DENOM=BL-PT(I-1)*BLM
IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E30
PT(I)=BLP/DENOM
QT(I)=(BLC+BLM*QT(I-1))/DENOM
11 END DO
BL=0.
DO 13 II=IST,L2
I=IT1-II
BL=BL*PT(I)+QT(I)
DO 13 J=JST,M2
13 F(I,J)=F(I,J)+BL
PT(JSTF)=0.
QT(JSTF)=0.
DO 21 J=JST,M2
BL=0.
BLP=0.
BLM=0.
BLC=0.
DO 22 I=IST,L2
BL=BL+AP(I,J)
IF(I.NE.L2) BL=BL-AIP(I,J)
IF(I.NE.IST) BL=BL-AIM(I,J)
BLP=BLP+AJP(I,J)
BLM=BLM+AJM(I,J)
BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J)+AIM(I,J)*F(I-1,J) &
& +AJP(I,J)*F(I,J+1)+AJM(I,J)*F(I,J-1)-AP(I,J)*F(I,J)
22 END DO
DENOM=BL-PT(J-1)*BLM
IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E30
PT(J)=BLP/DENOM
QT(J)=(BLC+BLM*QT(J-1))/DENOM
21 END DO
BL=0.
DO 23 JJ=JST,M2
J=JT1-JJ
BL=BL*PT(J)+QT(J)
DO 23 I=IST,L2
23 F(I,J)=F(I,J)+BL
10 CONTINUE
!----------------------------------------------------------------------
DO 90 J=JST,M2
PT(ISTF)=0.
QT(ISTF)=F(ISTF,J)
DO 70 I=IST,L2
DENOM=AP(I,J)-PT(I-1)*AIM(I,J)
PT(I)=AIP(I,J)/DENOM
TEMP=CON(I,J)+AJP(I,J)*F(I,J+1)+AJM(I,J)*F(I,J-1)
QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM
70 END DO
DO 80 II=IST,L2
I=IT1-II
80 F(I,J)=F(I+1,J)*PT(I)+QT(I)
90 END DO
!----------------------------------------------------------------------
DO 190 JJ=JST,M3
J=JT2-JJ
PT(ISTF)=0.
QT(ISTF)=F(ISTF,J)
DO 170 I=IST,L2
DENOM=AP(I,J)-PT(I-1)*AIM(I,J)
PT(I)=AIP(I,J)/DENOM
TEMP=CON(I,J)+AJP(I,J)*F(I,J+1)+AJM(I,J)*F(I,J-1)
QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM
170 END DO
DO 180 II=IST,L2
I=IT1-II
180 F(I,J)=F(I+1,J)*PT(I)+QT(I)
190 END DO
!----------------------------------------------------------------------
DO 290 I=IST,L2
PT(JSTF)=0.
QT(JSTF)=F(I,JSTF)
DO 270 J=JST,M2
DENOM=AP(I,J)-PT(J-1)*AJM(I,J)
PT(J)=AJP(I,J)/DENOM
TEMP=CON(I,J)+AIP(I,J)*F(I+1,J)+AIM(I,J)*F(I-1,J)
QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM
270 END DO
DO 280 JJ=JST,M2
J=JT1-JJ
280 F(I,J)=F(I,J+1)*PT(J)+QT(J)
290 END DO
!----------------------------------------------------------------------
DO 390 II=IST,L3
I=IT2-II
PT(JSTF)=0.
QT(JSTF)=F(I,JSTF)
DO 370 J=JST,M2
DENOM=AP(I,J)-PT(J-1)*AJM(I,J)
PT(J)=AJP(I,J)/DENOM
TEMP=CON(I,J)+AIP(I,J)*F(I+1,J)+AIM(I,J)*F(I-1,J)
QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM
370 END DO
DO 380 JJ=JST,M2
J=JT1-JJ
380 F(I,J)=F(I,J+1)*PT(J)+QT(J)
390 END DO
!----------------------------------------------------------------------
999 END DO
SELECT CASE(NF)

CASE(1)
u = f
CASE(2)
v = f
CASE(3)
pc = f
CASE(4)
t = f
CASE(np)
p = f
CASE(6)
conc = f
END SELECT
CALL RESET
END SUBROUTINE SOLVE

!================================================= =====================

[siddiquesil]

!================================================= =====================
!THIS SUBROUTINE RESETS 'ap' AND 'con' TO zero
!************************************************* *********************
SUBROUTINE RESET
USE GLOBAL
IMPLICIT NONE
INTEGER :: I,J

DO J=2,M2
DO I=2,L2
CON(I,J)=0.
AP(I,J)=0.
END DO
END DO
END SUBROUTINE RESET

!================================================= =====================

[siddiquesil]

!================================================= =====================
! THIS SUBROUTINE GENERATES THE CONTROL VOLUMES&VARIOUS RELATED GEOMETRICAL PARAMETERS
!************************************************* *********************
SUBROUTINE GEOMET
USE GLOBAL
IMPLICIT NONE
INTEGER :: I,J,IM4,IM5

1 FORMAT(//15X,'COMPUTATION IN CARTESIAN COORDINATES')
2 FORMAT(//15X,'COMPUTATION FOR AXISYMMETRIC SITUATION')
3 FORMAT(//15X,'COMPUTATION IN POLAR COORDINATES')
4 FORMAT(14X,40(1H*),//)
!-----------------------------------------------------------------------
NRHO=NP+1
NGAM=NRHO+1
L2=L1-1
L3=L2-1
M2=M1-1
M3=M2-1
X(1)=XU(2)
DO 5 I=2,L2
5 X(I)=0.5*(XU(I+1)+XU(I))
X(L1)=XU(L1)
Y(1)=YV(2)
DO 10 J=2,M2
10 Y(J)=0.5*(YV(J+1)+YV(J))
Y(M1)=YV(M1)
DO 15 I=2,L1
15 XDIF(I)=X(I)-X(I-1)
DO 18 I=2,L2
18 XCV(I)=XU(I+1)-XU(I)
DO 20 I=3,L2
20 XCVS(I)=XDIF(I)
XCVS(3)=XCVS(3)+XDIF(2)
XCVS(L2)=XCVS(L2)+XDIF(L1)
DO 22 I=3,L3
XCVI(I)=0.5*XCV(I)
22 XCVIP(I)=XCVI(I)
XCVIP(2)=XCV(2)
XCVI(L2)=XCV(L2)
DO 35 J=2,M1
35 YDIF(J)=Y(J)-Y(J-1)
DO 40 J=2,M2
40 YCV(J)=YV(J+1)-YV(J)
DO 45 J=3,M2
45 YCVS(J)=YDIF(J)
YCVS(3)=YCVS(3)+YDIF(2)
YCVS(M2)=YCVS(M2)+YDIF(M1)
IF(MODE.NE.1) GO TO 55
DO 52 J=1,M1
RMN(J)=1.0
52 R(J)=1.0
GO TO 56
55 DO 50 J=2,M1
50 R(J)=R(J-1)+YDIF(J)
RMN(2)=R(1)
DO 60 J=3,M2
60 RMN(J)=RMN(J-1)+YCV(J-1)
RMN(M1)=R(M1)
56 CONTINUE
DO 57 J=1,M1
SX(J)=1.
SXMN(J)=1.
IF(MODE.NE.3) GO TO 57
SX(J)=R(J)
IF(J.NE.1) SXMN(J)=RMN(J)
57 END DO
DO 62 J=2,M2
YCVR(J)=R(J)*YCV(J)
ARX(J)=YCVR(J)
IF(MODE.NE.3) GO TO 62
ARX(J)=YCV(J)
62 END DO
DO 64 J=4,M3
64 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J)
YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3)
YCVRS(M2)=0.5*(R(M1)+R(M3))*YCVS(M2)
IF(MODE.NE.2) GO TO 67
DO 65 J=3,M3
ARXJ(J)=0.25*(1.+RMN(J)/R(J))*ARX(J)
65 ARXJP(J)=ARX(J)-ARXJ(J)
GO TO 68
67 DO 66 J=3,M3
ARXJ(J)=0.5*ARX(J)
66 ARXJP(J)=ARXJ(J)
68 ARXJP(2)=ARX(2)
ARXJ(M2)=ARX(M2)
DO 70 J=3,M3
FV(J)=ARXJP(J)/ARX(J)
70 FVP(J)=1.-FV(J)
DO 85 I=3,L2
FX(I)=0.5*XCV(I-1)/XDIF(I)
85 FXM(I)=1.-FX(I)
FX(2)=0.
FXM(2)=1.
FX(L1)=1.
FXM(L1)=0.
DO 90 J=3,M2
FY(J)=0.5*YCV(J-1)/YDIF(J)
90 FYM(J)=1.-FY(J)
FY(2)=0.
FYM(2)=1.
FY(M1)=1.
FYM(M1)=0.
!----------------------------------------------------------------------
IM4=M1-3
IM5=M1-4
117 FORMAT(5E10.4)
IF(MODE.EQ.1) PRINT 1
IF(MODE.EQ.2) PRINT 2
IF(MODE.EQ.3) PRINT 3
PRINT 4
open(unit=20,file="grids.dat")
write(20,*)l1,m1
write(20,*)x
write(20,*)y
close(20)

END SUBROUTINE GEOMET

!================================================= =====================

[siddiquesil]

!================================================= =====================
!THIS SUBROUTINE FORMS COEFFS. FOR DISCRETISATION EQNS.
!************************************************* *********************
SUBROUTINE COEFF
USE GLOBAL
IMPLICIT NONE
INTEGER :: I,J,K,N
REAL :: ZERO,REL,FL,FLM,FLP,GM,GMM,VOL,APT,AREA
REAL :: SXT,SXB,ARHO,TMX,UMX,VMX,DELT,DELU,DELV,BTIME,BBK
CHARACTER(LEN=80) :: fname
!**************** COEFFICIENTS FOR THE U EQUATION **************
ITERL=1
ZERO=0.0
lconv = .false.
99 continue
CALL RESET
!----------------------------------------------------------------------
NF=1
IF(.NOT.LSOLVE(NF)) GO TO 500
IST=3
JST=2
CALL GAMSOR
REL=1.-RELAX(NF)
DO 102 I=3,L2
FL=XCVI(I)*V(I,2)*RHO(I,1)
FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1)
FLOW=R(1)*(FL+FLM)
DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2)
CALL DIFLOW
102 AJM(I,2)=ACOF+AMAX1(ZERO,FLOW)
DO 103 J=2,M2
FLOW=ARX(J)*U(2,J)*RHO(1,J)
DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J))
CALL DIFLOW
AIM(3,J)=ACOF+AMAX1(ZERO,FLOW)
DO 103 I=3,L2
IF(I.EQ.L2) GO TO 104
FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J))
FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J) )
FLOW=ARX(J)*0.5*(FL+FLP)
DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J))
GO TO 105
104 FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J))
105 CALL DIFLOW
AIM(I+1,J)=ACOF+AMAX1(ZERO,FLOW)
AIP(I,J)=AIM(I+1,J)-FLOW
IF(J.EQ.M2) GO TO 106
FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*R HO(I,J))
FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)* &
& RHO(I-1,J))
GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+ &
& 1.0E-30)*XCVI(I)
GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)* &
& GAM(I-1,J)+1.E-30)*XCVIP(I-1)
DIFF=RMN(J+1)*2.*(GM+GMM)
GO TO 107
106 FL=XCVI(I)*V(I,M1)*RHO(I,M1)
FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1)
DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1)
107 FLOW=RMN(J+1)*(FL+FLM)
CALL DIFLOW
AJM(I,J+1)=ACOF+AMAX1(ZERO,FLOW)
AJP(I,J)=AJM(I,J+1)-FLOW
VOL=YCVR(J)*XCVS(I)
APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1)) &
&/(XCVS(I)*DT)
AP(I,J)=AP(I,J)-APT
CON(I,J)=CON(I,J)+APT*UO(I,J)
AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) &
&/RELAX(NF)
CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*U(I,J)
DU(I,J)=VOL/(XDIF(I)*SX(J))
DU(I,J)=DU(I,J)/AP(I,J)
103 CONTINUE
!----------------------------------------------------------------------
cofu(:,:,1) = con(:,: )
cofu(:,:,2) = aip(:,: )
cofu(:,:,3) = aim(:,: )
cofu(:,:,4) = ajp(:,: )
cofu(:,:,5) = ajm(:,: )
cofu(:,:,6) = ap(:,: )

!***************** COEFFICIENTS FOR THE V EQUATION ***********
NF=2
CALL RESET
IST=2
JST=3
CALL GAMSOR
REL=1.-RELAX(NF)
DO 202 I=2,L2
AREA=R(1)*XCV(I)
FLOW=AREA*V(I,2)*RHO(I,1)
DIFF=AREA*GAM(I,1)/YCV(2)
CALL DIFLOW
202 AJM(I,3)=ACOF+AMAX1(ZERO,FLOW)
DO 203 J=3,M2
FL=ARXJ(J)*U(2,J)*RHO(1,J)
FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1)
FLOW=FL+FLM
DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J))
CALL DIFLOW
AIM(2,J)=ACOF+AMAX1(ZERO,FLOW)
DO 203 I=2,L2
IF(I.EQ.L2) GO TO 204
FL=ARXJ(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*R HO(I,J))
FLM=ARXJP(J-1)*U(I+1,J-1)*(FX(I+1)*RHO(I+1,J-1)+FXM(I+1)* &
& RHO(I,J-1))
GM=GAM(I,J)*GAM(I+1,J)/(XCV(I)*GAM(I+1,J)+XCV(I+1)*GAM(I,J)+ &
& 1.E-30)*ARXJ(J)
GMM=GAM(I,J-1)*GAM(I+1,J-1)/(XCV(I)*GAM(I+1,J-1)+XCV(I+1)* &
& GAM(I,J-1)+1.0E-30)*ARXJP(J-1)
DIFF=2.*(GM+GMM)/SXMN(J)
GO TO 205
204 FL=ARXJ(J)*U(L1,J)*RHO(L1,J)
FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1)
DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J))
205 FLOW=FL+FLM
CALL DIFLOW
AIM(I+1,J)=ACOF+AMAX1(ZERO,FLOW)
AIP(I,J)=AIM(I+1,J)-FLOW
IF(J.EQ.M2) GO TO 206
AREA=R(J)*XCV(I)
FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J)
FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J) )*RMN(J+1)
FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I)
DIFF=AREA*GAM(I,J)/YCV(J)
GO TO 207
206 AREA=R(M1)*XCV(I)
FLOW=AREA*V(I,M1)*RHO(I,M1)
DIFF=AREA*GAM(I,M1)/YCV(M2)
207 CALL DIFLOW
AJM(I,J+1)=ACOF+AMAX1(ZERO,FLOW)
AJP(I,J)=AJM(I,J+1)-FLOW
VOL=YCVRS(J)*XCV(I)
SXT=SX(J)
IF(J.EQ.M2) SXT=SX(M1)
SXB=SX(J-1)
IF(J.EQ.3) SXB=SX(1)
APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*RHO(I,J-1)* &
&0.5*(SXB+SXMN(J)))/(YCVRS(J)*DT)
AP(I,J)=AP(I,J)-APT
CON(I,J)=CON(I,J)+APT*VO(I,J)
AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) &
&/RELAX(NF)
CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*V(I,J)
DV(I,J)=VOL/YDIF(J)
DV(I,J)=DV(I,J)/AP(I,J)
203 CONTINUE
!----------------------------------------------------------------------
cofv(:,:,1) = con(:,: )
cofv(:,:,2) = aip(:,: )
cofv(:,:,3) = aim(:,: )
cofv(:,:,4) = ajp(:,: )
cofv(:,:,5) = ajm(:,: )
cofv(:,:,6) = ap(:,: )
!******************* CALCULATE UHAT AND VHAT ********************
DO J=2,M2
DO I=3,L2
UHAT(I,J)=(COFU(I,J,2)*U(I+1,J)+COFU(I,J,3)*U(I-1,J)+COFU(I,J,4) &
& *U(I,J+1)+COFU(I,J,5)*U(I,J-1)+COFU(I,J,1))/COFU(I,J,6)
END DO
END DO
DO J=3,M2
DO I=2,L2
VHAT(I,J)=(COFV(I,J,2)*V(I+1,J)+COFV(I,J,3)*V(I-1,J)+COFV(I,J,4) &
& *V(I,J+1)+COFV(I,J,5)*V(I,J-1)+COFV(I,J,1))/COFV(I,J,6)
END DO
END DO
!******************** COEFFICIENTS FOR THE PRESSURE EQUATION *****
NF=3
CALL RESET
IST=2
JST=2
CALL GAMSOR
DO 410 J=2,M2
DO 410 I=2,L2
VOL=YCVR(J)*XCV(I)
410 CON(I,J)=CON(I,J)*VOL
DO 402 I=2,L2
ARHO=R(1)*XCV(I)*RHO(I,1)
CON(I,2)=CON(I,2)+ARHO*V(I,2)
402 AJM(I,2)=0.
DO 403 J=2,M2
ARHO=ARX(J)*RHO(1,J)
CON(2,J)=CON(2,J)+ARHO*U(2,J)
AIM(2,J)=0.
DO 403 I=2,L2
IF(I.EQ.L2) GO TO 404
ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
FLOW=ARHO*UHAT(I+1,J)
CON(I,J)=CON(I,J)-FLOW
CON(I+1,J)=CON(I+1,J)+FLOW
AIP(I,J)=ARHO*DU(I+1,J)
AIM(I+1,J)=AIP(I,J)
GO TO 405
404 ARHO=ARX(J)*RHO(L1,J)
CON(I,J)=CON(I,J)-ARHO*U(L1,J)
AIP(I,J)=0.
405 IF(J.EQ.M2) GO TO 406
ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)* RHO(I,J))
FLOW=ARHO*VHAT(I,J+1)
CON(I,J)=CON(I,J)-FLOW
CON(I,J+1)=CON(I,J+1)+FLOW
AJP(I,J)=ARHO*DV(I,J+1)
AJM(I,J+1)=AJP(I,J)
GO TO 407
406 ARHO=RMN(M1)*XCV(I)*RHO(I,M1)
CON(I,J)=CON(I,J)-ARHO*V(I,M1)
AJP(I,J)=0.
407 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)
403 continue
DO 421 J=2,M2
DO 421 I=2,L2
AP(I,J)=AP(I,J)/RELAX(NP)
421 CON(I,J)=CON(I,J)+(1.0-RELAX(NP))*AP(I,J)*P(I,J)
!----------------------------------------------------------------------
cofp(:,:,2) = aip(:,: )
cofp(:,:,3) = aim(:,: )
cofp(:,:,4) = ajp(:,: )
cofp(:,:,5) = ajm(:,: )

NF=NP
CALL SOLVE

!******************** COMPUTE U* AND V* **************************
NF=1
IST=3
JST=2
DO 550 K=1,6
DO 550 J=JST,M2
DO 550 I=IST,L2
550 COF(I,J,K)=COFU(I,J,K)
!----------------------------------------------------------------------
con(:,: ) = cof(:,:,1)
aip(:,: ) = cof(:,:,2)
aim(:,: ) = cof(:,:,3)
ajp(:,: ) = cof(:,:,4)
ajm(:,: ) = cof(:,:,5)
ap(:,: ) = cof(:,:,6)
!----------------------------------------------------------------------
DO 551 J=JST,M2
DO 551 I=IST,L2
551 CON(I,J)=CON(I,J)+DU(I,J)*AP(I,J)*(P(I-1,J)-P(I,J))
CALL SOLVE

NF=2
IST=2
JST=3
DO 552 K=1,6
DO 552 J=JST,M2
DO 552 I=IST,L2
552 COF(I,J,K)=COFV(I,J,K)
!----------------------------------------------------------------------
con(:,: ) = cof(:,:,1)
aip(:,: ) = cof(:,:,2)
aim(:,: ) = cof(:,:,3)
ajp(:,: ) = cof(:,:,4)
ajm(:,: ) = cof(:,:,5)
ap(:,: ) = cof(:,:,6)
!----------------------------------------------------------------------
DO 553 J=JST,M2
DO 553 I=IST,L2
553 CON(I,J)=CON(I,J)+DV(I,J)*AP(I,J)*(P(I,J-1)-P(I,J))
CALL SOLVE

!*************** COEFFICIENTS FOR THE PRESSURE CORRECTION EQUATION ***
NF=3
CALL RESET
IST=2
JST=2
DO 554 K=2,5
DO 554 J=JST,M2
DO 554 I=IST,L2
554 COF(I,J,K)=COFP(I,J,K)
!----------------------------------------------------------------------
aip(:,: ) = cof(:,:,2)
aim(:,: ) = cof(:,:,3)
ajp(:,: ) = cof(:,:,4)
ajm(:,: ) = cof(:,:,5)
!----------------------------------------------------------------------
CALL GAMSOR
SMAX=0.
SSUM=0.
DO 510 J=2,M2
DO 510 I=2,L2
VOL=YCVR(J)*XCV(I)
510 CON(I,J)=CON(I,J)*VOL
DO 502 I=2,L2
ARHO=R(1)*XCV(I)*RHO(I,1)
CON(I,2)=CON(I,2)+ARHO*V(I,2)
502 END DO
DO 503 J=2,M2
ARHO=ARX(J)*RHO(1,J)
CON(2,J)=CON(2,J)+ARHO*U(2,J)
DO 503 I=2,L2
IF(I.EQ.L2) GO TO 504
ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
FLOW=ARHO*U(I+1,J)
CON(I,J)=CON(I,J)-FLOW
CON(I+1,J)=CON(I+1,J)+FLOW
GO TO 505
504 ARHO=ARX(J)*RHO(L1,J)
CON(I,J)=CON(I,J)-ARHO*U(L1,J)
505 IF(J.EQ.M2) GO TO 506
ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)* RHO(I,J))
FLOW=ARHO*V(I,J+1)
CON(I,J)=CON(I,J)-FLOW
CON(I,J+1)=CON(I,J+1)+FLOW
GO TO 507
506 ARHO=RMN(M1)*XCV(I)*RHO(I,M1)
CON(I,J)=CON(I,J)-ARHO*V(I,M1)
507 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)
PC(I,J)=0.
SMAX=AMAX1(SMAX,ABS(CON(I,J)))
SSUM=SSUM+CON(I,J)
503 CONTINUE
CALL SOLVE

!*************** COME HERE TO CORRECT THE VELOCITIES ************
DO 521 J=2,M2
DO 521 I=2,L2
IF(I.NE.2) U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J))
IF(J.NE.2) V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J))
521 CONTINUE
500 CONTINUE

!************** COEFFICIENTS FOR OTHER EQUATIONS ****************
CALL RESET
IST=2
JST=2
DO 600 N=4,NFMAX
NF=N
IF(.NOT.LSOLVE(NF)) GO TO 600
CALL GAMSOR
REL=1.-RELAX(NF)
DO 602 I=2,L2
AREA=R(1)*XCV(I)
FLOW=AREA*V(I,2)*RHO(I,1)
DIFF=AREA*GAM(I,1)/YDIF(2)
CALL DIFLOW
602 AJM(I,2)=ACOF+AMAX1(ZERO,FLOW)
DO 603 J=2,M2
FLOW=ARX(J)*U(2,J)*RHO(1,J)
DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J))
CALL DIFLOW
AIM(2,J)=ACOF+AMAX1(ZERO,FLOW)
DO 603 I=2,L2
IF(I.EQ.L2) GO TO 604
FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)* RHO(I,J))
DIFF=ARX(J)*2.*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+ &
& XCV(I+1)*GAM(I,J)+1.0E-30)*SX(J))
GO TO 605
604 FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J))
605 CALL DIFLOW
AIM(I+1,J)=ACOF+AMAX1(ZERO,FLOW)
AIP(I,J)=AIM(I+1,J)-FLOW
AREA=RMN(J+1)*XCV(I)
IF(J.EQ.M2) GO TO 606
FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RH O(I,J))
DIFF=AREA*2.*GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+ &
& YCV(J+1)*GAM(I,J)+1.0E-30)
GO TO 607
606 FLOW=AREA*V(I,M1)*RHO(I,M1)
DIFF=AREA*GAM(I,M1)/YDIF(M1)
607 CALL DIFLOW
AJM(I,J+1)=ACOF+AMAX1(ZERO,FLOW)
AJP(I,J)=AJM(I,J+1)-FLOW
603 CONTINUE
!----------------------------------------------------------------------
DO 349 I=1,L1
DO 349 J=1,M1
VOL=YCVR(J)*XCV(I)
APT=RHO(I,J)/DT
AP(I,J)=AP(I,J)-APT
if(nf.eq.4) CON(I,J)=CON(I,J)+APT*TO(I,J)

AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) &
&/RELAX(NF)
if(nf.eq.4) CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*t(I,J)
if(nf.eq.np) CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*p(I,J)

!----------------------------------------------------------------------
if(nf.eq.4) TOLD(i,j)=T(i,j)
!----------------------------------------------------------------------

349 CONTINUE
CALL SOLVE
!----------------------------------------------------------------------

600 END DO

!*****************---------------CONVERGENCE CRITERIA----------------------------
tmx = maxval(t)
umx = maxval(u)
vmx = maxval(v)

delt = maxval(abs(t-told))
delu = maxval(abs(u-uold))
delv = maxval(abs(v-vold))

if(tmx.gt.0.0) then
deltmx = delt/tmx
else
deltmx=0.0
endif
if(umx.gt.0.0) then
delumx = delu/umx
else
delumx=0.0
endif
if(vmx.gt.0.0) then
delvmx = delv/vmx
else
delvmx=0.0
endif


if(time.lt.120.0) then
WRITE(*,*)'\n',iterl,'tmx=',tmx,'delt=',deltmx,'um x=',umx,'delu=',delumx,'vmx=',vmx,'delv=',delvmx
endif

if(iterl.ge.istp) then
if((delt.lt.0.15).and.(delu.lt.0.005)) then
lconv = .true.
endif
endif

if(lconv) goto 59
IF(DELTMX.GT.ERT) GOTO 51
IF(DELUMX.GT.ERU) GOTO 51
IF(DELVMX.GT.ERV) GOTO 51

write(*,*)'converged'
goto 59


51 IF(ITERL.EQ.ISTP) THEN
WRITE(*,*)'exceeded max. no. of iteration'
WRITE(*,*)'\n',iterl,'tmx=',tmx,'delt=',deltmx,'um x=',umx,'delu=',delumx,'vmx=',vmx,'delv=',delvmx
GOTO 59
STOP
ENDIF
40 FORMAT(I5,6F10.6,2X,2F12.4)

ITERL=ITERL+1

DO 765 I=2,L1
DO 765 J=2,M1
UOLD(I,J)=U(I,J)
VOLD(I,J)=V(I,J)
765 CONTINUE

CALL BOUND
GOTO 99
59 ITLL=ITLL+1
BTIME=TIME+DT
write(*,*)'time=',time,'end of inner iterations'
fname="velocity.dat"
open(UNIT=10,file=fname)
write(10,*)u
write(10,*)v
close(10)

TIME=TIME+DT


!plotting every 50 timestep
IF(MOD(INT(TIME/DT),50).eq.0) THEN
CALL INTERPOLATION

ENDIF
CALL ENRBAL

BBK=AMOD(TIME,real(NTIME))
IF(BBK.LE.0.00001)THEN
ITER=ITER+1
CALL BOUND
END IF

IF(TIME.GE.TLAST) LSTOP=.TRUE.
DO 775 I=2,L1
DO 775 J=2,M1
UOLD(I,J)=U(I,J)
VOLD(I,J)=V(I,J)
775 CONTINUE
END SUBROUTINE COEFF

!================================================= =====================

[siddiquesil]

!================================================= =====================
!this subroutine stores initial values for the coming timestep
!************************************************* *********************
subroutine oldval
use global
implicit none
integer :: I,j

do i=1,L1
do j=1,m1
to(i,j)=t(i,j)
uo(i,j)=u(i,j)
vo(i,j)=v(i,j)
end do
end do
end subroutine oldval

!================================================= =====================

[siddiquesil]

!================================================= =====================
!this subroutine generates uniform grid
!************************************************* *********************
subroutine ugrid
use global
implicit none
integer :: I,j
real :: Dx,dy

xu(2)=0.
Dx=xl/float(l1-2)
do 1 i=3,l1
1 xu(i)=xu(i-1)+dx
yv(2)=0.
Dy=yl/float(m1-2)
do 2 j=3,m1
2 yv(j)=yv(j-1)+dy
end subroutine ugrid

[siddiquesil]

!================================================= =====================
!THIS SUBROTINE CHECKS OVERALL ENERGY BALANCE
!************************************************* *********************
SUBROUTINE ENRBAL
USE GLOBAL
IMPLICIT NONE
END SUBROUTINE ENRBAL
!================================================= =====================
!THIS SUBROUTINE GIVES INTERPOLATED VELOCITIES OF PROBLEM VARIABLES
!************************************************* *********************
SUBROUTINE INTERPOLATION
USE GLOBAL
do j=1,m1
u(1,j)=u(2,j)
end do
do i=1,l1
v(i,1)=v(i,2)
end do
do i=2,l2
do j=1,m1
ur(i,j)=0.5*(u(i+1,j)+u(i,j))
end do
end do
do i=1,l1
do j=2,m2
vr(i,j)=0.5*(v(i,j+1)+v(i,j))
end do
end do
CALL PRINTOUT
END SUBROUTINE INTERPOLATION
!================================================= =====================
!THIS SUBROUTINE DIRECTLY GIVES MATLAB PLOTS OF PROBLEM VARIABLES
!************************************************* *********************
SUBROUTINE PRINTOUT
USE GLOBAL
character::filename*20,fname*20,fname1*20,fname2*2 0
integer:: ileunit,ile,ile1
!print*,'************test',MOD(INT(TIME/DT),20)
!plotting every 50 timestep
IF(MOD(INT(TIME/DT),50).ne.0) return
!CALCULATE THE STREAM FUNCTION---------------------------------------
pc(2,2)=0.
DO I=2,L1
IF(I.NE.2) pc(I,2)=pc(I-1,2)-RHO(I-1,1)*V(I-1,2)*R(1)*XCV(I-1)
DO J=3,M1
RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1)
pc(I,J)=pc(I,J-1)+RHOM*U(I,J-1)*ARX(J-1)
ENDDO
ENDDO

!if(time.le.2500) goto700
!goto 711
if(int(time).eq.500) filename='out500.m'

open(unit=55,file=FileName)
write(55,FMT='(A3,<l1>(E15.8,1X),A2)') "x=[" , x , "];"
write(55,FMT='(A3,<m1>(E15.8,1X),A2)') "y=[" , y , "];"


write(55,*) "t=["
do j=1,m1
write(55,FMT='(<l1>(E15.8,1X),A1)') t(:,j), ";"
end do
write(55,*) "];"
write(55,*) "pc=["
do j=1,m1
write(55,FMT='(<l1>(E15.8,1X),A1)') pc(:,j), ";"
end do
write(55,*) "];"

write(55,*) "u=["
do j=1,m1
write(55,FMT='(<l1>(E15.8,1X),A1)') ur(:,j), ";"
end do
write(55,*) "];"
write(55,*) "v=["
do j=1,m1
write(55,FMT='(<l1>(E15.8,1X),A1)') vr(:,j), ";"
end do
write(55,*) "];"

write(55,*) "figure(1)"
write(55,*) "contour(x,y,t);"
write(55,*) "figure(2)"
write(55,*) "quiver(x,y,u,v)"
write(55,*) "figure(3)"
write(55,*) "contour(x,y,pc)"
write(55,*) "title('Pictures at time=",time,"seconds');"




close(55)
! 711 continue
if(ltecplot) then
!techplot plotting
ileunit=int(time)
ile=int(time)
write(fname,1000)ileunit
1000 format('temp.',i5.5)

open(UNIT=ileunit,FILE=fname) !directs to an output file
write(ileunit,*) "VARIABLES = 'X','Y','T'"
write(ileunit,*) "ZONE I=",l1,"J=",m1
do j=1,m1
do i=1,l1
write(ileunit,702)(x(i)),(y(j)),(t(i,j))
702 format(1x,f8.4,2x,f8.4,2x,f8.4)
enddo
enddo
close(ileunit)
!techplot plotting of velocity
write(fname1,1001)ile
1001 format('vel.',i5.5)


OPEN(UNIT=ile,FILE=fname1)

write(ile,*) "VARIABLES = 'X','Y','U','V'"
write(ile,*) "ZONE I=",l1,"J=",m1
do j=1,m1
do i=1,l1
write(ile,701)(x(i)),(y(j)),(ur(i,j)),(vr(i,j))
701 format(1x,f8.4,2x,f8.4,2x,f8.4,2x,f8.4)
! 701 format((1x,4f10.2),A1)
enddo
enddo


close(ile)

ile1=int(time)
write(fname2,1002)ile1
1002 format('stream.',i5.5)

open(UNIT=ile1,FILE=fname2) !directs to an output file
write(ile1,*) "VARIABLES = 'X','Y','PC'"
write(ile1,*) "ZONE I=",l1,"J=",m1
do j=1,m1
do i=1,l1
write(ile1,703)(x(i)),(y(j)),(pc(i,j))
703 format(1x,f8.4,2x,f8.4,2x,f8.4)
enddo
enddo
close(ile1)

endif





END SUBROUTINE PRINTOUT

!================================================= =====================
!================= USER PART BEGINS HERE ==========================
!================================================= =====================

!================================================= =====================
!THIS SUBROUTINE INITIALISES DIFFERENT VARIABLES AND READS INPUT DATA
!************************************************* *********************
SUBROUTINE GRID
USE GLOBAL
IMPLICIT NONE
INTEGER :: I,J

IPREF=1
JPREF=1
MODE=1
RELAX =0.6
LSOLVE = .FALSE.
LSOLVE(1)=.TRUE.
LSOLVE(2)=.TRUE.
LSOLVE(3)=.TRUE.
LSOLVE(4)=.TRUE.
LSOLVE(11)=.TRUE.
LSTOP=.FALSE.

DO I=1,NFX3
NTIMES(I)=1
LBLK(I)=.TRUE.
END DO

TIME=0.
TITLE(1)='UVEL'
TITLE(2)='VVEL'
TITLE(4)='TEMP'
TITLE(6)='CONC'

deltmx = 1.0

!================================================= =======================
!=========== meaning of symbols:
!This is a sample data file for free convection in a rect. cavity
!Symbol meanings are:
!xl : x-direction domain length
!yl : y-direction domain length
!l1 :no. of grids in x-dirn
!m1 :no. of grids in y-dirn
!ngrid : switch for grid type; if ngrid=1 then uniform grid, else non-uniform
!tk : thermal conductivity
!g : acceleration due to gravity
!beta : volumetric expansion coefficient
!cp : specific heat at const. press.
!tref : ref. temp. for buoyancy source term
!amiu : dynamic viscosity
!rhocon : density
!flux : heat flux at boundary
!relax(nf) : relaxation factor for a particular nf
!dt : time step size to start with
!tlast : total time domain
!time2 : time till the time step value of dt is valid
!dt2 : time step value after time t2; similarly dt3.....
!eru,erv,ert: convergence tolerances for u,v,t (relative error)
!istp : maximum no. of iterations in each timestep
!ntimes(nf) : no. of times corr. to which eqn corr. to a particular nf is solved
!tinitial :initial value of temp.for a transient problem
!rhoinitial :initial value of density for a transient problem

!***************** INPUT DATA HERE ***************************

XL=0.22 ; YL=0.11 ; L1=40;M1=40
NTIME=1 ; NGRID=1 ; ISTP=500
tk=0.62 ; cp=4179.0
rhocon=997.0 ;RHOINITIAL=997.0
DT=10.0 ; TLAST=510.0
ERU=0.00001; ERV=0.00001 ; ERT=0.000001
tinitial=23.0 ;TREF=23.0; BETA=0.0002761 ; AMIU=0.000855 ;FLUX=500.0
g=9.81; relax(1)=0.5; relax(2)=0.5; relax(3)=0.5; relax(4)=0.5
ltecplot=.false.
!----------------------------------------------------------------------
ITER=0
R(1)=0.0
CALL ugrid

!************ CON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE **
DO J=1,M1
DO I=1,L1
U(I,J)=0.
V(I,J)=0.
PC(I,J)=0.
CONC(i,j)=0.
P(I,J)=0.
RHO(I,J)= rhocon
DU(I,J)=0.
DV(I,J)=0.
UOLD(I,J)=0.
VOLD(I,J)=0.
AIP(I,J)=0.
AIM(I,J)=0.
AJM(I,J)=0.
AJP(I,J)=0.
CON(I,J)=0.
AP(I,J)=0.
AP0(I,J)=0.
told(i,j)=0.
END DO
END DO

cofu = 0.0
cofv = 0.0
cof = 0.0
cofp = 0.0

END SUBROUTINE grid

!======================================================================
!THIS SUBROUTINE GIVES INITIAL CONDITIONS FOR THE PROBLEM
!************************************************* *********************
SUBROUTINE START
USE GLOBAL
IMPLICIT NONE
INTEGER :: I,J
!real:;temp
DO I=1,L1
DO J=1,M1
U(I,J)=0.0
V(I,J)=0.0
T(I,J)=TINITIAL
TOLD(I,J)=tinitial
rho(i,j)=rhocon
END DO
END DO
END SUBROUTINE start

!================================================= =====================
!THIS SUBROUTINE CALCULATES THE DENSITY
!************************************************* *********************
SUBROUTINE dense
USE GLOBAL
IMPLICIT NONE
INTEGER :: I,J

DO I=1,L1
DO J=1,M1
rho(i,j)=rhoinitial
END DO
END DO
END SUBROUTINE DENSE

!================================================= =====================
!THIS SUBROUTINE GIVES BOUNDARY CONDITIONS FOR THE PROBLEM
!************************************************* *********************
SUBROUTINE bound
USE GLOBAL
IMPLICIT NONE
INTEGER::I,J

!U-VELOCITY BOUNDARY CONDITION
DO I=2,L1
DO J=1,M1
!TOP & BOTTOM
u(i,m1)=0.0
U(I,1)=0.0
!SIDES
U(2,J)=0.0
U(L1,J)=0.0
END DO
END DO
DO I=1,L1
DO J=2,M1
V(1,J)=0.0
V(L1,J)=0.0
V(I,2)=0.0
V(I,M1)=0.0
END DO
END DO
!TEMPERATURE BOUNDARY CONDITION

DO J=1,M1
!LEFT
T(1,J)=T(2,J)+FLUX*XDIF(2)/TK
!SIDES
T(L1,J)=T(L2,J)+(FLUX*XDIF(L1)/TK)
ENDDO
DO I=1,L1
!BOTTOM AND TOP
T(I,1)=T(I,2)
t(i,m1)=t(i,m2)

END DO


END SUBROUTINE bound

!================================================= =====================
!THIS SUBROUTINE GIVES DIFF. COEFF. & SOURCE TERMS FOR DISCRETISED EQNS.
!************************************************* *********************
SUBROUTINE gamsor
USE GLOBAL
IMPLICIT NONE
INTEGER :: I,J
REAL:: TAV

DO I=1,L1
DO J=1,M1
IF(nf.EQ.4) gam(i,j)=tk/cp
IF(NF/=4) GAM(I,J)=AMIU
END DO
END DO
IF(NF==2) THEN
DO I=2,L2
DO J=3,M2
Tav=FY(J)*T(I,J)+FYM(J)*T(I,J-1)
CON(I,J)=CON(I,J)+rhoinitial*g*BETA*(Tav-tref)

END DO
END DO
ENDIF

if(nf.eq.4) then
do i=1,l1
gam(i,m1)=0.0
enddo
endif

END SUBROUTINE gamsor

!================================================= =====================
!******************** END OF THE CODE ************************
!================================================= =====================

[vesp]

wtf? what is the purpose of this post?

[Hayder Mohammed]

Dear sir I need a fortran subroutine to solve melting and solidification by using enthalpy method ...finite volume or finite difference method ..thanks

[siddiquesil]

You can do melting and solidification using this code.
U just need to add source term and change the boundary condition here in this code.

The rest all are same.

Read Patankar Book to understand this code.

best of luck

[praveen]

Please do not post code like this. It is not a good way to share code.

You can put the code here (assuming you have the rights to do it)

https://www.cfd-online.com/Wiki/Source_code_archive

and make a new post and link to it.

[siddiquesil]

Ok.
Thank you sir. I will do it

Sir, I have converted this code for spherical coordinate.
Work title"Flow over a sphere under moderate Reynolds number"
My U-velocity is converging but V-velocity is not converging.
What would be the reason?

[Hayder Mohammed]

which code do you mean sir ?....actually I want to simulate a heat transfer through wall subjected to a solar radiation and heat convection.
the wall consists of three layers : the outer layer is a mortar followed by a layer of brick +phase change material and the inner layer is a mortar.
I want to do this simulation using FORTRAN 90 code using enthalpy -finite volume method
so if you help me I'll be very grateful to you.

[siddiquesil]

See, first u try to simulate a cavity filled with phased change material.
Upper and lower side of the cavity are insulated. Left side temperature is higher than the melting point temp and right side is kept at zero degree.

[amcic1996]

Originally Posted by siddiquesil

Ok.
Thank you sir. I will do it

Sir, I have converted this code for spherical coordinate.
Work title"Flow over a sphere under moderate Reynolds number"
My U-velocity is converging but V-velocity is not converging.
What would be th pornjk porn800 redtubee reason?

Please go through this brief instruction before attempting anything with this code

[xiejin]

Hi Sir,
Do you know how to convert Fluent mesh file format to Fortran FVM use?