c
c filename GXCONVEC.HTM 050315
SUBROUTINE GXCONV
C**********************************************************************
C This subroutine applies corrections to the convective terms of the
C momentum equation. The corrections are required if shocks or
C density gradients due to other causes are present. The effect of
C these corrections is to make the convective fluxes used in the
C momentum equations consistent with those used in the continuity
C equation. It is activated when UCONV=T and NAMGRD=CONV.
C
C If the static enthalpy or static temperature equations are
C solved, the correction applied to the momentum equations means
C that total enthalpy is no longer preserved. This is corrected
C by deactivating the built-in pressure term in the H1 equation,
C and replacing it by a modifed version supplied in this routine.
C The modified source is introduced via the following Q1 statements:
C
C PATCH(ENTSOR,CELL,1,NX,1,NY,1,NZ,1,LSTEP)
C COVAL(ENTSOR,H1 [or TEM1] ,FIXFLU,GRND1)
C
C Transient cases require
C PATCH(DPDT,VOLUME,1,NX,1,NY,1,NZ,1,LSTEP)
C COVAL(DPDT, H1[or TEM1], FIXFLU, GRND1)
C
C If very steep density gradients are present, convergence may be
C improved by setting DENPCO=T, which fully accounts for density
C gradients in the pressure correction equation.
C
C The routine is called from GREX, and contains active code in
C Group 1, Section 1, and Group 8, Section 8.
C
C The present coding makes the following assumptions:
C 1. BFC cases require STORE(UCMP,VCMP,WCMP) in Q1.
C 2. Two phase BFC cases require STORE(U2CM,V2CM,W2CM)
C
C The modfications are described in detail in CHAM TR/147, and in
C the Enyclopaedia entry 'Compressible Flow Corrections'.
C
C Authors: Mike Malin, John Ludwig, CHAM UK
C Date : 24/12/91
C Modification date: 1/07/94
C**********************************************************************
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
LOGICAL ERRO, HSOR, SLD
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB,ERRMES*40
COMMON/GENI/IGNI1(9),NFM,IGNI11(50)
1 /F0/IF1(258),I1E,I2E,I1N,I2N,I1H,I2H,IF2(39)
COMMON/LSG/LSGD1(63),GCV,LSGD2(36)
LOGICAL LSGD1,GCV,LSGD2
COMMON/STOGEO/STORGEO
LOGICAL STORGEO
INTEGER HPORH
SAVE IU1,IU2,IV1,IV2,IW1,IW2,ITEM
SAVE L0CH,L0CL,L0AL,L0FE,L0FW,L0FN,L0FS,L0FH,L0FL,L0CLL
SAVE HSOR
C
NAMSUB = 'GXCONV'
C*****************************************************************
C
C---- GROUP 1. Run title and other preliminaries
C
IF(IGR==1.AND.ISC==1) THEN
IF(UCONV) THEN
C User may here change message transmitted to the VDU screen or
C batch-run log file.
IF(.NOT.NULLPR) THEN
CALL WRYT40('GXCONVEC OF 050315 HAS BEEN CALLED ')
CALL WRIT40('GXCONVEC OF 050315 HAS BEEN CALLED ')
ENDIF
C
C--- Reserve local storage
CALL GXMAKE(L0CH,NX*NY,'CH ')
CALL GXMAKE(L0CL,NX*NY,'CL ')
IF(.NOT.STORGEO.AND.NZ>1) CALL GXMAKE(L0AL,NX*NY,'ALOW')
HSOR=.FALSE.
ITEM=LBNAME('TEM1')
C--- Stores for enthalpy pressure-term correction
IF(SOLVE(H1).OR.SOLVE(ITEM)) THEN
HSOR=.TRUE.
IF(NX>1) THEN
CALL GXMAKE0(L0FE,NX*NY*NZ,'FE ')
CALL GXMAKE0(L0FW,NX*NY*NZ,'FW ')
ENDIF
IF(NY>1) THEN
CALL GXMAKE0(L0FN,NX*NY*NZ,'FN ')
CALL GXMAKE0(L0FS,NX*NY*NZ,'FS ')
ENDIF
IF(NZ>1) THEN
CALL GXMAKE0(L0FH,NX*NY*NZ,'FH ')
CALL GXMAKE0(L0FL,NX*NY*NZ,'FL ')
CALL GXMAKE(L0CLL,NX*NY,'CLL ')
CALL FN1(-(L0FL+(NZ-1)*NX*NY),1.0)
ENDIF
ENDIF
C--- Select velocity stores - use components for BFC
IF(BFC.AND..NOT.GCV) THEN
IU1=LBNAME('UCMP'); IV1=LBNAME('VCMP'); IW1=LBNAME('WCMP')
ERRO=(NX>1.AND.IU1==0).OR.(NY>1.AND.IV1==0).OR.
1 (NZ>1.AND.IW1==0)
IF(ERRO) THEN
ERRMES='STORE(UCMP,VCMP,WCMP) in Q1. '
GO TO 10
ENDIF
IF(.NOT.ONEPHS) THEN
IU2=LBNAME('U2CM'); IV2=LBNAME('V2CM'); IW2=LBNAME('W2CM')
ERRO=(NX>1.AND.IU2==0).OR.(NY>1.AND.IV2==0).OR.
1 (NZ>1.AND.IW2==0)
IF(ERRO) THEN
ERRMES='STORE(U2CM,V2CM,W2CM) in Q1. '
GO TO 10
ENDIF
ENDIF
ELSE
IU1=U1; IV1=V1; IW1=W1
IF(.NOT.ONEPHS) THEN
IU2=U2; IV2=V2; IW2=W2
ENDIF
ENDIF
RETURN
C--- Error exit
10 CONTINUE
CALL WRIT40('Shock correction requires: ')
CALL WRIT40(ERRMES)
CALL WRIT40('Please correct Q1 and run again. ')
CALL SET_ERR(214,'Error. See result file.',2)
ENDIF
C*****************************************************************
C
C--- GROUP 8. Terms (in differential equations) & devices
C
ELSEIF(IGR==8.AND.ISC==8) THEN
C
C * -----------GROUP 8 SECTION 8 --- CONVECTION COEFFICIENTS
C--- Entered when UCONV =.TRUE.
C
C--- Corrections applicable to velocities in compressible flow
IF(INDVAR==W1.AND.(NDIREC==5.OR.NDIREC==6)) THEN
C--- Correct W1 fluxes
CALL GXCNWH(NDIREC,I1H,L0FH,L0FL,L0CLL,HSOR)
ELSEIF(INDVAR==W2.AND.(NDIREC==5.OR.NDIREC==6)) THEN
C--- Correct W2 fluxes
CALL GXCNWH(NDIREC,I2H,L0FH,L0FL,L0CLL,HSOR)
ELSEIF(INDVAR==V1.AND.NDIREC==1) THEN
C--- Correct V1 fluxes
CALL GXCNVN(I1N,L0CH,L0CL,L0FN,L0FS,HSOR)
ELSEIF(INDVAR==V2.AND.NDIREC==1) THEN
C--- Correct V2 fluxes
CALL GXCNVN(I2N,L0CH,L0CL,L0FN,L0FL,HSOR)
ELSEIF(INDVAR==U1.AND.NDIREC==3) THEN
C--- Correct U1 fluxes
CALL GXCNUE(I1E,L0CH,L0CL,L0FE,L0FW,HSOR)
ELSEIF(INDVAR==U2.AND.NDIREC==3) THEN
C--- Correct U2 fluxes
CALL GXCNUE(I2E,L0CH,L0CL,L0FE,L0FW,HSOR)
ENDIF
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C Index for Coefficient - CO
C Index for Value - VAL
ELSEIF(IGR==13.AND.ISC==13) THEN
IF(INDVAR==H1.OR.INDVAR==ITEM) THEN
C---Enthalpy/Temperature pressure source term consistent with momentum
C modification.
IF(NPATCH(1:4)=='DPDT') THEN
C--- Add (P - Pold)/dt
CALL FN1(VAL,0.0)
IF(.NOT.STEADY) THEN
L0P1=L0F(P1); L0PO=L0F(OLD(P1)); L0VAL=L0F(VAL)
DO I = 1,NX*NY
IF(.NOT.SLD(I)) THEN
F(L0VAL+I)=HUNIT*(F(L0P1+I)-F(L0PO+I))/DT
ENDIF
ENDDO
ENDIF
ELSEIF(NPATCH(1:6)=='ENTSOR') THEN
CALL FN1(VAL,0.0)
C--- Add spatial contributions
IF(NX>1) CALL GXENTX(L0FE,L0FW,IU1,DEN1)
IF(NY>1) CALL GXENTY(L0FN,L0FS,IV1,DEN1)
IF(NZ>1) CALL GXENTZ(L0FH,L0FL,IW1,DEN1,L0AL)
ENDIF
ENDIF
ENDIF
END
C**********************************************************************
c
SUBROUTINE GXCNWH(NDIR,L0CON,L0FH,L0FL,L0CLL,HSOR)
C
C This subroutine operates on the High/Low fluxes for W1
C Input arguments:
C NDIR - Direction indicator; 5 = High, 6 = Low
C L0CON - Index for 3D scalar convection store (at IZ=1)
C HSOR - Logical switch for enthalpy pressure factors
C
C LD2 - Convection flux index for High and Low
C**********************************************************************
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
LOGICAL HSOR,SLD
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXCNWH'
C
C LOW CONVECTIONS
C ===============
C
IF(NDIR==6) THEN
IF(IZ>1) THEN
C--- Convection through low face of scalar cell
C--- Compute convection CP for low face of w-cell by averaging scalar convections
L0CL=L0CON+(IZ-2)*NX*NY; L0CH=L0CON+(IZ-1)*NX*NY
L0EASP20=L0F(EASP20)
DO I=1,NX*NY
IF(SLD(I).OR.SLD(I-NX*NY)) THEN ! current or low solid
F(L0EASP20+I)=0.0
ELSE ! average current and low
F(L0EASP20+I)=0.5*(F(L0CL+I)+F(L0CH+I))
ENDIF
ENDDO
ENDIF
C... Enthalpy pressure source factors
IF(HSOR) THEN
LLFL=L0FL+(IZ-1)*NX*NY
IF(IZ==1) THEN
CALL FN1(-LLFL,0.0)
ELSEIF(IZ==2) THEN
CALL FN1(-LLFL,1.0)
ELSE
CALL FN15(-LLFL,-L0CL,-L0CLL,0.0,1.0)
ENDIF
CALL FN0(-L0CLL,EASP20)
ENDIF
C... Finally remove outflows
CALL FN22(EASP20,0.0)
C
C HIGH CONVECTIONS
C ================
C
ELSEIF(NDIR==5) THEN
IF(IZ
SUBROUTINE GXCNVN(L0CON,L0CN,L0CS,L0FN,L0FS,HSOR)
C
C L0CON - Store for convection at NORTH face of scalar cell
C This subroutine operates on the North/South fluxes for V1
C Input arguments:
C L0CN - Store for convection at NORTH face of velocity cell
C L0CS - Store for convection at SOUTH face of velocity cell
C L0FS - Store for SOUTH face enthalpy pressure factor
C L0FN - Store for NORTH face enthalpy pressure factor
C HSOR - Logical switch for enthalpy pressure factors
C
C LD11 - Convection flux index for North
C LD12 - Convection flux index for South
C**********************************************************************
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
C
LOGICAL HSOR,SLD
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXCNVN'
C
C NORTH CONVECTIONS
C =================
C
C... Compute convection for v through north face by averaging
L0CONN=L0CON+(IZ-1)*NX*NY
DO 92 IX=1,NX
DO 90 IY=1,NY-2
I=(IX-1)*NY+IY
IF(SLD(I).OR.SLD(I+1)) THEN ! solid either side
F(L0CN+I)=0.0
ELSE
F(L0CN+I)=0.5*(F(L0CONN+I)+F(L0CONN+I+1))
ENDIF
90 CONTINUE
DO 91 IY=NY-1,NY
I=(IX-1)*NY+IY
F(L0CN+I)=0.0
91 CONTINUE
92 CONTINUE
C... Finally remove outflows
CALL FN66(LD11,-L0CN,0.0,-1.0)
C
C SOUTH CONVECTIONS
C =================
C
C--- Convection through south face of scalar cell
DO 101 IX=1,NX
II=1+NY*(IX-1)
F(L0CS+II)=0.0
DO 101 IY=2,NY
I=IY+NY*(IX-1)
F(L0CS+I)=F(L0CN+I-1)
101 CONTINUE
CALL FN66(LD12,-L0CS,0.0,1.0)
C... Enthalpy pressure term factors
IF(HSOR) THEN
LLFN=L0FN+(IZ-1)*NX*NY
LLFS=L0FS+(IZ-1)*NX*NY
DO 102 IX=1,NX
DO 102 IY=1,NY
I=(IX-1)*NY+IY
IF(IY<=2) THEN
F(LLFS+I)=1.0
ELSE
F(LLFS+I)=F(L0CONN+I-1)/(F(L0CN+I-2)+TINY)
ENDIF
IF(IY>=NY-1) THEN
F(LLFN+I)=1.0
ELSE
F(LLFN+I)=F(L0CONN+I)/(F(L0CN+I)+TINY)
ENDIF
102 CONTINUE
ENDIF
END
C**********************************************************************
c
SUBROUTINE GXCNUE(L0CON,L0CE,L0CW,L0FE,L0FW,HSOR)
C
C This subroutine operates on the East/West fluxes for U1
C Input arguments:
C L0CON - Store for convection at EAST face of scalar cell
C L0CE - Store for convection at EAST face of velocity cell
C L0CW - Store for convection at WEST face of velocity cell
C L0FE - Store for EAST face enthalpy pressure factor
C L0FW - Store for WEST face enthalpy pressure factor
C HSOR - Logical switch for enthalpy pressure factors
C
C LD11 - Convection flux index for East
C LD12 - Convection flux index for West
C**********************************************************************
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
C
LOGICAL HSOR,SLD
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXCNUE'
C
C EAST CONVECTIONS
C =================
C
C--- Convection through east face of scalar cell
L0CONE=L0CON+(IZ-1)*NX*NY
LIM=ITWO(NX-1,NX-2,XCYCLE)
C... Compute convection for U through east face by averaging
DO 90 I=1,NY*LIM ! common to XCYCLE T/F
IF(SLD(I).OR.SLD(I+NY)) THEN ! solid either side
F(L0CE+I)=0.0
ELSE
F(L0CE+I)=0.5*(F(L0CONE+I)+F(L0CONE+I+NY))
ENDIF
90 CONTINUE
IF(XCYCLE) THEN ! deal with XCYCLE cells at IX=NX
IADD=(NX-1)*NY
DO 91 IY=1,NY
IF(SLD(IY+IADD).OR.SLD(IY)) THEN ! solid either side
F(L0CE+IY+IADD)=0.0
ELSE
F(L0CE+IY+IADD)=0.5*(F(L0CONE+IY+IADD)+F(L0CONE+IY))
ENDIF
91 CONTINUE
ELSE ! for XCYCLE=F set edge convection to 0
DO 92 IX=NX-1,NX
DO 92 IY=1,NY
I=(IX-1)*NY+IY
F(L0CE+I)=0.0
92 CONTINUE
ENDIF
C... Finally remove outflows
CALL FN66(LD11,-L0CE,0.0,-1.0)
C
C WEST CONVECTIONS
C =================
C
C--- Convection through west face of velocity cell
IF(.NOT.XCYCLE) THEN
c... If not cyclic, set Cw(1) to 0.0
DO 100 IY=1,NY
F(L0CW+IY)=0.0
100 CONTINUE
ELSE
C... If cyclic, shift Ce(NX) to Cw(1)
IADD=(NX-1)*NY
DO 101 IY=1,NY
F(L0CW+IY)=F(L0CE+IY+IADD)
101 CONTINUE
ENDIF
C... Now shift Ce(IX-1) to Cw(IX)
DO 102 IX=2,NX
DO 102 IY=1,NY
I=IY+NY*(IX-1)
F(L0CW+I)=F(L0CE+I-NY)
102 CONTINUE
CALL FN66(LD12,-L0CW,0.0,1.0)
C... Enthalpy pressure source term factors
IF(HSOR) THEN
LLFE=L0FE+(IZ-1)*NX*NY
LLFW=L0FW+(IZ-1)*NX*NY
DO 104 IX=1,NX
DO 104 IY=1,NY
I=(IX-1)*NY+IY
IF(IX<=2) THEN
F(LLFW+I)=1.0
ELSE
F(LLFW+I)=F(L0CONE+I-NY)/(F(L0CE+I-NY*2)+TINY)
ENDIF
IF(IX>=NX-1) THEN
F(LLFE+I)=1.0
ELSE
F(LLFE+I)=F(L0CONE+I)/(F(L0CE+I)+TINY)
ENDIF
104 CONTINUE
ENDIF
END
C**********************************************************************
c
SUBROUTINE GXENTX(L0FE,L0FW,IUVEL,IRHO)
C
C This subroutine creates the East/West Enthalpy sources
C Input arguments:
C L0FE - Index for East face correction factor
C L0FW - Index for West face correction factor
C IUVEL - U velocity index
C IRHO - Phase density index
C**********************************************************************
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
logical sld, ef, wf
C
NAMSUB = 'GXENTX'
C
L0VAL=L0F(VAL); L0U1=L0F(IUVEL); L0D1=L0F(IRHO); L0AE=L0F(AEAST)
L0P1=L0F(P1); LLFE=L0FE+(IZ-1)*NX*NY; LLFW=L0FW+(IZ-1)*NX*NY
DO 131 I=1,NY*NX
IF(SLD(I)) GO TO 131
IX=(I-1)/NY+1
C... East face term
IF(IX1) THEN
IF(IX==2) THEN
IF(.NOT.XCYCLE) THEN
GW=F(L0D1+I-NY)*F(L0AE+I-NY)*F(L0U1+I-NY)
ELSE
IADD=(NX-1)*NY
IY=I-(IX-1)*NY
GW=F(LLFW+I)*F(L0D1+I-NY)*F(L0AE+I-NY)*0.5*(F(L0U1+I-NY)
1 +F(L0U1+IADD+IY))
ENDIF
ELSE
GW=F(LLFW+I)*F(L0D1+I-NY)*F(L0AE+I-NY)*0.5*(F(L0U1+I-NY)
1 +F(L0U1+I-2*NY))
ENDIF
ELSE
IF(.NOT.XCYCLE) THEN
GW=0.0
ELSE
IADD=(NX-1)*NY
GW=F(LLFW+I)*F(L0D1+IADD)*F(L0AE+IADD)*0.5*(F(L0U1+IADD)
1 +F(L0U1+IADD-NY))
ENDIF
ENDIF
C... Assemble E/W source
IF(.NOT.WF(I)) THEN
IF(IX>1) THEN
F(L0VAL+I)=F(L0VAL+I) +
1 AMAX1( GW,0.0)*(F(L0P1+I)-F(L0P1+I-NY))*HUNIT/F(L0D1+I-NY)
ELSEIF(XCYCLE) THEN
IADD=(NX-1)*NY
F(L0VAL+I)=F(L0VAL+I) +
1 AMAX1( GW,0.0)*(F(L0P1+I)-F(L0P1+I+IADD))*HUNIT/
2 F(L0D1+I+IADD)
ENDIF
ENDIF
IF(.NOT.EF(I)) THEN
IF(IX
SUBROUTINE GXENTY(L0FN,L0FS,IVVEL,IRHO)
C
C This subroutine creates the North/South Enthalpy sources
C Input arguments:
C L0FN - Index for North face correction factor
C L0FS - Index for South face correction factor
C IVVEL - V velocity index
C IRHO - Phase density index
C**********************************************************************
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
logical sld, nf, sf
C
NAMSUB = 'GXENTY'
C
L0VAL=L0F(VAL); L0V1=L0F(IVVEL); L0D1=L0F(IRHO); L0AN=L0F(ANORTH)
L0P1=L0F(P1); LLFN=L0FN+(IZ-1)*NX*NY; LLFS=L0FS+(IZ-1)*NX*NY
DO 131 IX=1,NX
DO 131 IY=1,NY
I=(IX-1)*NY+IY
IF(SLD(I)) GO TO 131
C... North face term
IF(IY1) THEN
IF(IY==2) THEN
GS=F(L0D1+I-1)*F(L0AN+I-1)*F(L0V1+I-1)
ELSE
GS=F(LLFS+I)*F(L0D1+I-1)*F(L0AN+I-1)*0.5*(F(L0V1+I-1)
1 +F(L0V1+I-2))
ENDIF
ELSE
GS=0.0
ENDIF
C... Assemble source
IF(.NOT.SF(I)) THEN
IF(IY>1) THEN
F(L0VAL+I)=F(L0VAL+I) +
1 AMAX1( GS,0.0)*(F(L0P1+I)-F(L0P1+I-1))*HUNIT/F(L0D1+I-1)
ENDIF
ENDIF
IF(.NOT.NF(I)) THEN
IF(IY
SUBROUTINE GXENTZ(L0FH,L0FL,IWVEL,IRHO,L0AL)
C
C This subroutine creates the High/Low Enthalpy sources
C Input arguments:
C L0FH - Index for High face correction factor
C L0FL - Index for Low face correction factor
C IWVEL - U velocity index
C IRHO - Phase density index
C**********************************************************************
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
COMMON/STOGEO/STORGEO
LOGICAL STORGEO
LOGICAL SLD, HF, LF
C
NAMSUB = 'GXENTZ'
C
L0VAL=L0F(VAL); L0W1=L0F(IWVEL); L0WH=L0F(HIGH(IWVEL))
L0WL=L0F(LOW(IWVEL)); L0WLL=L0F(ANYZ(IWVEL,IZ-2))
L0D1=L0F(IRHO); L0DH=L0F(HIGH(IRHO)); L0DL=L0F(LOW(IRHO))
L0AH=L0F(AHIGH); L0P1=L0F(P1)
L0PH=L0F(HIGH(P1)); L0PL=L0F(LOW(P1))
IF(STORGEO) THEN
L0AL=L0AH-NX*NY
ELSE
IF(BFC) THEN
CALL FNGBFC(-L0AL,3,IZ-1)
ELSE
IF(STORE(HPOR)) THEN
CALL FN21(-L0AL,LOW(HPOR),LXYAL,0.0,1.0)
ELSE
CALL FN0(-L0AL,AHIGH)
ENDIF
ENDIF
ENDIF
LLFH=L0FH+(IZ-1)*NX*NY; LLFL=L0FL+(IZ-1)*NX*NY
CALL FN1(-LLFL,1.0); CALL FN1(-LLFH,1.0)
DO 131 I=1,NX*NY
IF(SLD(I)) GO TO 131
C... High Face
IF(IZ1) THEN
IF(IZ==2) THEN
GL=F(L0DL+I)*F(L0AL+I)*F(L0WL+I)
ELSE
GL=F(LLFL+I)*F(L0DL+I)*F(L0AL+I)*0.5*(F(L0WL+I)
1 +F(L0WLL+I))
ENDIF
ELSE
GL=0.0
ENDIF
C... Assemble source
IF(.NOT.HF(I)) THEN
IF(IZ1) THEN
F(L0VAL+I)=F(L0VAL+I) +
1 AMAX1( GL,0.0)*(F(L0P1+I)-F(L0PL+I))*HUNIT/F(L0DL+I)
ENDIF
ENDIF
131 CONTINUE
END
C**********************************************************************
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXMACH(IMACH,IPTOT,IVABS,GAMA,IPH)
C.... This subroutine calculates the Mach Number, the total pressure
C and the absolute velocity. These are stored in the indices
C IMACH,IPTOT and IVABS respectively.
C
C The ratio of specific heats, gamma, is passed in via the
C argument GAMA.
C
C For the built in options RHO1/RHO2 = GRND3/GRND5, it is set
C automatically. Air is assumed if GAMA is zero.
C
C The argument IPH sets the phase to 1 or 2.
C
C If MACH is not STOREd, PTOT is calculated as P1+PRESS0+0.5*RHO*VABS^2
C
C It is called from Group 19.6 of GROUND.
C
C Author : J C Ludwig, CHAM
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'fgemcb'
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
LOGICAL QEQ,EQZ,SLD
SAVE IRUN0,L0VSQ,L0CSQ
SAVE GAMMA
DATA IRUN0 /0/
C
NAMSUB = 'GXMACH'
C
IF(IMACH/=0.OR.IPTOT/=0) THEN
IF(IRUN0/=IRUN) THEN
CALL GXMAKE(L0VSQ,NXYMAX,'VLSQ')
CALL GXMAKE(L0CSQ,NXYMAX,'SNSQ')
IRUN0=IRUN
ENDIF
GAMMA=GAMA
C... Compute VABS squared and store in L0VSQ
C... Get VABS if store present
IF(STORE(IVABS)) THEN
CALL FN71(-L0VSQ,IVABS,1.0,2.0)
ELSE
CALL FNVLSQ(-L0VSQ,IPH)
ENDIF
C... Select density store
IF(IPH==1) THEN
IRHO=DEN1
IF(QEQ(RHO1,GRND3)) GAMMA=1/RHO1B
IF(QEQ(RHO1,GRND5)) GAMMA=1/RHO1C
ELSE
IRHO=DEN2
IF(QEQ(RHO2,GRND3)) GAMMA=1/RHO2B
IF(QEQ(RHO2,GRND5)) GAMMA=1/RHO2C
ENDIF
C... If GAMMA not set, assume air
IF(EQZ(GAMMA)) GAMMA=1.4
IF(IMACH/=0) THEN ! MACH is stored, use compressible formula
C... Calculate square of sonic velocity from Csq = GAMMA.P1/RHO
CALL FN15(-L0CSQ,P1,IRHO,GAMMA*PRESS0,GAMMA)
C... Get Mach squared - Msq = VABSsq/Csq
CALL FN15(IMACH,-L0VSQ,-L0CSQ,0.0,1.0)
C... Compute total pressure if P0 is STOREd
IF(STORE(IPTOT)) THEN
L0P0=L0F(IPTOT); L0P=L0F(P1); L0M=L0F(IMACH)
ACON=GAMMA/(GAMMA-1.0); BCON=(GAMMA-1.0)/2.0
DO I=1,NX*NY
IF(SLD(I)) CYCLE
AMACH=AMAX1(0.0,AMIN1(100.0,F(L0M+I)))
F(L0P0+I)=(PRESS0+F(L0P+I))*ABS((1.0 + BCON*AMACH))**ACON
ENDDO
ENDIF
C... Now get square root of Mach Number
CALL FN30(IMACH)
ELSEIF(IPTOT/=0) THEN ! no MACH but PTOT present
L0P0=L0F(IPTOT); L0P=L0F(P1); L0RHO=L0F(IRHO)
DO I=1,NX*NY
IF(SLD(I)) CYCLE
F(L0P0+I)=PRESS0+F(L0P+I)+0.5*F(L0RHO+I)*F(L0VSQ+I)
ENDDO
ENDIF
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C SUBROUTINE GXCMPR is called from section 8, group 8 of GREX3, and
C is entered, for velocities only, when UCONV = T and COMPRS = T.
C
C.... The argument L0CNV is the F-array index for the convection
C coefficient; L0P and L0NEXP are P1 and P1(next), where 'next'
C is north, east or high; GA is a constant, set in GREX3 .
C
SUBROUTINE GXCMPR(LCNV,LP,LNEXTP,GA)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
COMMON /IDATA/NX,NY,NZ,LUPR1,IGFILL(116)
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXCMPR'
G1 = (GA-1.0)/GA
G2 = (1.0+1.0/GA)*0.3333333
G3 = 2.0*GA
L0CNV = L0F(LCNV)
L0P = L0F(LP)
L0NEXP = L0F(LNEXTP)
IXLAST=IXL
IF(INDVAR==3.OR.INDVAR==4) IXLAST=MIN0(IXL,NX-1)
IYLAST=IYL
IF(INDVAR==5.OR.INDVAR==6) IYLAST=MIN0(IYL,NY-1)
DO 10 IX = IXF,IXLAST
IADD = NY * (IX-1)
DO 10 IY = IYF,IYLAST
I = IY + IADD
INCONV = L0CNV + I
IF(F(INCONV)>0.0) THEN
PRAT = (F(L0NEXP+I) + PRESS0)/ABSPRES(I)
PRATM1 = 1.0 - PRAT
IF(ABS(PRATM1)<0.03) THEN
TOP = PRATM1*(1.0+G2*PRATM1)
F(INCONV) = F(INCONV)*(1.0-TOP/(G3+TOP))
ELSE
PRATG1=PRAT**G1
F(INCONV) = F(INCONV)*G1*PRATM1/(1.0-PRATG1)
ENDIF
ENDIF
10 CONTINUE
NAMSUB = 'gxcmpr'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c