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