PROGRAM IGMSS C----------------------------------------------------------------------- C A 3-D, STEADY-STATE, SUPERSONIC,SUPERALFVENIC SOLAR WIND MODEL IS C DEVELOPED BY NUMERICALLY INTEGRATING 8-SIMULTANEOUS MHD EQUATIONS. C A MODIFIED LAX-WENDROFF SCHEME IS USED. UNITS ARE IN KG,KM,SEC AND C KELVIN. C----------------------------------------------------------------------- C Input - unit 5 - file ss3din, note name is lower case in OPEN stmt C Output - unit 8 - file (or link) g3bssg - this is an BINARY file C The original version of igmss made an ASCII grid file C There is a separate program (named BING) to C convert those to binary format for use in IGM3D. C----------------------------------------------------------------------- C Changes: C 2002/09/25: Get Bfudge from 'Params.trec' (in ssmapbc.f) C 2001/05/23: Changed GAM from 1.67 to 1.46 [Totten, et al., JGR 1996] C 1989/02/22: CHANGE DTH, DPH INPUT UNITS TO DEGREES C 1990/11/05: INCREASE IX, INSTALL CHECK IM < OR = IX C IRUN=2, CHANGE TITLE -> 18-330 Rs C 1991/05/02: OUTPUT binary grid file g3bssg, mod input also C----------------------------------------------------------------------- PARAMETER ( PI = 3.141593 ) PARAMETER ( IX = 9480*3, JX = 30, KX = 72 ) COMMON/BLOC1/RO(JX,KX),VR(JX,KX),VT(JX,KX),VP(JX,KX), & BR(JX,KX),BT(JX,KX),BP(JX,KX),TP(JX,KX) COMMON/BLOC2/DTH,DPH,DR COMMON/BLOC3/TH(JX),PHI(KX),RAD(IX) COMMON/BLOC4/JM,KM,JP1,KP1,JM1,KM1,IM COMMON/BLOC5/ROP,VRP,VTP,VPP,BRP,BTP,BPP,TPP COMMON/BLOC6/FA(8),FB(8),FC(8),FD(8),FP(8) COMMON/BLOC7/SM,SG,RS,GAM,R COMMON/BLOC8/FE(8),FW(8),FN(8),FS(8) COMMON/BLOC9/GE(8),GW(8),GN(8),GS(8),GEN(8),GES(8),GWN(8),GWS(8),GP(8) COMMON/BLOC10/HE(8),HW(8),HN(8),HS(8),HEN(8),HES(8),HWN(8),HWS(8),HP(8) COMMON/BLOC11/SE(8),SW(8),SN(8),SS(8),SEN(8),SES(8),SWN(8),SWS(8),SP(8) COMMON/BLOC12/GA(8),GC(8),HB(8),HD(8),SA(8),SB(8),SC(8),SD(8) COMMON/BLOC13/F(8),G(8),H(8),S(8) DIMENSION RON(JX,KX),VRN(JX,KX),VTN(JX,KX),VPN(JX,KX),BRN(JX,KX), & BTN(JX,KX),BPN(JX,KX),TPN(JX,KX) DIMENSION PVR(120,31),PVT(120,31),PVP(120,31),PRO(120,31), & PTP(120,31),PBR(120,31),PBT(120,31),PBP(120,31) DATA GAM/1.46/,R/0.824E-2/,SM/1.988E30/,SG/6.668E-20/,RS/6.95E5/ CHARACTER*80 TITLE C --- INPUT FILE ss3bin SETS GRID GEOMETRY C --- Physical parameters at lower boundary are set in S.R. LBOUNDC C --- Description of inputs & operation instructions C TITLE C DPH,DTH,DRX - grid spacings ( 2*degrees, Rs ) (of output grid) C KMX,JMX,IMX - label dimensions of output grid C IRUN, ISEG - labels C KM,JM - internal (compute) dimensions of internal grid C IDPRT - internal subdivision of DRX: DR = DRX/IDPRT C --- The number of internal radial grid pts, IM = (IMX+1)*IDPRT C *** Special NOTE!!! IGM3D has special input routine for grid w/ T=0 C *** if T=0 it expects grid is one symmetric half plane j=1,16, k=1. C *** IGM3D uses S.R. FILL to duplicate half plane according to KM & JM C --- READ-IN INPUT DATA, AND ECHO IT OPEN(5,FILE='ss3bin') READ (5,'(A)') TITLE READ (5,*) DPX,DTX,DRX READ (5,*) KMX,JMX,IMX READ (5,*) IRUN, ISEG READ (5,*) KM,JM READ (5,*) IDPRT WRITE (6,'(A)') TITLE WRITE (6,*) DPX,DTX,DRX WRITE (6,*) KMX,JMX,IMX WRITE (6,*) IRUN, ISEG WRITE (6,*) KM,JM WRITE (6,*) IDPRT DTH = DTX*PI/180. DPH = DPX*PI/180. DR = DRX*RS/IDPRT IM = (IMX+1)*IDPRT IF( IM .GT. IX ) STOP 'EXPAND GRID ARRAY (IX)' JP1=JM+1 KP1=KM+1 JM1=JM-1 KM1=KM-1 CALL GRID WRITE (6,*) '--- radial output grid ---' WRITE (6,2000) (RAD(IK)/RS,IK=1,IM,IDPRT) I=1 IU=1 c --- Specify lower boundary condition CALL SSMAPBC ( ITIME) C --- Start new grid file OPEN (8,FILE='g3b.ss',FORM='UNFORMATTED') WRITE (8) TITLE(1:80) WRITE (8) ITIME, IRUN, ISEG, IMX, JMX, KMX, DRX, DTX, DPX DO 7 J=1,JM WRITE (8) ( RO(J,K),VR(J,K),VT(J,K),VP(J,K), & BR(J,K),BT(J,K),BP(J,K),TP(J,K), K=1,KM ) 7 CONTINUE WRITE (6,*) 'I=', IU, I WRITE (6,2002) ((J,K, RO(J,K),VR(J,K),VT(J,K),VP(J,K),BR(J,K) & ,BT(J,K),BP(J,K),TP(J,K),J=1,10),K=1,72,71) DO 8 J=1,JM WRITE (8) ( RO(J,K),VR(J,K),VT(J,K),VP(J,K), & BR(J,K),BT(J,K),BP(J,K),TP(J,K), K=1,KM ) 8 CONTINUE DO 11 J=1,JM PVR(IU,J)=VR(J,1) PVT(IU,J)=VT(J,1) PVP(IU,J)=VP(J,1) PRO(IU,J)=RO(J,1) PTP(IU,J)=TP(J,1) PBR(IU,J)=BR(J,1) PBT(IU,J)=BT(J,1) 11 PBP(IU,J)=BP(J,1) VrMinX = VR(1,1) VrMaxX = VR(1,1) TpMinX = Tp(1,1) DO 14 J=2,JM1 DO 14 K=2,KM1 VrMinX = MIN(VR(J,K),VrMinX) VrMaxX = MAX(VR(J,K),VrMaxX) TpMinX = MIN(Tp(J,K),TpMinX) 14 CONTINUE IPRINT=IDPRT +1 C --- Enter main loop - march outward in radius WRITE (6,*) '--- beginning radial march ---' 20 CONTINUE DO 100 J=2,JM1 DO 100 K=2,KM1 CALL FLUX1 (I,J,K) DO 200 IVAR=1,8 CALL EQN1 (I,J,K,IVAR) 200 CONTINUE CALL SIMULT (FA,1,I,J,K) CALL FLUX2 (1,I,J,K) CALL SIMULT (FC,2,I,J,K) CALL FLUX2 (2,I,J,K) CALL SIMULT (FB,3,I,J,K) CALL FLUX2 (3,I,J,K) CALL SIMULT (FD,4,I,J,K) CALL FLUX2 (4,I,J,K) DO 300 IVAR=1,8 CALL EQN2 (I,J,K,IVAR) 300 CONTINUE CALL SIMULT (FP,5,I,J,K) RON(J,K)=ROP VRN(J,K)=VRP VTN(J,K)=VTP VPN(J,K)=VPP BRN(J,K)=BRP BTN(J,K)=BTP BPN(J,K)=BPP TPN(J,K)=TPP 100 CONTINUE VrMin = VR(1,1) VrMax = VR(1,1) TpMin = Tp(1,1) DO 400 J=2,JM1 DO 400 K=2,KM1 RO(J,K) = RON(J,K) VR(J,K) = VRN(J,K) VrMin = MIN(VR(J,K),VrMin) VrMax = MAX(VR(J,K),VrMax) TpMin = MIN(Tp(J,K),TpMin) VT(J,K) = VTN(J,K) VP(J,K) = VPN(J,K) BR(J,K) = BRN(J,K) BT(J,K) = BTN(J,K) BP(J,K) = BPN(J,K) TP(J,K) = TPN(J,K) 400 CONTINUE DVrMin = VrMin - VrMinX DVrMax = VrMax - VrMaxX VrMinX = VrMin VrMaxX = VrMax TpMinX = TpMin CALL ARVIS(I) CALL BOUND IF (I.LT.IPRINT) GO TO 450 IU=IU+1 DO 401 J=1,JM WRITE (8) ( RO(J,K),VR(J,K),VT(J,K),VP(J,K), & BR(J,K),BT(J,K),BP(J,K),TP(J,K), K=1,KM ) 401 CONTINUE IF (VrMin .LT. 0.0) STOP 'NEGATIVE Vr' IF (TpMin .LT. 0.0) STOP 'NEGATIVE Tp' DO 430 J=1,JM PVR(IU,J)=VR(J,1) PVT(IU,J)=VT(J,1) PVP(IU,J)=VP(J,1) PRO(IU,J)=RO(J,1) PTP(IU,J)=TP(J,1) PBR(IU,J)=BR(J,1) PBT(IU,J)=BT(J,1) PBP(IU,J)=BP(J,1) 430 CONTINUE IPRINT = IPRINT+IDPRT 450 CONTINUE WRITE (6,2003) 'Iu=', IU,MOD(I,IDPRT), VrMin,DVrMin,VrMax, TpMin I = I+1 C WRITE (6,2002) ((J,K, RO(J,K),VR(J,K),VT(J,K),VP(J,K),BR(J,K) C & ,BT(J,K),BP(J,K),TP(J,K),J=1,10),K=1,72,71) IF (I.GE.IM) GO TO 500 GO TO 20 1000 FORMAT(1H ) 2000 FORMAT(8E11.3) 2002 FORMAT(2I3,8E11.3) 2003 FORMAT(A,2I3,F8.1,F8.3,F8.1,F12.1) 500 CONTINUE STOP END C --- ------------------------------------------------------------------ SUBROUTINE GRID PARAMETER ( IX = 9480*3, JX = 30, KX = 72 ) PARAMETER ( PI = 3.141593 ) COMMON/BLOC2/DTH,DPH,DR COMMON/BLOC3/TH(JX),PHI(KX),RAD(IX) COMMON/BLOC4/JM,KM,JP1,KP1,JM1,KM1,IM COMMON/BLOC7/SM,SG,RS,GAM,R C TH(1) = 17.5*PI/180. PHI(1) = 2.5*PI/180. RAD(1) = 21.5*RS C DO 10 J=2,JM TH(J)=TH(J-1)+DTH 10 CONTINUE DO 20 K=2,KM PHI(K)=PHI(K-1)+DPH 20 CONTINUE DO 30 II=2,IM RAD(II)=RAD(II-1)+DR 30 CONTINUE RETURN END C --- ------------------------------------------------------------------ SUBROUTINE FLUX1 (I,J,K) COMMON/BLOC6/FA(8),FB(8),FC(8),FD(8),FP(8) COMMON/BLOC8/FE(8),FW(8),FN(8),FS(8) COMMON/BLOC9/GE(8),GW(8),GN(8),GS(8),GEN(8),GES(8),GWN(8),GWS(8),GP(8) COMMON/BLOC10/HE(8),HW(8),HN(8),HS(8),HEN(8),HES(8),HWN(8),HWS(8),HP(8) COMMON/BLOC11/SE(8),SW(8),SN(8),SS(8),SEN(8),SES(8),SWN(8),SWS(8),SP(8) COMMON/BLOC13/F(8),G(8),H(8),S(8) C CALL FLARG(I,J,K,0) DO 10 II=1,8 FP(II)=F(II) GP(II)=G(II) HP(II)=H(II) SP(II)=S(II) 10 CONTINUE C CALL FLARG (I,J+1,K,0) DO 20 II=1,8 FE(II)=F(II) GE(II)=G(II) HE(II)=H(II) SE(II)=S(II) 20 CONTINUE C CALL FLARG (I,J-1,K,0) DO 30 II=1,8 FW(II)=F(II) GW(II)=G(II) HW(II)=H(II) SW(II)=S(II) 30 CONTINUE C CALL FLARG (I,J,K+1,0) DO 40 II=1,8 FN(II)=F(II) GN(II)=G(II) HN(II)=H(II) SN(II)=S(II) 40 CONTINUE C CALL FLARG (I,J,K-1,0) DO 50 II=1,8 FS(II)=F(II) GS(II)=G(II) HS(II)=H(II) SS(II)=S(II) 50 CONTINUE C CALL FLARG (I,J+1,K+1,0) DO 60 II=1,8 GEN(II)=G(II) HEN(II)=H(II) SEN(II)=S(II) 60 CONTINUE C CALL FLARG (I,J+1,K-1,0) DO 70 II=1,8 GES(II)=G(II) HES(II)=H(II) SES(II)=S(II) 70 CONTINUE C CALL FLARG (I,J-1,K+1,0) DO 80 II=1,8 GWN(II)=G(II) HWN(II)=H(II) SWN(II)=S(II) 80 CONTINUE C CALL FLARG (I,J-1,K-1,0) DO 90 II=1,8 GWS(II)=G(II) HWS(II)=H(II) SWS(II)=S(II) 90 CONTINUE C RETURN END C --- ------------------------------------------------------------------ SUBROUTINE FLARG (I,M,N,JFLAG) PARAMETER ( IX = 9480*3, JX = 30, KX = 72 ) COMMON/BLOC1/RO(JX,KX),VR(JX,KX),VT(JX,KX),VP(JX,KX), & BR(JX,KX),BT(JX,KX),BP(JX,KX),TP(JX,KX) COMMON/BLOC2/DTH,DPH,DR COMMON/BLOC3/TH(JX),PHI(KX),RAD(IX) COMMON/BLOC5/ROP,VRP,VTP,VPP,BRP,BTP,BPP,TPP COMMON/BLOC7/SM,SG,RS,GAM,R COMMON/BLOC13/F(8),G(8),H(8),S(8) C IF (JFLAG.GT.0) GO TO 10 ROP=RO(M,N) VRP=VR(M,N) VTP=VT(M,N) VPP=VP(M,N) BRP=BR(M,N) BTP=BT(M,N) BPP=BP(M,N) TPP=TP(M,N) 10 CONTINUE T1=RAD(I) THM=TH(M) IF (JFLAG.GT.0) T1=RAD(I+1) IF (JFLAG.EQ.1) THM=TH(M)+0.25*DTH IF (JFLAG.EQ.2) THM=TH(M)-0.5*DTH T0=T1**2 XMU1=4.*3.1416 XMU2=2.*XMU1 F(1)=T0*ROP*VRP TEMP1=2.*ROP*R*TPP+ROP*VRP**2 TEMP2=(-BRP**2+BTP**2+BPP**2)/XMU2 F(2)=T0*(TEMP1+TEMP2) F(3)=T0*(ROP*VRP*VTP-BRP*BTP/XMU1) F(4)=T0*(ROP*VRP*VPP-BRP*BPP/XMU1) F(5)=T0*BRP F(6)=T1*(VRP*BTP-VTP*BRP) F(7)=T1*(VRP*BPP-VPP*BRP) T2=(GAM/(GAM-1.))*2.*R TEMP1=VRP*ROP*(T2*TPP+0.5*(VRP**2+VTP**2+VPP**2)) TEMP2=BTP*(VRP*BTP-VTP*BRP)/XMU1 TEMP3=BPP*(VPP*BRP-VRP*BPP)/XMU1 F(8)=T0*(TEMP1+TEMP2-TEMP3) C G(1)=T0*SIN(THM)*ROP*VTP G(2)=T0*SIN(THM)*(ROP*VRP*VTP-BRP*BTP/XMU1) G(3)=T0*SIN(THM)*(2.*ROP*R*TPP+ROP*VTP**2+(BRP**2-BTP**2+BPP**2)/XMU2) G(4)=T0*SIN(THM)*(ROP*VTP*VPP-BTP*BPP/XMU1) G(5)=T0*SIN(THM)*BTP G(6)=0. G(7)=T1*SIN(THM)*(VTP*BPP-VPP*BTP) TEMP1=VTP*ROP*(T2*TPP+0.5*(VRP**2+VTP**2+VPP**2)) TEMP2=BRP*(VRP*BTP-VTP*BRP)/XMU1 TEMP3=BPP*(VTP*BPP-VPP*BTP)/XMU1 G(8)=T0*SIN(THM)*(TEMP1-TEMP2+TEMP3) C H(1)=T0*ROP*VPP H(2)=T0*(ROP*VRP*VPP-BRP*BPP/XMU1) H(3)=T0*(ROP*VTP*VPP-BTP*BPP/XMU1) H(4)=T0*(2.*ROP*R*TPP+ROP*VPP**2+(BRP**2+BTP**2-BPP**2)/XMU2) H(5)=T0*BPP H(6)=T1*(VPP*BTP-VTP*BPP) H(7)=0. TEMP1=VPP*ROP*(T2*TPP+.5*(VRP**2+VTP**2+VPP**2)) TEMP2=BRP*(VPP*BRP-VRP*BPP)/XMU1 TEMP3=BTP*(VTP*BPP-VPP*BTP)/XMU1 H(8)=T0*(TEMP1+TEMP2-TEMP3) C S(1)=0. TEMP1=2.*T1*2.*ROP*R*TPP TEMP2=-ROP*SM*SG+T1*BRP**2/XMU1 TEMP3=ROP*T1*(VTP**2+VPP**2) S(2)=TEMP1+TEMP2+TEMP3 TEMP1=T1*1./TAN(THM)*2.*ROP*R*TPP TEMP2=(BRP**2+BTP**2-BPP**2)/XMU2*T1*1./TAN(THM) TEMP3=T1*BRP*BTP/XMU1-ROP*T1*(VRP*VTP-VPP**2*1./TAN(THM)) S(3)=TEMP1+TEMP2+TEMP3 TEMP1=BTP*BPP*T1*1./TAN(THM)/XMU1 TEMP2=T1*BRP*BPP/XMU1 TEMP3=-ROP*T1*(VRP*VPP+VTP*VPP*1./TAN(THM)) S(4)=TEMP1+TEMP2+TEMP3 S(5)=0. S(6)=0. S(7)=(VTP*BPP-VPP*BTP)*1./TAN(THM) S(8)=-ROP*VRP*SM*SG C RETURN END C --- ------------------------------------------------------------------ SUBROUTINE EQN1 (I,J,K,L) PARAMETER ( IX = 9480*3, JX = 30, KX = 72 ) COMMON/BLOC2/DTH,DPH,DR COMMON/BLOC3/TH(JX),PHI(KX),RAD(IX) COMMON/BLOC6/FA(8),FB(8),FC(8),FD(8),FP(8) COMMON/BLOC8/FE(8),FW(8),FN(8),FS(8) COMMON/BLOC9/GE(8),GW(8),GN(8),GS(8),GEN(8),GES(8),GWN(8),GWS(8),GP(8) COMMON/BLOC10/HE(8),HW(8),HN(8),HS(8),HEN(8),HES(8),HWN(8),HWS(8),HP(8) COMMON/BLOC11/SE(8),SW(8),SN(8),SS(8),SEN(8),SES(8),SWN(8),SWS(8),SP(8) C T3=RAD(I) C T3=RAD(I)+0.5*DR T2=(TH(J-1)+TH(J))/2. T1=(TH(J)+TH(J+1))/2. FA(L)=0.5*(FP(L)+FE(L))-(GE(L)-GP(L))*DR/(T3*SIN(T1)*DTH) 1 -((HEN(L)-HES(L))/SIN(TH(J+1))+(HN(L)-HS(L))/SIN(TH(J)))*DR/ 2 (4.*T3*DPH)+0.5*DR*(SE(L)+SP(L)) FC(L)=0.5*(FP(L)+FW(L))-(GP(L)-GW(L))*DR/(T3*SIN(T2)*DTH) 1 -((HN(L)-HS(L))/SIN(TH(J))+(HWN(L)-HWS(L))/SIN(TH(J-1)))*DR/ 2 (4.*T3*DPH)+0.5*DR*(SP(L)+SW(L)) FB(L)=0.5*(FP(L)+FN(L))-(GEN(L)-GWN(L)+GE(L)-GW(L))*DR/ 1 (4.*T3*SIN(TH(J))*DTH)-(HN(L)-HP(L))*DR/(T3*SIN(TH(J))* 2 DPH)+0.5*DR*(SN(L)+SP(L)) FD(L)=0.5*(FP(L)+FS(L))-(GE(L)-GW(L)+GES(L)-GWS(L))*DR/ 1 (4.*T3*SIN(TH(J))*DTH)-(HP(L)-HS(L))*DR/(T3*SIN(TH(J))* 2 DPH)+0.5*DR*(SP(L)+SS(L)) C RETURN END C --- ------------------------------------------------------------------ SUBROUTINE SIMULT (FL,IFLAG,I,J,K) PARAMETER ( IX = 9480*3, JX = 30, KX = 72 ) COMMON/BLOC1/RO(JX,KX),VR(JX,KX),VT(JX,KX),VP(JX,KX), & BR(JX,KX),BT(JX,KX),BP(JX,KX),TP(JX,KX) COMMON/BLOC2/DTH,DPH,DR COMMON/BLOC3/TH(JX),PHI(KX),RAD(IX) COMMON/BLOC5/ROP,VRP,VTP,VPP,BRP,BTP,BPP,TPP COMMON/BLOC7/SM,SG,RS,GAM,R DIMENSION FL(8) DATA ERRORV/0.00001/ C TEMP0=RAD(I+1) IF (IFLAG.EQ.1) GO TO 10 IF (IFLAG.EQ.2) GO TO 20 IF (IFLAG.EQ.3) GO TO 30 IF (IFLAG.EQ.4) GO TO 40 IF (IFLAG.EQ.5) GO TO 50 C 10 CONTINUE VR0=0.5*(VR(J+1,K)+VR(J,K)) BR0=FL(5)/TEMP0**2 GO TO 100 20 CONTINUE VR0=0.5*(VR(J,K)+VR(J-1,K)) BR0=FL(5)/TEMP0**2 GO TO 100 30 CONTINUE VR0=0.5*(VR(J,K+1)+VR(J,K)) BR0=FL(5)/TEMP0**2 GO TO 100 40 CONTINUE VR0=0.5*(VR(J,K)+VR(J,K-1)) BR0=FL(5)/TEMP0**2 GO TO 100 50 CONTINUE VR0=VR(J,K) BR0=FL(5)/TEMP0**2 C 100 CONTINUE KITER=1 150 CONTINUE TEMP=FL(6)+TEMP0*FL(3)/FL(1)*BR0 TEMP2=TEMP0*VR0-TEMP0**3*BR0**2/(4.*3.14*FL(1)) BTP=TEMP/TEMP2 TEMP1=FL(7)+TEMP0*FL(4)/FL(1)*BR0 BPP=TEMP1/TEMP2 A=-0.5*(GAM+1.)/(GAM-1.)*FL(1) B=GAM*TEMP0**2*BR0**2/((GAM-1.)*4.*3.14*2.)+(GAM-2.)*TEMP0 1 **2*(BTP**2+BPP**2)/((GAM-1.)*4.*3.14*2.)+GAM*FL(2)/(GAM-1.) C=(FL(3)**2+FL(4)**2)*0.5/FL(1)-FL(8) 1 -TEMP0**4*BR0**2*(BTP**2+BPP**2)/(2.*(4.*3.14)**2*FL(1)) TEMPA=B**2-4.*A*C IF (TEMPA.LT.0.) GO TO 140 TEMP=SQRT(TEMPA) VR1=(-B+TEMP)*0.5/A VR2=(-B-TEMP)*0.5/A IF (VR1.GE.VR0.AND.VR2.GE.VR0) GO TO 190 VRP=AMAX1(VR1,VR2) GO TO 195 190 VRP=AMAX1(VR1,VR2) 195 CONTINUE ERROR=ABS((VR0-VRP)/VRP) IF (ERROR.LT.ERRORV) GO TO 200 IF (KITER.GT.40) GO TO 200 KITER=KITER+1 VR0=VRP GO TO 150 200 CONTINUE GO TO 250 140 VRP=VR0 250 CONTINUE ROP=FL(1)/(TEMP0**2*VRP) VTP=FL(3)/FL(1)+TEMP0**2*BR0*BTP/(4.*3.14*FL(1)) VPP=FL(4)/FL(1)+TEMP0**2*BR0*BPP/(4.*3.14*FL(1)) BRP=BR0 TEMP1=FL(8)/(TEMP0**2*VRP)-0.5*ROP*(VRP**2+VTP**2+VPP**2) TEMP2=-(BTP*FL(6)+BPP*FL(7))/(4.*3.14*TEMP0*VRP) P=(GAM-1.)/GAM*(TEMP1+TEMP2) TPP=P/(2.*ROP*R) 2000 FORMAT (1X,8E9.3) RETURN END C --- ------------------------------------------------------------------ SUBROUTINE FLUX2 (IFLAG,I,J,K) COMMON/BLOC5/ROP,VRP,VTP,VPP,BRP,BTP,BPP,TPP COMMON/BLOC7/SM,SG,RS,GAM,R COMMON/BLOC12/GA(8),GC(8),HB(8),HD(8),SA(8),SB(8),SC(8),SD(8) COMMON/BLOC13/F(8),G(8),H(8),S(8) C IF (IFLAG.EQ.1) GO TO 10 IF (IFLAG.EQ.2) GO TO 20 IF (IFLAG.EQ.3) GO TO 30 IF (IFLAG.EQ.4) GO TO 40 10 CONTINUE CALL FLARG (I,J,K,1) DO 50 II=1,8 GA(II)=G(II) SA(II)=S(II) 50 CONTINUE GO TO 500 20 CONTINUE CALL FLARG (I,J,K,2) DO 60 II=1,8 GC(II)=G(II) SC(II)=S(II) 60 CONTINUE GO TO 500 30 CONTINUE CALL FLARG (I,J,K,3) DO 70 II=1,8 HB(II)=H(II) SB(II)=S(II) 70 CONTINUE GO TO 500 40 CONTINUE CALL FLARG (I,J,K,4) DO 80 II=1,8 HD(II)=H(II) SD(II)=S(II) 80 CONTINUE 500 CONTINUE RETURN END C --- ------------------------------------------------------------------ SUBROUTINE EQN2 (I,J,K,L) PARAMETER ( IX = 9480*3, JX = 30, KX = 72 ) COMMON/BLOC2/DTH,DPH,DR COMMON/BLOC3/TH(JX),PHI(KX),RAD(IX) COMMON/BLOC6/FA(8),FB(8),FC(8),FD(8),FP(8) COMMON/BLOC9/GE(8),GW(8),GN(8),GS(8),GEN(8),GES(8),GWN(8),GWS(8),GP(8) COMMON/BLOC10/HE(8),HW(8),HN(8),HS(8),HEN(8),HES(8),HWN(8),HWS(8),HP(8) COMMON/BLOC11/SE(8),SW(8),SN(8),SS(8),SEN(8),SES(8),SWN(8),SWS(8),SP(8) COMMON/BLOC12/GA(8),GC(8),HB(8),HD(8),SA(8),SB(8),SC(8),SD(8) C T2=RAD(I) C T2=RAD(I)+0.5*DR T3=RAD(I+1) C T3=T2 FP(L)=FP(L)-.5*DR*((GE(L)-GW(L))/(2.*T2*SIN(TH(J))*DTH) 1 +(GA(L)-GC(L))/(T3*SIN(TH(J))*DTH)) 2 -0.5*DR*((HN(L)-HS(L))/(2.*T2*SIN(TH(J))*DPH) 3 +(HB(L)-HD(L))/(T3*SIN(TH(J))*DPH)) 4 +0.5*DR*(SP(L)+0.25*(SA(L)+SB(L)+SC(L)+SD(L))) C RETURN END C --- ------------------------------------------------------------------ SUBROUTINE BOUND C TREAT COMPUTATIONAL BOUNDARIES EXCEPT LOWER. PARAMETER ( IX = 9480*3, JX = 30, KX = 72 ) COMMON/BLOC1/RO(JX,KX),VR(JX,KX),VT(JX,KX),VP(JX,KX), & BR(JX,KX),BT(JX,KX),BP(JX,KX),TP(JX,KX) COMMON/BLOC4/JM,KM,JP1,KP1,JM1,KM1,IM C K = 1 & K = KM BOUNDARYS: DO 10 J=1,JM RO(J, 1)=RO(J, 2) RO(J,KM)=RO(J,KM1) VR(J, 1)=VR(J, 2) VR(J,KM)=VR(J,KM1) VT(J, 1)=VT(J, 2) VT(J,KM)=VT(J,KM1) VP(J, 1)=VP(J, 2) VP(J,KM)=VP(J,KM1) BR(J, 1)=BR(J, 2) BR(J,KM)=BR(J,KM1) BT(J, 1)=BT(J, 2) BT(J,KM)=BT(J,KM1) BP(J, 1)=BP(J, 2) BP(J,KM)=BP(J,KM1) TP(J, 1)=TP(J, 2) TP(J,KM)=TP(J,KM1) 10 CONTINUE C J = 1 & J = JM BOUNDARYS: JM2=JM1-1 DO 20 K=1,KM RO(1,K) = 1.1*RO(2,K) - 0.1*RO(3,K) VR(1,K) = 1.1*VR(2,K) - 0.1*VR(3,K) VT(1,K) = 1.1*VT(2,K) - 0.1*VT(3,K) VP(1,K) = 1.1*VP(2,K) - 0.1*VP(3,K) BR(1,K) = 1.1*BR(2,K) - 0.1*BR(3,K) BT(1,K) = 1.1*BT(2,K) - 0.1*BT(3,K) BP(1,K) = 1.1*BP(2,K) - 0.1*BP(3,K) TP(1,K) = 1.1*TP(2,K) - 0.1*TP(3,K) RO(JM,K) = 1.1*RO(JM1,K) - 0.1*RO(JM2,K) VR(JM,K) = 1.1*VR(JM1,K) - 0.1*VR(JM2,K) VT(JM,K) = 1.1*VT(JM1,K) - 0.1*VT(JM2,K) VP(JM,K) = 1.1*VP(JM1,K) - 0.1*VP(JM2,K) BR(JM,K) = 1.1*BR(JM1,K) - 0.1*BR(JM2,K) BT(JM,K) = 1.1*BT(JM1,K) - 0.1*BT(JM2,K) BP(JM,K) = 1.1*BP(JM1,K) - 0.1*BP(JM2,K) TP(JM,K) = 1.1*TP(JM1,K) - 0.1*TP(JM2,K) 20 CONTINUE RETURN END