c *****************************************************************
c SPLITTING IN THE Y-DIRECTION
c *****************************************************************
SUBROUTINE CHARY
IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N)
COMMON/CELL/ ROCELL(2001,2001),UCELL(2001,2001),VCELL(2001,2001),
& ETCELL(2001,2001),EICELL(2001,2001)
COMMON/GXN/UXN(2001,2001),VXN(2001,2001),ROXN(2001,2001)
& ,PXN(2001,2001)
COMMON/GX/UX(2001,2001),VX(2001,2001),ROX(2001,2001),PX(2001,2001)
COMMON/GY/UY(2001,2001),VY(2001,2001),ROY(2001,2001),PY(2001,2001)
COMMON/XY/X(2001),Y(2001)
COMMON/D1X/PN(2001),UN(2001),RON(2001),VN(2001)
COMMON/BAL/NX,NY,CFL,EPS,GAM,DT
COMMON/FON/ROFONL,EIFONL,PFONL,UFONL,VFONL,ROFONR,EIFONR
& ,PFONR,UFONR,VFONR
RELAX=0.000
c *****************************************************************
DO 1 I=1,NX-1
HX=X(I+1)-X(I)
J=1
UR=UX(I+1,J)
UL=UX(I,J)
VR=VX(I+1,J)
VL=VX(I,J)
PR=PX(I+1,J)
PL=PX(I,J)
ROR=ROX(I+1,J)
ROL=ROX(I,J)
EIR=PR/(ROR*(GAM-1))
EIL=PL/(ROL*(GAM-1))
c ------------------------------------------------
UT=UY(I,J+1)
UB=UY(I,J)
VT=VY(I,J+1)
VB=VY(I,J)
PT=PY(I,J+1)
PB=PY(I,J)
ROT=ROY(I,J+1)
ROB=ROY(I,J)
EIT=PT/(ROT*(GAM-1))
EIB=PB/(ROB*(GAM-1))
c ------------------------------------------------
ROC=ROCELL(I,J)
UC=UCELL(I,J)
VC=VCELL(I,J)
ETC=ETCELL(I,J)
EIC=EICELL(I,J)
PC=(GAM-1)*ROC*EIC
SOUND=DSQRT(GAM*PC/ROC)
c -------------------------------------------------
c PRIMARY VAR, AX
AX11=UC
AX12=ROC
AX13=0
AX14=0
AX21=0
AX22=AX11
AX23=0
AX24=1/AX12
AX31=0
AX32=0
AX33=AX11
AX34=0
AX41=0
AX42=GAM*PC
AX43=0
AX44=AX11
DFIDX1=(ROR-ROL)/HX
DFIDX2=(UR-UL)/HX
DFIDX3=(VR-VL)/HX
DFIDX4=(PR-PL)/HX
FY1=-(AX11*DFIDX1+AX12*DFIDX2+AX13*DFIDX3+AX14*DFIDX4)
FY2=-(AX21*DFIDX1+AX22*DFIDX2+AX23*DFIDX3+AX24*DFIDX4)
FY3=-(AX31*DFIDX1+AX32*DFIDX2+AX33*DFIDX3+AX34*DFIDX4)
FY4=-(AX41*DFIDX1+AX42*DFIDX2+AX43*DFIDX3+AX44*DFIDX4)
c CHAR VAR, OMEGAX
SCW=PC/ROC**GAM
G=2*DSQRT(GAM)/(GAM-1)*SCW**(1/(2*GAM))
c G=2./DSQRT(GAM*ROC)
OY11=0
OY12=0
OY13=1
OY14=2./DSQRT(GAM*ROC)/(2*DSQRT(PC))
OY21=0
OY22=1
OY23=0
OY24=0
OY31=0
OY32=0
OY33=1
OY34=-OY14
OY41=-GAM/ROC
OY42=0
OY43=0
OY44=1/PC
GY1=(OY11*FY1+OY12*FY2+OY13*FY3+OY14*FY4)
GY2=(OY21*FY1+OY22*FY2+OY23*FY3+OY24*FY4)
GY3=(OY31*FY1+OY32*FY2+OY33*FY3+OY34*FY4)
GY4=(OY41*FY1+OY42*FY2+OY43*FY3+OY44*FY4)
RT=VT+G*(PT)**((GAM-1)/(2*GAM))
RB=VB+G*(PB)**((GAM-1)/(2*GAM))
RC=VC+G*(PC)**((GAM-1)/(2*GAM))
QT=VT-G*(PT)**((GAM-1)/(2*GAM))
QB=VB-G*(PB)**((GAM-1)/(2*GAM))
QC=VC-G*(PC)**((GAM-1)/(2*GAM))
ST=DLOG(PT/ROT**GAM)
SB=DLOG(PB/ROB**GAM)
SC=DLOG(PC/ROC**GAM)
RMAX=max(RT,RB,RC)+DT*GY1
RMIN=min(RT,RB,RC)+DT*GY1
UMAX=max(UT,UB,UC)+DT*GY2
UMIN=min(UT,UB,UC)+DT*GY2
QMAX=max(QT,QB,QC)+DT*GY3
QMIN=min(QT,QB,QC)+DT*GY3
SMAX=max(ST,SB,SC)+DT*GY4
SMIN=min(ST,SB,SC)+DT*GY4
RTN=2*RC-RB
RBN=2*RC-RT
QTN=2*QC-QB
QBN=2*QC-QT
STN=2*SC-SB
SBN=2*SC-ST
UTN=2*UC-UB
UBN=2*UC-UT
IF(RTN.GT.RMAX) RTN=RMAX
IF(RTN.LT.RMIN) RTN=RMIN
IF(QTN.GT.QMAX) QTN=QMAX
IF(QTN.LT.QMIN) QTN=QMIN
IF(UTN.GT.UMAX) UTN=UMAX
IF(UTN.LT.UMIN) UTN=UMIN
IF(STN.GT.SMAX) STN=SMAX
IF(STN.LT.SMIN) STN=SMIN
c ------------------------------------------
IF(RBN.GT.RMAX) RBN=RMAX
IF(RBN.LT.RMIN) RBN=RMIN
IF(QBN.GT.QMAX) QBN=QMAX
IF(QBN.LT.QMIN) QBN=QMIN
IF(UBN.GT.UMAX) UBN=UMAX
IF(UBN.LT.UMIN) UBN=UMIN
IF(SBN.GT.SMAX) SBN=SMAX
IF(SBN.LT.SMIN) SBN=SMIN
RTNB=RTN
QTNB=QTN
UTNB=UTN
STNB=STN
GB=G
SOUNDB=SOUND
VCB=VC
DO 2 J=2,NY-1
c -------------------------------------------
UR=UX(I+1,J)
UL=UX(I,J)
VR=VX(I+1,J)
VL=VX(I,J)
PR=PX(I+1,J)
PL=PX(I,J)
ROR=ROX(I+1,J)
ROL=ROX(I,J)
EIR=PR/(ROR*(GAM-1))
EIL=PL/(ROL*(GAM-1))
c ------------------------------------------------
UT=UY(I,J+1)
UB=UY(I,J)
VT=VY(I,J+1)
VB=VY(I,J)
PT=PY(I,J+1)
PB=PY(I,J)
ROT=ROY(I,J+1)
ROB=ROY(I,J)
EIT=PT/(ROT*(GAM-1))
EIB=PB/(ROB*(GAM-1))
c ------------------------------------------------
ROC=ROCELL(I,J)
UC=UCELL(I,J)
VC=VCELL(I,J)
ETC=ETCELL(I,J)
EIC=EICELL(I,J)
PC=(GAM-1)*ROC*EIC
SOUND=DSQRT(GAM*PC/ROC)
c -------------------------------------------------
HY=Y(J+1)-Y(J)
c -------------------------------------------------
c PRIMARY VAR, AX
AX11=UC
AX12=ROC
AX13=0
AX14=0
AX21=0
AX22=AX11
AX23=0
AX24=1/AX12
AX31=0
AX32=0
AX33=AX11
AX34=0
AX41=0
AX42=GAM*PC
AX43=0
AX44=AX11
DFIDX1=(ROR-ROL)/HX
DFIDX2=(UR-UL)/HX
DFIDX3=(VR-VL)/HX
DFIDX4=(PR-PL)/HX
FY1=-(AX11*DFIDX1+AX12*DFIDX2+AX13*DFIDX3+AX14*DFIDX4)
FY2=-(AX21*DFIDX1+AX22*DFIDX2+AX23*DFIDX3+AX24*DFIDX4)
FY3=-(AX31*DFIDX1+AX32*DFIDX2+AX33*DFIDX3+AX34*DFIDX4)
FY4=-(AX41*DFIDX1+AX42*DFIDX2+AX43*DFIDX3+AX44*DFIDX4)
c CHAR VAR, OMEGAX
SCW=PC/ROC**GAM
G=2*DSQRT(GAM)/(GAM-1)*SCW**(1/(2*GAM))
c G=2./DSQRT(GAM*ROC)
OY11=0
OY12=0
OY13=1
OY14=2./DSQRT(GAM*ROC)/(2*DSQRT(PC))
OY21=0
OY22=1
OY23=0
OY24=0
OY31=0
OY32=0
OY33=1
OY34=-OY14
OY41=-GAM/ROC
OY42=0
OY43=0
OY44=1/PC
GY1=(OY11*FY1+OY12*FY2+OY13*FY3+OY14*FY4)
GY2=(OY21*FY1+OY22*FY2+OY23*FY3+OY24*FY4)
GY3=(OY31*FY1+OY32*FY2+OY33*FY3+OY34*FY4)
GY4=(OY41*FY1+OY42*FY2+OY43*FY3+OY44*FY4)
RT=VT+G*(PT)**((GAM-1)/(2*GAM))
RB=VB+G*(PB)**((GAM-1)/(2*GAM))
RC=VC+G*(PC)**((GAM-1)/(2*GAM))
QT=VT-G*(PT)**((GAM-1)/(2*GAM))
QB=VB-G*(PB)**((GAM-1)/(2*GAM))
QC=VC-G*(PC)**((GAM-1)/(2*GAM))
ST=DLOG(PT/ROT**GAM)
SB=DLOG(PB/ROB**GAM)
SC=DLOG(PC/ROC**GAM)
RMAX=(1+RELAX)*max(RT,RB,RC)+DT*GY1
RMIN=(1-RELAX)*min(RT,RB,RC)+DT*GY1
UMAX=max(UT,UB,UC)+DT*GY2
UMIN=min(UT,UB,UC)+DT*GY2
QMAX=(1-RELAX)*max(QT,QB,QC)+DT*GY3
QMIN=(1+RELAX)*min(QT,QB,QC)+DT*GY3
SMAX=max(ST,SB,SC)+DT*GY4
SMIN=min(ST,SB,SC)+DT*GY4
RTN=2*RC-RB
RBN=2*RC-RT
QTN=2*QC-QB
QBN=2*QC-QT
STN=2*SC-SB
SBN=2*SC-ST
UTN=2*UC-UB
UBN=2*UC-UT
IF(RTN.GT.RMAX) RTN=RMAX
IF(RTN.LT.RMIN) RTN=RMIN
IF(QTN.GT.QMAX) QTN=QMAX
IF(QTN.LT.QMIN) QTN=QMIN
IF(UTN.GT.UMAX) UTN=UMAX
IF(UTN.LT.UMIN) UTN=UMIN
IF(STN.GT.SMAX) STN=SMAX
IF(STN.LT.SMIN) STN=SMIN
c ------------------------------------------
IF(RBN.GT.RMAX) RBN=RMAX
IF(RBN.LT.RMIN) RBN=RMIN
IF(QBN.GT.QMAX) QBN=QMAX
IF(QBN.LT.QMIN) QBN=QMIN
IF(UBN.GT.UMAX) UBN=UMAX
IF(UBN.LT.UMIN) UBN=UMIN
IF(SBN.GT.SMAX) SBN=SMAX
IF(SBN.LT.SMIN) SBN=SMIN
c CHAR DIRECTIONS
VCH=VC+VCB !!-VB
SCH=SOUND+SOUNDB !!-DSQRT(GAM*PB/ROB)
CHAR3=VCH
CHAR4=VCH
CHAR1=VCH+SCH
CHAR2=VCH-SCH
c UPDATE NODAL VARS
IF(CHAR1.GT.0.AND.CHAR2.LE.0.AND.CHAR3.GT.0) THEN
PN(J)=((RTNB-QBN)/(G+GB))**(2*GAM/(GAM-1))
VN(J)=(G*RTNB+GB*QBN)/(G+GB)
RON(J)=(PN(J)/DEXP(STNB))**(1/GAM)
UN(J)=UTNB
GOTO 20
ENDIF
IF(CHAR1.GT.0.AND.CHAR2.LE.0.AND.CHAR3.LE.0) THEN
PN(J)=((RTNB-QBN)/(G+GB))**(2*GAM/(GAM-1))
VN(J)=(G*RTNB+GB*QBN)/(G+GB)
RON(J)=(PN(J)/DEXP(SBN))**(1/GAM)
UN(J)=UBN
GOTO 20
ENDIF
IF(CHAR1.GT.0.AND.CHAR2.GT.0.AND.CHAR3.GT.0) THEN
PN(J)=((RTNB-QTNB)/(GB+GB))**(2*GAM/(GAM-1))
VN(J)=(GB*RTNB+GB*QTNB)/(GB+GB)
RON(J)=(PN(J)/DEXP(STNB))**(1/GAM)
UN(J)=UTNB
GOTO 20
ENDIF
IF(CHAR1.LE.0.AND.CHAR2.LE.0.AND.CHAR3.LE.0) THEN
PN(J)=((RBN-QBN)/(G+G))**(2*GAM/(GAM-1))
VN(J)=(G*RBN+G*QBN)/(G+G)
RON(J)=(PN(J)/DEXP(SBN))**(1/GAM)
UN(J)=UBN
GOTO 20
ENDIF
20 CONTINUE
RTNB=RTN
QTNB=QTN
UTNB=UTN
STNB=STN
GB=G
SOUNDB=SOUND
VCB=VC
2 CONTINUE
DO 3 J=1,NY-1
PY(I,J)=PN(J)
VY(I,J)=VN(J)
ROY(I,J)=RON(J)
IF(ROY(I,J).LT.EPS) ROY(I,J)=EPS
UY(I,J)=UN(J)
PN(J)=0
UN(J)=0
RON(J)=0
VN(J)=0
3 CONTINUE
c -----------------------------------------------------------------
1 CONTINUE
DO 1000 I=1,NX
DO 1000 J=1,NY
UX(I,J)=UXN(I,J)
VX(I,J)=VXN(I,J)
ROX(I,J)=ROXN(I,J)
PX(I,J)=PXN(I,J)
UXN(I,J)=0
VXN(I,J)=0
ROXN(I,J)=0
PXN(I,J)=0
1000 CONTINUE
RETURN
END