c
C File name ............. GXHOCS.FOR ....................... 110817
C File includes the following subroutines and functions:
C GXHOCS, GXHOEN, GXHOHL, BLKSLD, GXFLPS
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
SUBROUTINE GXHOCS
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C SUBROUTINE GXHOCS introduces
c
c higher-order convection schemes for
C use with the staggered-mesh option in PHOENICS. These are
C discretization schemes for the convection terms which have
C truncation errors of order greater than one (i.e. higher than
C First-Order Upwind). These schemes allow more accurate resolution
C of convection-dominated flows, especially when these contain
C gradients of the flow variables which are not aligned with the
C grid directions (e.g. recirculating flows).
C
C Two classes of schemes have been implemented:
C
C LINEAR SCHEMES: - These include the Central-Differencing scheme
C (CDS), the Linear-Upwind scheme (LUS), the Quadratic-Upwind
C scheme (QUICK), the Cubic-Upwind scheme (CUS) and Fromm's scheme.
C These offer generally good resolution, but do not guarantee
C boundedness and may therefore be prone to spurious oscillations
C such as "wiggles" around shock waves or negative values of KE or
C EP.
C
C NON-LINEAR FLUX-LIMITED SCHEMES: - These schemes also offer
C improved resolution but are designed to limit the higher-order
C correction so as to guarantee positive coefficients and therefore
C bounded behaviour. This limiting may reduce the accuracy of the
C discretization to some extent. These schemes modify themselves
C according to the local flow conditions and are therefore also
C termed "non-linear". The flux limiters that have been implemented
C are:
C SMART (bounded QUICK)
C Koren (bounded CUS, TVD)
C MUSCL (van Leer: bounded Fromm, TVD)
C HQUICK (based on QUICK, smooth/harmonic)
C OSPRE (smooth/continuous)
C van Leer harmonic (smooth/harmonic, TVD)
C van Albada (smooth/continuous)
C Minmod (= Zhu/Rodi SOUCUP scheme, TVD)
C Superbee (compressive limiter for sharp gradients, TVD)
C UMIST (bounded QUICK, TVD )
C HCUS (based on CUS, smooth/harmonic)
C CHARM (based on QUICK and UDS, polymonial ratio)
C
C These two classes, as implemented, offer a choice of 17
C different schemes, though not all of these are in fact recom-
C mended for general use. Those recommended are:
C
C (i) Linear: either CUS or QUICK. LUS is less accurate, but much
C more stable.
C
C (ii) Non-linear: SMART and Koren give good accuracy with strong
C relaxation. MUSCL and H-QUICK are less accurate but more
C stable. OSPRE is the most stable but again a little less
C accurate. van Leer harmonic, van Albada, Minmod and Superbee
C are classical limiters that are included for completeness.
C
C GXHOCS is called from Group 1, Sections 1 and 2, and also from
C Group 13 of GREX3:
C
C ISC=13 Section 13 : COVAL(HOCS,PHI,FIXFLU,GRND1)
C Linear - KAPPA schemes
C
C ISC=14 Section 14 : COVAL(HOCS,PHI,FIXFLU,GRND2)
C Non-Linear - Flux-Limited Positive schemes
C
C The integer array INLCS(1) in GXHOCS contains the scheme number
C for each variable PHI. The LINEAR schemes available are:
C INLCS(phi) = 1 -> LUS = 2 -> FROMM = 3 -> CUS
C = 4 -> QUICK = 5 -> CDS
C The NON-LINEAR schemes available are:
C INLCS(phi) = 6 -> SMART = 7 -> KOREN = 8 -> MUSCL
C = 9 -> H-QUICK = 10 -> OSPRE
C = 11 -> VAN LEER HARMONIC/HLPA
C = 12 -> VAN ALBADA = 13 -> MINMOD/SOUCUP
C = 14 -> SUPERBEE = 15 -> UMIST
C = 16 -> H-CUS = 17 -> CHARM
C
C Q1-file activation:
C
C A particular scheme is assigned to one or more dependent variables
C by the SCHEME command, whose syntax is:
C
C SCHEME(NAME,variable name 1,variable name 2,...etc.)
C
C The first argument NAME identifies the required discretisation
C scheme, as follows:
C NAME = LUS, FROMM, CUS, QUICK, CDS, SMART, KOREN, MUSCL, HQUICK,
C OSPRE, VANLH, VANALB, MINMOD, SUPBEE, UMIST, HCUS or CHARM.
C
C The second argument permits the user to specify those SOLVEd
C variables which will use the selected scheme. If ALL is entered as
C the second argument, then the selected scheme is applied to all
C SOLVEd-for variables. For example,
C
C SCHEME(QUICK,U1,V1);SCHEME(SMART,H1,C1,C2)
C
C selects QUICK for U1 and V1, and SMART for H1,C1 and C2, and
C Upstream differencing (as DIFCUT=0) for any SOLVEd variables which
C do not appear in a SCHEME command. The INLCS array is determined
C by the SCHEME command, and then transmitted to GXHOCS via the
C SPEDAT facility. The schemes may also be activated via the MENU.
C
C The schemes are described under the Encyclopaedia entry 'SCHEMES
C FOR CONVECTION DISCRETISATION' ( see also the HELP entry for the PIL
C command SCHEME ).
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'satear'
COMMON /IDA10/ INLCS(150)
COMMON /LB/ LB1(16),LZDZ,LB2(8),LZDZG,LB3(3),DXU2D,DYV2D,
1 LB4(8),DXG2D,DYG2D,LB5(96),LRHO1,LRHO1H,LB6,
1 LRHO2,LRHO2H,LB7(3),AEAST,ANORTH,AHIGH,LB8(46),
1 LCEA,LCNA,LCHA,LCEA2,LCNA2,LCHA2,LB9(58),CON1E,
1 CON2E,CON1N,CON2N,CON1H,CON2H,LB10(36)
1 /GENI/ NXNY,IGN1(8),NFM,IGN2(38),NPHI4,IGN3(11)
1 /IGE/ IXF,IXL,IYF,IYL,IGE1(2),IGR,ISC,IRUN,IZ,
1 IGE2,ISWEEP,ISTEP,INDVAR,VAL,IGE4(10)
1 /CGE/ NPATCH
1 /BFCIN1/ L0B(106),LULAST
1 /BFCINT/ APROJE,APROJN,APROJH,IBFC1(7),DHXPE,IBFC2,DHYPN,
1 IBFC3,DHZPH,IBFC4(12),DZGPH,IBFC5(39),XP,YP,ZP,
1 IBFC6(36)
1 /GHOCS/ L0PHI,L0VAL,L0RUA,L0RHO,L0RHOH,L0R,L0FPR,L0PRP,
1 L0CVL,L0UCT,L0VCT,L0WCT,L0XP,L0YP,L0ZP,L0DXG,
1 L0DYG,L0DZG,KAPP(150),KAPM(150),U1ORU2,V1ORV2,
1 W1ORW2,UVW,STACON,ISCHEM,DEBHOC,debhcs
1 /IZLIM/ IZF,IZL
1 /RSG/ RGFIL1(19),RSG20,RGFIL2(180)
1 /ISG/ ISGF1(57),ISCHM1,ISGF2(42)
COMMON/LMAIN2/ DBFIL1,DBGSWP,DBGTZ,DBFIL2(8)
COMMON/BFICRT/ ICRT(6),IUCMP(2),IVCMP(2),IWCMP(2)
COMMON /NAMFN/ NAMFUN,NAMSUB
COMMON/MDFCNV/LIMCNV ! convection was modified by in-form
LOGICAL LIMCNV
C
CHARACTER NPATCH*8,NAMFUN*6,NAMSUB*6,CTEMP*8
logical debhcs,debhoc,dbgswp,dbgtz,dbfil1,dbfil2
LOGICAL INIRUA,U1ORU2,V1ORV2,W1ORW2,UVW,FLUID1,STACON,SOLVE
INTEGER VAL,AEAST,ANORTH,AHIGH,APROJN,
1 APROJH,DXU2D,DYV2D,DXG2D,DYG2D,DHXPE,DHYPN,DHZPH,DZGPH,
1 CON1E,CON1N,CON1H,CON2E,CON2N,CON2H,APROJE,XP,YP,ZP
REAL KAPP,KAPM,KAPPA
CHARACTER*6 SCNAM(17)
DATA SCNAM /'LUS', 'FROMM', 'CUS', 'QUICK',
1 'CDS', 'SMART', 'KOREN', 'VANL1',
1 'HQUICK', 'OSPRE', 'VANL2', 'VANALB',
1 'MINMOD', 'SUPBEE', 'UMIST', 'HCUS',
1 'CHARM'/
SAVE INIRUA
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
NAMSUB='GXHOCS'
C
IF(IGR==1 .AND. ISC==1) THEN
C
C>>>>>>>>>>>>>>>>>>>>> Group 1 ------- Section 1 <<<<<<<<<<<<<<<<<<<<<
C
C.... Construct INLCS array from SPEDAT facility & then set parameter
C values for KAPPA schemes.
if(.not.nullpr) WRITE(LUPR1,101)
101 FORMAT(2X,'NO',3X,'VARIABLE',4X,'SCHEME')
ONED3=1./3.
DO MPHI=1,NPHI
INLCS(MPHI)=0
IF(MPHI>2 .AND. SOLVE(MPHI)) THEN
WRITE(CTEMP,'(I3.3)') MPHI
CALL GETSDI('SCHEME','INLCS'//CTEMP,INLCS(MPHI))
IF(INLCS(MPHI)/=0) THEN
IF(.NOT.NULLPR) WRITE(LUPR1,102)
1 MPHI,NAME(MPHI),SCNAM(INLCS(MPHI))
ELSE
CYCLE
ENDIF
IF(INLCS(MPHI)<6.OR.INLCS(MPHI)==10.OR.
1 INLCS(MPHI)==12) THEN
IF(.NOT.NULLPR) THEN
CALL WRIT40('Scheme '//SCNAM(INLCS(MPHI)))
CALL WRIT40('Use of this scheme is not recommended')
CALL WRIT40('See PHOENICS Encyclopaedia entry on ')
CALL WRIT40('"Schemes of discretization" ')
ENDIF
ENDIF
IF(INLCS(MPHI)==1) THEN
KAPPA=-1.
ELSEIF(INLCS(MPHI)==2) THEN
KAPPA=0.
ELSEIF(INLCS(MPHI)==3) THEN
KAPPA=ONED3
ELSEIF(INLCS(MPHI)==4) THEN
KAPPA=0.5
ELSEIF(INLCS(MPHI)==5) THEN
KAPPA=1.
ENDIF
IF(INLCS(MPHI)>=1.AND.INLCS(MPHI)<=5) THEN
KAPP(MPHI) = 0.5*(1.0 + KAPPA)
KAPM(MPHI) = 0.5*(1.0 - KAPPA)
ELSE
KAPP(MPHI) = 0.
KAPM(MPHI) = 0.
ENDIF
ENDIF
ENDDO
IF(ISCHM1>1) THEN
WRITE(LUPR1,'('' Schemes become active after sweep '',I6)')
1 ISCHM1
ELSE
WRITE(LUPR1,'('' Schemes are active for all sweeps '')')
ENDIF
CALL WRITST
102 FORMAT(1X,I3,5X,A4,5X,A6)
C
C.... STACON=T signifies that in-line convections for velocity
C will be computed by averaging scalar convection fluxes.
STACON=UCONV.AND.NAMGRD=='CONV'
CALL GETSDL('SCHEME','DEBHOC',DEBHOC)
ELSEIF(IGR==13) THEN
C
C>>>>>>>>>>>>>>>>>>>>> Group 13 ------- Section 13-14 <<<<<<<<<<<<<<<<
C
IF(NPATCH(1:4)/='HOCS') RETURN
C
ISCHEM=INLCS(INDVAR)
STACON = STACON.OR.LIMCNV
debhcs=debhoc.and.dbgswp.and.dbgtz.and.dbgphi(indvar)
C.... Initialize values of source term to zero on slab
CALL FN1(VAL,0.0)
IF(STEADY.AND.FIINIT(1)/=READFI) THEN ! always active for transient & restart
IF(ISWEEP1) L0RHOH = L0F(LRHO1H)
IF(.NOT.ONEPHS) L0R = L0F(9)
ELSE
L0RHO = L0F(DEN2)
IF(NZ>1) L0RHOH = L0F(LRHO2H)
IF(.NOT.ONEPHS) L0R = L0F(10)
ENDIF
C.... Set cell-geometry addresses
IF(NX>1) THEN
L0DXU=L0F(DXU2D); L0DXG=L0F(DXG2D); L0AE=L0F(AEAST)
IF(BFC) L0DXU = L0B(DHXPE)+IZNXNY
ENDIF
IF(NY>1) THEN
L0DYV=L0F(DYV2D); L0DYG=L0F(DYG2D); L0AN=L0F(ANORTH)
IF(BFC) L0DYV = L0B(DHYPN)+IZNXNY
ENDIF
IF(NZ>1) THEN
IF(.NOT.BFC) THEN
L0AH = L0F(AHIGH)
L0DZW=L0F(LZDZ); L0DZG=L0F(LZDZG)
ELSE
L0AH = L0B(APROJH)+IZNXNY
L0DZG = L0B(DZGPH)+IZNXNY
L0DZW = L0B(DHZPH)+IZNXNY
ENDIF
ENDIF
C
C.... N.B. Distance of cell centre from cell face in covariant
C components is not available: DHXPE etc. are NORMAL to cell faces.
C (DHXPE,DHYPN,DHZPH are accessed here but not used.)
C Hence in BFC case it is assumed faces are midway between centres.
IF(BFC) THEN
L0DXU=L0DXG; L0DYV=L0DYG; L0DZW=L0DZG
ENDIF
C
C.... Set offset indices for 1st/2nd-phase velocities and mass fluxes
IF(FLUID1(INDVAR)) THEN
IADD = 0; IAD2 = 0
ELSE
IADD = 1; IAD2 = 6
ENDIF
C.... Store original PATCH limits
C (Limits IZF,IZL are not defined and hence GETPTC required)
CALL GETPTC(NPATCH,PATYP,IXF,IXL,IYF,IYL,IZF,IZL,ITF,ITL)
IXFP=IXF; IXLP=IXL; IYFP=IYF; IYLP=IYL; IZFP=IZF; IZLP=IZL
C
C.... Check for momentum variables and reset patch limits appropriately
U1ORU2=.FALSE.
V1ORV2=.FALSE.
W1ORW2=.FALSE.
IF(INDVAR==3.OR.INDVAR==4) THEN
U1ORU2=.TRUE.
IXL = MIN(IXL,(NX-1))
ELSEIF(INDVAR==5.OR.INDVAR==6) THEN
V1ORV2=.TRUE.
IYL = MIN(IYL,(NY-1))
ELSEIF(INDVAR==7.OR.INDVAR==8) THEN
W1ORW2=.TRUE.
IZL = MIN(IZL,(NZ-1))
ENDIF
UVW = U1ORU2.OR.V1ORV2.OR.W1ORW2
IF(BFC.AND.UVW) THEN
IF(FLUID1(INDVAR)) THEN
IF(NX>1) JUCMP = IUCMP(1)
IF(NY>1) JVCMP = IVCMP(1)
IF(NZ>1) JWCMP = IWCMP(1)
IUCRT=ICRT(1); IVCRT=ICRT(3); IWCRT=ICRT(5)
ELSE
IF(NX>1) JUCMP = IUCMP(2)
IF(NY>1) JVCMP = IVCMP(2)
IF(NZ>1) JWCMP = IWCMP(2)
IUCRT=ICRT(2); IVCRT=ICRT(4); IWCRT=ICRT(6)
ENDIF
L0UCT=L0F(IUCRT); L0VCT=L0F(IVCRT); L0WCT=L0F(IWCRT)
L0XP=L0B(XP); L0YP=L0B(YP); L0ZP=L0B(ZP)
L0XP = L0XP + IZNXNY
L0YP = L0YP + IZNXNY
L0ZP = L0ZP + IZNXNY
ENDIF
C
INIRUA=(STEADY.AND.ISWEEP==FSWEEP).OR.(.NOT.STEADY.AND.
1 ISTEP==FSTEP.AND.ISWEEP==FSWEEP)
C... NORTH cell-face contribution to HOCS source term
IF(NY>1) THEN
IF(.NOT.(BFC.AND.UVW)) THEN
L0CVL = L0F(5+IADD)
ELSE
L0CVL = L0F(JVCMP)
ENDIF
IF(FLUID1(INDVAR)) THEN
LCONN = ITWO(LCNA,CON1N,NZ==1)
L0RUA = L0F(LCONN) + IZNXNY
ELSE
LCONN = ITWO(LCNA2,CON2N,NZ==1)
L0RUA = L0F(LCONN) + IZNXNY
ENDIF
IF(INIRUA.AND.NZ>1.AND.W1ORW2) THEN
DO J=1,NXNY
F(L0RUA+NXNY+J)=0.0
ENDDO
ENDIF
IYF=IYF+1; IYL=IYL-2
CALL GXHOEN(ISC,L0DYG,L0DYV,L0AN,L0DXG,V1ORV2,U1ORU2,1,NY,1)
IYF=IYFP; IYL=IYLP
ENDIF
C
C... HIGH/LOW cell-face contribution to HOCS source term
IF(NZ>1) THEN
IF(BFC.AND.W1ORW2) THEN
L0CVL = L0F(JWCMP)
ELSE
L0CVL = L0F(7+IADD)
ENDIF
C.... HIGH face
C ---------
IF(FLUID1(INDVAR)) THEN
L0RUA = L0F(CON1H) + IZNXNY
ELSE
L0RUA = L0F(CON2H) + IZNXNY
ENDIF
IF(INIRUA.AND.NZ>1.AND.W1ORW2) THEN
DO J=1,NXNY
F(L0RUA+NXNY+J)=0.0
ENDDO
ENDIF
IF(IZ>=IZF+1 .AND. IZ<=IZL-2) THEN
CALL GXHOHL(5,ISC,L0DZW,L0DZG,L0AH)
ENDIF
C.... LOW face
C --------
IF(FLUID1(INDVAR)) THEN
L0RUA = L0F(CON1H) + (IZ-2)*NXNY
ELSE
L0RUA = L0F(CON2H) + (IZ-2)*NXNY
ENDIF
c L0CVL = L0CVL-NFM
c L0PHI = L0PHI-NFM
c IF(BFC)
c 1 CALL SUB3(L0UCT,L0UCT-NFM,L0VCT,L0VCT-NFM,L0WCT,L0WCT-NFM)
IF(IZ>=IZF+2 .AND. IZ<=IZL-1) THEN
L0CVL = L0CVL-NFM
L0PHI = L0PHI-NFM
IF(BFC) THEN
L0UCT=L0UCT-NFM; L0VCT=L0VCT-NFM; L0WCT=L0WCT-NFM
ENDIF
CALL GXHOHL(6,ISC,L0DZW,L0DZG,L0AH)
L0CVL = L0CVL+NFM
L0PHI = L0PHI+NFM
IF(BFC) THEN
L0UCT=L0UCT+NFM; L0VCT=L0VCT+NFM; L0WCT=L0WCT+NFM
ENDIF
ENDIF
IZF=IZFP; IZL=IZLP
ENDIF
C
C... EAST cell-face contribution to HOCS source term
IF(NX>1) THEN
IF(.NOT.(BFC.AND.UVW)) THEN
L0CVL = L0F(3+IADD)
ELSE
L0CVL = L0F(JUCMP)
ENDIF
IF(FLUID1(INDVAR)) THEN
LCONE = ITWO(LCEA,CON1E,NZ==1)
L0RUA = L0F(LCONE) + IZNXNY
ELSE
LCONE = ITWO(LCEA2,CON2E,NZ==1)
L0RUA = L0F(LCONE) + IZNXNY
ENDIF
IF(INIRUA.AND.NZ>1.AND.W1ORW2) THEN
DO J=1,NXNY
F(L0RUA+NXNY+J)=0.0
ENDDO
ENDIF
IXF=IXF+1; IXL=IXL-2
CALL GXHOEN(ISC,L0DXG,L0DXU,L0AE,L0DYG,U1ORU2,V1ORV2,NY,1,3)
IXF=IXFP; IXL=IXLP
ENDIF
ENDIF
NAMSUB='gxhocs'
END
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXHOEN(ISC,L0DG,L0DV,L0AF,L0DGN,UVI,UVN,INDI,INDN,IDIR)
C *********************************************************************
C Purpose: to calculate higher-order correction for EAST or NORTH
C faces on a single slab of cells
C east face north face
C Input: L0DG = L0DXG L0DYG
C L0DV = L0DXU L0DYV
C L0AF = L0AE L0AN
C L0DGN = L0DYG L0DXG
C UVI = U1ORU2 V1ORV2
C UVN = V1ORV2 U1ORU2
C INDI = NY 1
C INDN = 1 NY
C IDIR = 3 1
C Called by: GXHOCS
C *********************************************************************
INCLUDE 'farray'
COMMON
1 /IDATA/ NX,NY,NZ,LUPR1,IDAT1(116)
1 /LDATA/ CARTES,LD1(2),ONEPHS,LD2(3),XCYCLE,LD3(10),STEADY,
1 BFC,LD4(44),NONORT,LD5(19)
1 /IGE/ IXF,IXL,IYF,IYL,IGE1(5),IZ,
1 IGE3,ISWEEP,IGE4,INDVAR,IGE5(11)
1 /GENI/ NXNY,IGN1(8),NFM,IGN2(38),NPHI4,IGN3(11)
1 /GHOCS/ L0PHI,L0VAL,L0RUA,L0RHO,L0RHOH,L0R,L0FPR,L0PRP,
1 L0CVL,L0UCT,L0VCT,L0WCT,L0XP,L0YP,L0ZP,LBDXG,
1 LBDYG,LBDZG,KAPP(150),KAPM(150),U1ORU2,V1ORV2,
1 W1ORW2,UVW,STACON,ISCHEM,DEBHOC,debhcs
1 /NAMFN/ NAMFUN,NAMSUB
1 /RDATA/ TINY,RDAT(84)
C
logical debhoc,debhcs
LOGICAL BLKSLD,DOCOR,U1ORU2,V1ORV2,W1ORW2,CARTES,UVI,UVN,UVW,LD1,
1 XCYCLE,LD2,STEADY,BFC,LD3,LD4,LD5,NONORT,QGE,ONEPHS,
1 STACON
REAL KAPP,KAPM
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXHOEN'
C
C.... Cycle over the patch.
debhcs=.false.
if(debhcs) call banner(1,namsub,110817)
IZNXNY=(IZ-1)*NXNY
DO 10 IX=IXF,IXL
IADD = (IX-1)*NY
DO 10 IY=IYF,IYL
if(debhcs) call writ2i('ix ',ix,'iy ',iy)
IXY = IY+IADD
C.... Check for blocked faces or solid cells within stencil
DOCOR=.NOT.BLKSLD(UVI,UVN,W1ORW2,INDI,INDN,NXNY,IXY+IZNXNY,
1 IDIR)
if(debhcs) call writ1l('docor ',docor)
C.... If no blocked/solid cells/faces within stencil then continue
IF(DOCOR) THEN
C
C.... Set mass-flux coefficient for face
RUA = F(L0RUA+IXY)
if(debhcs) call writ3l('uvi ',uvi,'uvn ',uvn,
1 'stacon ',stacon)
IF(UVI) THEN
IF(STACON) THEN
RUA = 0.5*(RUA + F(L0RUA+IXY+INDI))
ELSE
IF(F(L0CVL+IXY)>0.0) THEN
GRO = F(L0RHO+IXY)
IF(.NOT.ONEPHS) GVF = F(L0R+IXY)
ELSE
GRO = F(L0RHO+IXY+INDI)
IF(.NOT.ONEPHS) GVF = F(L0R+IXY+INDI)
ENDIF
GUF = 0.5*(F(L0CVL+IXY) + F(L0CVL+IXY+INDI))
RUA = GRO*GUF*F(L0AF+IXY)
IF(.NOT.ONEPHS) RUA=RUA*GVF
ENDIF
ELSEIF(UVN) THEN
RUA = 0.5*(RUA + F(L0RUA+IXY+INDN))
ELSEIF(W1ORW2) THEN
C N109 fails because f(l0rua+ixy+nxny) has 1e19 on 1st sweep
c see mrm 06.05.05 initialise con1n & con1e for 1st sweep
RUA = 0.5*(RUA + F(L0RUA+IXY+NXNY))
ENDIF
C
C.... Values of convected variable on each side of face
FP = F(L0PHI+IXY)
FEN = F(L0PHI+IXY+INDI)
C
C.... Calculate required cell-geometry data
IF(UVI) THEN
DXC = F(L0DV+IXY+INDI)
if(dxc<=0.0) then
call writ3i('indvar ',indvar,'ix ',ix,
1 'iy ',iy)
CALL SET_ERR(444, 'Error. See reslt file.',1)
C stop
endif
IF(QGE(RUA,0.0)) THEN
DXU = F(L0DV+IXY)
GDX = 0.5*F(L0DV+IXY+INDI)
ELSE
DXU = F(L0DV+IXY+2*INDI)
GDX = -0.5*F(L0DV+IXY+INDI)
ENDIF
ELSE
DXC = F(L0DG+IXY)
IF(QGE(RUA,0.0)) THEN
DXU = F(L0DG+IXY-INDI)
GDX = 0.5*F(L0DV+IXY)
ELSE
DXU = F(L0DG+IXY+INDI)
GDX = -0.5*F(L0DV+IXY+INDI)
ENDIF
ENDIF
C
C.... Calculate gradients of convected variable
IF(.NOT.(BFC.AND.NONORT.AND.UVW)) THEN
C
GRPHC = (FEN - FP)/(DXC+tiny)
IF(QGE(RUA,0.0)) THEN
fpnei = F(L0PHI+IXY-INDI)
GRPHU = (FP - F(L0PHI+IXY-INDI))/(DXU+tiny)
ELSE
fpnei = F(L0PHI+IXY+2*INDI)
GRPHU = (F(L0PHI+IXY+2*INDI) - FEN)/(DXU+tiny)
ENDIF
C
ELSE
C
IF(UVI) THEN
INDB = INDI
DP = F(L0DG+IXY)
DE = F(L0DG+IXY+INDI)
if(debhcs) call writ2r('dp 1',dp,'de ',de)
ELSEIF(UVN) THEN
INDB = INDN
DP = F(L0DGN+IXY)
DE = F(L0DGN+IXY+INDI)
if(debhcs) then
call writ2r('dp 2',dp,'de ',de)
call writ3i('l0dgn ',l0dgn,'ixy ',ixy,
1 'indi ',indi)
endif
ELSE
INDB = NXNY
DP = F(LBDZG+IXY)
DE = F(LBDZG+IXY+INDI)
if(debhcs) call writ2r('dp 3',dp,'de ',de)
ENDIF
C
if(de<=0.0.or.dp<=0) then
call writ3i('indvar ',indvar,'ix ',ix,
1 'iy ',iy)
call writ2r('dp ',dp,'de ',de)
CALL SET_ERR(511, 'Error. See result file.',1)
C call wayout(1)
endif
DP=DP+TINY
DE=DE+TINY
DPX = (F(L0XP+IXY+INDB)-F(L0XP+IXY))/DP
DPY = (F(L0YP+IXY+INDB)-F(L0YP+IXY))/DP
DPZ = (F(L0ZP+IXY+INDB)-F(L0ZP+IXY))/DP
C
DEX = (F(L0XP+IXY+INDI+INDB)-F(L0XP+IXY+INDI))/DE
DEY = (F(L0YP+IXY+INDI+INDB)-F(L0YP+IXY+INDI))/DE
DEZ = (F(L0ZP+IXY+INDI+INDB)-F(L0ZP+IXY+INDI))/DE
C
FENP = F(L0UCT+IXY+INDI)*DPX+F(L0VCT+IXY+INDI)*DPY+
1 F(L0WCT+IXY+INDI)*DPZ
C
GRPHC = (FENP - FP)/(DXC+tiny)
IF(QGE(RUA,0.0)) THEN
FUPW = F(L0UCT+IXY-INDI)*DPX+F(L0VCT+IXY-INDI)*DPY+
1 F(L0WCT+IXY-INDI)*DPZ
GRPHU = (FP - FUPW)/(DXU+tiny)
ELSE
FUPW = F(L0UCT+IXY+2*INDI)*DPX+F(L0VCT+IXY+2*INDI)*DPY+
1 F(L0WCT+IXY+2*INDI)*DPZ
GRPHU = (FUPW - FENP)/(DXU+tiny)
ENDIF
C
FPE = F(L0UCT+IXY)*DEX+F(L0VCT+IXY)*DEY+F(L0WCT+IXY)*DEZ
C
GRPHCI = (FEN - FPE)/(DXC+tiny)
IF(QGE(RUA,0.0)) THEN
FUPW = F(L0UCT+IXY-INDI)*DEX+F(L0VCT+IXY-INDI)*DEY+
1 F(L0WCT+IXY-INDI)*DEZ
GRPHUI = (FPE - FUPW)/(DXU+tiny)
ELSE
FUPW = F(L0UCT+IXY+2*INDI)*DEX+F(L0VCT+IXY+2*INDI)*DEY+
1 F(L0WCT+IXY+2*INDI)*DEZ
GRPHUI = (FUPW - FEN)/(DXU+tiny)
ENDIF
C
ENDIF
C
C.... Use gradients to calculate higher-order correction
IF(ISC==13) THEN
C.... Kappa schemes
FLUXO = RUA*(KAPP(INDVAR)*GRPHC+KAPM(INDVAR)*GRPHU)*GDX
FLUXI = FLUXO
IF(BFC.AND.NONORT.AND.UVW) FLUXI=RUA*
1 (KAPP(INDVAR)*GRPHCI+KAPM(INDVAR)*GRPHUI)*GDX
ELSE
C.... Flux-limited positive schemes
CALL GXFLPS(GRPHC,GRPHU,ISCHEM,PSI)
FLUXO = RUA*PSI*GRPHU*GDX
FLUXI = FLUXO
IF(BFC.AND.NONORT.AND.UVW) THEN
CALL GXFLPS(GRPHCI,GRPHUI,ISCHEM,PSI)
FLUXI = RUA*PSI*GRPHUI*GDX
ENDIF
ENDIF
C
F(L0VAL+IXY) = F(L0VAL+IXY) - FLUXO
F(L0VAL+IXY+INDI) = F(L0VAL+IXY+INDI) + FLUXI
ENDIF
10 CONTINUE
NAMSUB = 'gxhoen'
if(debhcs) call banner(2,namsub,110817)
END
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXHOHL(NDIR,ISC,L0DZW,L0DZG,L0AH)
C *********************************************************************
C Purpose: to calculate higher-order correction for HIGH/LOW faces
C on a single slab of cells
C low face high face
C Input: NDIR = 6 5
C Called by: GXHOCS
C *********************************************************************
INCLUDE 'farray'
COMMON
1 /IDATA/ NX,NY,NZ,LUPR1,IDAT1(116)
1 /LDATA/ CARTES,LD1(2),ONEPHS,LD2(3),XCYCLE,LD3(10),STEADY,
1 BFC,LD4(44),NONORT,LD5(19)
1 /IGE/ IXF,IXL,IYF,IYL,IGE1(2),IGR,ISC1,IRUN,IZ,
1 IGE2,ISWEEP,IGE3,INDVAR,IGE4(11)
1 /IZLIM/ IZF,IZL
1 /GENI/ NXNY,IGN1(8),NFM,IGN2(38),NPHI4,IGN3(11)
1 /GHOCS/ L0PHI,L0VAL,L0RUA,L0RHO,L0RHOH,L0R,L0FPR,L0PRP,
1 L0CVL,L0UCT,L0VCT,L0WCT,L0XP,L0YP,L0ZP,LBDXG,
1 LBDYG,LBDZG,KAPP(150),KAPM(150),U1ORU2,V1ORV2,
1 W1ORW2,UVW,STACON,ISCHEM,debhoc,debhcs
1 /NAMFN/ NAMFUN,NAMSUB
1 /RDATA/ TINY,RDAT(84)
C
LOGICAL BLKSLD,DOCOR,U1ORU2,V1ORV2,W1ORW2,CARTES,UVW,LD1,XCYCLE,
1 LD2,STEADY,BFC,LD3,NONORT,LD4,QGE,LD5,ONEPHS,STACON
logical debhoc,debhcs
REAL KAPP,KAPM
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXHOHL'
C
C.... Set indices for cell geometry
IF(.NOT.BFC) THEN
LDVXY = L0DZW+IZ
LDGXY = L0DZG+IZ
ELSE
LDZG = L0DZG
LDZ = L0DZW
ENDIF
C
C.... Adjust indices by one slab if LOW face
ISLAB=(IZ-1)*NXNY
IF(NDIR==6) THEN
ISLAB=(IZ-2)*NXNY
IF(.NOT.BFC) THEN
LDVXY = LDVXY-1
LDGXY = LDGXY-1
ELSE
LDZG = LDZG-NXNY
LDZ = LDZ-NXNY
L0AH = L0AH-NXNY
ENDIF
ENDIF
C
INDI = NFM
INDZ = 1
IF(BFC) INDZ = NXNY
C
IF(U1ORU2) THEN
IXLP= IXL
IXL = MIN(IXL,(NX-1))
ELSEIF(V1ORV2) THEN
IYLP= IYL
IYL = MIN(IYL,(NY-1))
ENDIF
C.... Cycle over the patch.
DO 10 IX=IXF,IXL
IADD = (IX-1)*NY
DO 10 IY=IYF,IYL
IXY = IY+IADD
C
C.... Check for blocked faces or solid cells within stencil
DOCOR=.NOT.BLKSLD(W1ORW2,U1ORU2,V1ORV2,NXNY,NY,1,IXY+ISLAB,5)
C.... If no blocked/solid cells/faces within stencil then continue
IF(DOCOR) THEN
C
C.... Set mass-flux coefficient for face
RUA = F(L0RUA+IXY)
IF(W1ORW2) THEN
IF(STACON) THEN
RUA = 0.5*(RUA + F(L0RUA+IXY+NXNY))
ELSE
IF(F(L0CVL+IXY)>0.0) THEN
GRO = F(L0RHO+IXY)
IF(.NOT.ONEPHS) GVF = F(L0R+IXY)
ELSE
GRO = F(L0RHOH+IXY)
IF(.NOT.ONEPHS) GVF = F(L0R+IXY+INDI)
ENDIF
GUF = 0.5*(F(L0CVL+IXY) + F(L0CVL+IXY+INDI))
GAF = F(L0AH+IXY)
RUA = GRO*GUF*GAF
IF(.NOT.ONEPHS) RUA=RUA*GVF
ENDIF
ELSEIF(U1ORU2) THEN
RUA = 0.5*(RUA + F(L0RUA+IXY+NY))
ELSEIF(V1ORV2) THEN
RUA = 0.5*(RUA + F(L0RUA+IXY+1))
ENDIF
C
C.... Values of convected variable on each side of face
FP = F(L0PHI+IXY)
FH = F(L0PHI+IXY+INDI)
C
C.... Calculate required cell-geometry data
IF(BFC) THEN
LDGXY = LDZG+IXY
LDVXY = LDZ +IXY
ENDIF
IF(W1ORW2) THEN
DZC = F(LDVXY+INDZ)
IF(QGE(RUA,0.0)) THEN
DZU = F(LDVXY)
GDZ = 0.5*F(LDVXY+INDZ)
ELSE
DZU = F(LDVXY+2*INDZ)
GDZ = -0.5*F(LDVXY+INDZ)
ENDIF
ELSE
DZC = F(LDGXY)
IF(QGE(RUA,0.0)) THEN
DZU = F(LDGXY-INDZ)
GDZ = 0.5*F(LDVXY)
ELSE
DZU = F(LDGXY+INDZ)
GDZ = -0.5*F(LDVXY+INDZ)
ENDIF
ENDIF
C
C.... Calculate gradients of convected variable
IF(.NOT.(BFC.AND.NONORT.AND.UVW)) THEN
C
dzc=dzc+tiny
dzu=dzu+tiny
GRPHC = (FH-FP)/DZC
IF(QGE(RUA,0.0)) THEN
fpnei=F(L0PHI+IXY-INDI)
GRPHU = (FP-F(L0PHI+IXY-INDI))/DZU
ELSE
fpnei= F(L0PHI+IXY+2*INDI)
GRPHU = (F(L0PHI+IXY+2*INDI)-FH)/DZU
ENDIF
C
ELSE
C
IF(W1ORW2) THEN
INDB = NXNY
DP = F(LBDZG+IXY)
ELSEIF(U1ORU2) THEN
INDB = NY
DP = F(LBDXG+IXY)
ELSEIF(V1ORV2) THEN
INDB = 1
DP = F(LBDYG+IXY)
ENDIF
C
dp=dp+tiny
DPX = (F(L0XP+IXY+INDB)-F(L0XP+IXY))/DP
DPY = (F(L0YP+IXY+INDB)-F(L0YP+IXY))/DP
DPZ = (F(L0ZP+IXY+INDB)-F(L0ZP+IXY))/DP
C
IF(NDIR==5) THEN
FH = F(L0UCT+IXY+INDI)*DPX+F(L0VCT+IXY+INDI)*DPY+
1 F(L0WCT+IXY+INDI)*DPZ
ELSE
FP = F(L0UCT+IXY)*DPX+F(L0VCT+IXY)*DPY+F(L0WCT+IXY)*DPZ
ENDIF
C
GRPHC = (FH - FP)/DZC
IF(QGE(RUA,0.0)) THEN
FUPW = F(L0UCT+IXY-INDI)*DPX+F(L0VCT+IXY-INDI)*DPY+
1 F(L0WCT+IXY-INDI)*DPZ
GRPHU = (FP - FUPW)/DZU
ELSE
FUPW = F(L0UCT+IXY+2*INDI)*DPX+F(L0VCT+IXY+2*INDI)*DPY+
1 F(L0WCT+IXY+2*INDI)*DPZ
GRPHU = (FUPW - FH)/DZU
ENDIF
C
ENDIF
C
C.... Use gradients to calculate higher-order correction
IF(ISC==13) THEN
C.... Kappa schemes
FLUX = RUA*(KAPP(INDVAR)*GRPHC+KAPM(INDVAR)*GRPHU)*GDZ
ELSE
C.... Flux-limited positive schemes
CALL GXFLPS(GRPHC,GRPHU,ISCHEM,PSI)
FLUX = RUA*PSI*GRPHU*GDZ
ENDIF
C
IF(NDIR==5) THEN
F(L0VAL+IXY) = F(L0VAL+IXY) - FLUX
ELSE
F(L0VAL+IXY) = F(L0VAL+IXY) + FLUX
ENDIF
ENDIF
C
10 CONTINUE
C
IF(U1ORU2) THEN
IXL = IXLP
ELSEIF(V1ORV2) THEN
IYL = IYLP
ENDIF
C
NAMSUB = 'gxhohl'
C
END
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXFLPS(GRC,GRU,ISC,PSI)
C *********************************************************************
C Purpose: to calculate flux-limiter functions for Flux-Limited
C Positive Schemes. Ten limiters have been implemented:
C MUSCL, van Leer harmonic, van Albada, Koren, SMART,
C Minmod, OSPRE, H-Quick, Superbee and UMIST. Others may
C be implemented by imitation.
C
C Input: GRC - gradient of convected variable across face
C GRU - gradient of convected variable upwind of face
C ISC - scheme number which indicates scheme selection
C
C Output: PSI - Flux-limiter function
C Called by: GXHOEN, GXHOHL
C *********************************************************************
C
c.... initialisation in case of inappropriate entry
PSI=0.0
RP=(GRC+1.E-10)/(GRU+1.E-10)
C
IF(ISC==6) THEN
C
C..SMART: PLNS limiter based on QUICK (Gaskell & Lau [1988])
C
R1=AMIN1(0.25+0.75*RP,4.0)
R1=AMIN1(R1,2.0*RP)
PSI=AMAX1(0.0,R1)
C
ELSEIF(ISC==7) THEN
C
C..Koren: PLNS limiter based on CUS (Koren [1993])
C
R1=AMIN1(0.33333+0.66666*RP,2.0)
R1=AMIN1(R1,2.0*RP)
PSI=AMAX1(1.E-20,R1)
C
ELSEIF(ISC==8) THEN
C
C..van Leer1/MUSCL : PLS limiter based on Fromm & identical to Noll'S
C MLU scheme (van Leer[1979])
C
R1=AMIN1(0.5+0.5*RP,2.0)
R1=AMIN1(R1,2.0*RP)
PSI=AMAX1(1.E-20,R1)
C
ELSEIF(ISC==9) THEN
C
C..H-QUICK: HNS limiter based on QUICK (Deconinck & Waterson [1995])
C
IF(RP>0.0) THEN
PSI=4.0*RP/(RP+3.0)
ELSE
PSI=1.0E-20
ENDIF
C
ELSEIF(ISC==10) THEN
C
C..OSPRE: PRS limiter (Deconinck & Waterson [1995])
C
RP2=RP*RP
PSI=1.5*(RP2+RP)/(RP2+RP+1.0)
C
ELSEIF(ISC==11) THEN
C
C..van Leer2: PRS harmonic limiter based on Fromm, which is
C identical to Zhu HLPA scheme (see Hirsch[1990])
C
IF(RP>0.0) THEN
PSI=2.0*RP/(RP+1.0)
ELSE
PSI=1.0E-20
ENDIF
C
ELSEIF(ISC==12) THEN
C
C..van Albada: 'classical' limiter (van Albada et al [1982])
C
PSI=RP*(RP+1.0)/(RP*RP+1.0)
C
ELSEIF(ISC==13) THEN
C
C..Minmod: 'classical' limiter, which is identical to Zhu/Rodi SOUCUP
C scheme (see Roe [1986])
C
RP=AMIN1(RP,1.0)
PSI=AMAX1(0.0,RP)
C
ELSEIF(ISC==14) THEN
C
C..Superbee: 'classical' limiter (see Roe [1986])
C
R1=AMIN1(2*RP,1.0)
R2=AMIN1(RP,2.0)
R2=AMAX1(R1,R2)
PSI=AMAX1(0.0,R2)
C
ELSEIF(ISC==15) THEN
C
C..UMIST: PLS limiter based mainly on QUICK Lien & Leschziner [1994]
C
R1=AMIN1(0.75+0.25*RP,2.0)
R2=AMIN1(0.25+0.75*RP,2.0*RP)
R1=AMIN1(R1,R2)
PSI=AMAX1(0.0,R1)
C
ELSEIF(ISC==16) THEN
C
C..H-CUS: HNS limiter based on CUS (Deconinck & Waterson [1995])
C
IF(RP>0.0) THEN
PSI=3.0*RP/(RP+2.0)
ELSE
PSI=1.0E-20
ENDIF
C
ELSEIF(ISC==17) THEN
C
C..CHARM: PRNS limiter based on QUICK (Zhou [1995])
C
IF(RP>0.0) THEN
PSI=RP*(3.*RP+1.)/(RP*RP+2.0*RP+1.0)
ELSE
PSI=1.0E-20
ENDIF
ENDIF
C
END
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
FUNCTION BLKSLD(VELI,VELN1,VELN2,INDI,INDN1,INDN2,IJK,IDIR)
C *********************************************************************
C Purpose: to detect presence of blocked faces or solid cells
C within stencil used to calculate higher-order correction
C for a single cell face. If BLKSLD is .TRUE. then no
C higher-order correction will be calculated.
C
C Input: VELI - .TRUE. if variable is "in-line" velocity
C VELN1 - .TRUE. if variable is NOT "in-line" velocity
C VELN2 - .TRUE. if variable is NOT "in-line" velocity
C INDI - Offset index in dirn. normal to face
C INDN1 - Offset index corresponding to VELN1
C INDN2 - Offset index corresponding to VELN2
C IJK - IY+(IX-1)*NY + (IZ-1)*NXNY
C IDIR -
C
C Output: BLKSLD = .TRUE. if blockage/solid detected
C Called by: GXHOEN, GXHOHL
C *********************************************************************
C
INCLUDE 'farray'
C
LOGICAL BLKSLD,VELI,VELN1,VELN2,LSLDIR,LSOLID
C
BLKSLD = LSOLID(IJK)
IF(BLKSLD) RETURN
C.... Check for blocked faces
BLKSLD = LSLDIR(IJK,IDIR)
BLKSLD = BLKSLD.OR.LSLDIR(IJK+INDI,IDIR)
BLKSLD = BLKSLD.OR.LSLDIR(IJK-INDI,IDIR)
IF(BLKSLD) RETURN
C
IF(VELI) THEN
BLKSLD = LSLDIR(IJK+2*INDI,IDIR)
ELSEIF(VELN1) THEN
BLKSLD = LSLDIR(IJK+INDN1,IDIR)
BLKSLD = BLKSLD.OR.LSLDIR(IJK+INDI+INDN1,IDIR)
BLKSLD = BLKSLD.OR.LSLDIR(IJK-INDI+INDN1,IDIR)
ELSEIF(VELN2) THEN
BLKSLD = LSLDIR(IJK+INDN2,IDIR)
BLKSLD = BLKSLD.OR.LSLDIR(IJK+INDI+INDN2,IDIR)
BLKSLD = BLKSLD.OR.LSLDIR(IJK-INDI+INDN2,IDIR)
ENDIF
C
END
c