From CFD-Wiki
!Sample program for solving Smith-Hutton Test using different schemes
!of covective terms approximation - Coefficients computing modul
!Copyright (C) 2005 Michail Kiričkov
!This program is free software; you can redistribute it and/or
!modify it under the terms of the GNU General Public License
!as published by the Free Software Foundation; either version 2
!of the License, or (at your option) any later version.
!This program is distributed in the hope that it will be useful,
!but WITHOUT ANY WARRANTY; without even the implied warranty of
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!GNU General Public License for more details.
!You should have received a copy of the GNU General Public License
!along with this program; if not, write to the Free Software
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
!*********************************************************************
Subroutine Coef_1(nf)
include 'icomm_1.f90'
Dimension F_out(nx,ny),CheckFlux(nx,ny)
Character Filename*10
! calculation of fluxes
! all geometry has rectangular 2D notation
Do 100 I= 2,NXmax
Do 100 J= 2,NYmax
Gam_e = ( Gam(i+1,j ) + Gam(i ,j ) ) * 0.5
Gam_w = ( Gam(i-1,j ) + Gam(i ,j ) ) * 0.5
Gam_s = ( Gam(i ,j-1) + Gam(i ,j ) ) * 0.5
Gam_n = ( Gam(i ,j+1) + Gam(i ,j ) ) * 0.5
Ro_e = ( Ro(i+1,j ) + Ro(i ,j ) ) * 0.5
Ro_w = ( Ro(i-1,j ) + Ro(i ,j ) ) * 0.5
Ro_s = ( Ro(i ,j-1) + Ro(i ,j ) ) * 0.5
Ro_n = ( Ro(i ,j+1) + Ro(i ,j ) ) * 0.5
Area_w = Y_xi(i-1,j-1)
Area_e = Y_xi(i ,j-1)
Area_s = X_et(i-1,j-1)
Area_n = X_et(i-1,j )
Del_e = Del_X_et(i ,j )
Del_w = Del_X_et(i-1,j )
Del_n = Del_Y_xi(i ,j )
Del_s = Del_Y_xi(i ,j-1)
! upwind differencing (all other will be included into the source term)
Con_e = Area_e * ( F(i,j,1) + F(i+1,j ,1) ) * 0.5
Con_w = Area_w * ( F(i,j,1) + F(i-1,j ,1) ) * 0.5
Con_s = Area_s * ( F(i,j,2) + F(i ,j-1,2) ) * 0.5
Con_n = Area_n * ( F(i,j,2) + F(i ,j+1,2) ) * 0.5
Diff_e = Area_e * Gam_e / Del_e
Diff_w = Area_w * Gam_w / Del_w
Diff_s = Area_s * Gam_s / Del_s
Diff_n = Area_n * Gam_n / Del_n
Flux_e = Area_e * Con_e * Ro_e
Flux_w = Area_w * Con_w * Ro_w
Flux_s = Area_s * Con_s * Ro_s
Flux_n = Area_n * Con_n * Ro_n
Aw(i,j) = Diff_w + max( Con_w,0.)
Ae(i,j) = Diff_e + max(-1.* Con_e,0.)
As(i,j) = Diff_s + max( Con_s,0.)
An(i,j) = Diff_n + max(-1.* Con_n,0.)
CheckFlux(i,j) = Flux_e - Flux_w + Flux_s - Flux_n
Ap(i,j)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) + CheckFlux(i,j)
Sp(i,j)= 0.
!--------------------------------QUICK SCHEME-------------------------
go to 700 (now QUICK is "off")
if( (i.GT.2).AND.(i.LT.NXmax).and.(j.GT.4).AND.(j.LT.NYmax) ) then
if(Con_e.GT.0.) Sp(i,j) = Sp(i,j) + Con_e * (-1.) * (-0.125 * F(i-1,j ,nf) - 0.25 * F(i ,j ,nf) + 0.375 * F(i-1,j ,nf) )
if(Con_w.GT.0.) Sp(i,j) = Sp(i,j) + Con_w * (-1.) * (-0.125 * F(i-2,j ,nf) - 0.25 * F(i-1,j ,nf) + 0.375 * F(i-2,j ,nf) )
if(Con_n.GT.0.) Sp(i,j) = Sp(i,j) + Con_n * (-1.) * (-0.125 * F(i ,j-1,nf) - 0.25 * F(i ,j ,nf) + 0.375 * F(i ,j-1,nf) )
if(Con_s.GT.0.) Sp(i,j) = Sp(i,j) + Con_s * (-1.) * (-0.125 * F(i ,j-2,nf) - 0.25 * F(i ,j-1,nf) + 0.375 * F(i ,j-2,nf) )
if(Con_e.LT.0.) Sp(i,j) = Sp(i,j) + Con_e * (-1.) * ( 0.375 * F(i ,j ,nf) - 0.25 * F(i+1,j ,nf) - 0.125 * F(i+2,j ,nf) )
if(Con_w.LT.0.) Sp(i,j) = Sp(i,j) + Con_w * (-1.) * ( 0.375 * F(i-1,j ,nf) - 0.25 * F(i ,j ,nf) - 0.125 * F(i+1,j ,nf) )
if(Con_n.LT.0.) Sp(i,j) = Sp(i,j) + Con_n * (-1.) * ( 0.375 * F(i ,j ,nf) - 0.25 * F(i ,j+1,nf) - 0.125 * F(i ,j+2,nf) )
! if(Con_s.LT.0.) Sp(i,j) = Sp(i,j) + Con_s * (-1.) * ( 0.375 * F(i ,j-1,nf) - 0.25 * F(i ,j+0,nf) - 0.125 * F(i ,j+1,nf) )
end if
700 continue
!-------------------------------- HLPA SCHEME----------------------------
go to 600 (now HLPA is "off")
! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f)
if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then
!------------------ w face -------------------
Fww = F(i-2,j,nf)
Fw = F(i-1,j,nf)
Fp = F(i ,j,nf)
Fe = F(i+1,j,nf)
call HLPA(Con_w,Fww,Fw,Fp,Fe,Delta_f)
Sp(i,j) = Sp(i,j) + Con_w * Delta_f
!------------------ e face--------------------
Fww = F(i-1,j,nf)
Fw = F(i ,j,nf)
Fp = F(i+1,j,nf)
Fe = F(i+2,j,nf)
call HLPA(Con_e,Fww,Fw,Fp,Fe,Delta_f)
Sp(i,j) = Sp(i,j) + Con_e * Delta_f * (-1.)
!------------------ s face--------------------
Fww = F(i ,j-2,nf)
Fw = F(i ,j-1,nf)
Fp = F(i ,j ,nf)
Fe = F(i ,j+1,nf)
call HLPA(Con_s,Fww,Fw,Fp,Fe,Delta_f)
Sp(i,j) = Sp(i,j) + Con_s * Delta_f
!------------------ n face--------------------
Fww = F(i ,j-1,nf)
Fw = F(i ,j ,nf)
Fp = F(i ,j+1,nf)
Fe = F(i ,j+2,nf)
call HLPA(Con_n,Fww,Fw,Fp,Fe,Delta_f)
Sp(i,j) = Sp(i,j) + Con_n * Delta_f *(-1.)
end if
600 continue
!-------------------------------- HLPA SCHEME----------------------------
!------------------------------------------------------------------------
100 continue
!----------------------------------------------------------------
NImax = NXmaxp
NJmax = NYmaxp
F_out = CheckFlux
Filename ='0_Flx.txt'
Call Out_array(F_out,NImax,NJmax,Filename)
!-------------------------------------------------------------------
Return
End
