PROGRAM RIEMANN C C This program computes the solution of a 1D c relativistic Riemann problem with C initial data UL if X<0.5 and UR if X>0.5 C in the whole spatial domain [0, 1] C CG GLOBAL DATA (COMMON BLOCKS) REFERENCED: CG /GAMMA/, /STATES/, /LS/, /RS/ C CF FILE ACCESS: CF solution.dat C CM MODULES CALLED: CM GEP, RAREF C PROGRAM RIEMAN IMPLICIT NONE C ------- C COMMON BLOCKS C ------- DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, WL, & RHOR, PR, UR, HR, CSR, VELR, WR COMMON /STATES/ RHOL, PL, UL, HL, CSL, VELL, WL, & RHOR, PR, UR, HR, CSR, VELR, WR DOUBLE PRECISION RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL COMMON /LS/ RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL DOUBLE PRECISION RHORS, URS, HRS, CSRS, VELRS, VSHOCKR COMMON /RS/ RHORS, URS, HRS, CSRS, VELRS, VSHOCKR DOUBLE PRECISION GAMMA COMMON /GAMMA/ GAMMA C --------- C INTERNAL VARIABLES C --------- INTEGER MN, N, I, ILOOP PARAMETER (MN = 400) DOUBLE PRECISION TOL, PMIN, PMAX, DVEL1, DVEL2, CHECK DOUBLE PRECISION PS, VELS DOUBLE PRECISION RHOA(MN), PA(MN), VELA(MN), UA(MN) DOUBLE PRECISION XI DOUBLE PRECISION RAD(MN), X1, X2, X3, X4, X5, T C ------- C INITIAL STATES C ------- WRITE(*,*) ' ADIABATIC INDEX OF THE GAS: ' READ (*,*) GAMMA WRITE(*,*) ' TIME FOR THE SOLUTION: ' READ (*,*) T C ----- C LEFT STATE C ----- WRITE(*,*) ' -- LEFT STATE -- ' WRITE(*,*) ' PRESSURE : ' READ (*,*) PL WRITE(*,*) ' DENSITY : ' READ (*,*) RHOL WRITE(*,*) ' FLOW VELOCITY: ' READ (*,*) VELL C ------ C RIGHT STATE C ------ WRITE(*,*) ' -- RIGHT STATE -- ' WRITE(*,*) ' PRESSURE : ' READ (*,*) PR WRITE(*,*) ' DENSITY : ' READ (*,*) RHOR WRITE(*,*) ' FLOW VELOCITY: ' READ (*,*) VELR C ------------------------------ C SPECIFIC INTERNAL ENERGY, SPECIFIC ENTHALPY, SOUND SPEED AND C FLOW LORENTZ FACTORS IN THE INITIAL STATES C ------------------------------ UL = PL/(GAMMA-1.D0)/RHOL UR = PR/(GAMMA-1.D0)/RHOR HL = 1.D0+UL+PL/RHOL HR = 1.D0+UR+PR/RHOR CSL = DSQRT(GAMMA*PL/RHOL/HL) CSR = DSQRT(GAMMA*PR/RHOR/HR) WL = 1.D0/DSQRT(1.D0-VELL**2) WR = 1.D0/DSQRT(1.D0-VELR**2) C -------- C NUMBER OF POINTS C -------- N = 400 C ------------- C TOLERANCE FOR THE SOLUTION C ------------- TOL = 0.D0 C ILOOP = 0 PMIN = (PL + PR)/2.D0 PMAX = PMIN 5 ILOOP = ILOOP + 1 PMIN = 0.5D0*MAX(PMIN,0.D0) PMAX = 2.D0*PMAX CALL GETDVEL(PMIN, DVEL1) CALL GETDVEL(PMAX, DVEL2) CHECK = DVEL1*DVEL2 IF (CHECK.GT.0.D0) GOTO 5 C --------------------------- C PRESSURE AND FLOW VELOCITY IN THE INTERMEDIATE STATES C --------------------------- CALL GETP(PMIN, PMAX, TOL, PS) VELS = 0.5D0*(VELLS + VELRS) C --------------- C SOLUTION ON THE NUMERICAL MESH C --------------- C ----------- C POSITIONS OF THE WAVES C ----------- IF (PL.GE.PS) THEN X1 = 0.5D0 + (VELL - CSL )/(1.D0 - VELL*CSL )*T X2 = 0.5D0 + (VELS - CSLS)/(1.D0 - VELS*CSLS)*T ELSE X1 = 0.5D0 + VSHOCKL*T X2 = X1 END IF X3 = 0.5D0 + VELS*T IF (PR.GE.PS) THEN X4 = 0.5D0 + (VELS + CSRS)/(1.D0 + VELS*CSRS)*T X5 = 0.5D0 + (VELR + CSR )/(1.D0 + VELR*CSR )*T ELSE X4 = 0.5D0 + VSHOCKR*T X5 = X4 END IF C ---------- C SOLUTION ON THE MESH C ---------- DO 100 I=1,N RAD(I) = DFLOAT(I)/DFLOAT(N) 100 CONTINUE DO 120 I=1,N IF (RAD(I).LE.X1) THEN PA(I) = PL RHOA(I) = RHOL VELA(I) = VELL UA(I) = UL ELSE IF (RAD(I).LE.X2) THEN XI = (RAD(I) - 0.5D0)/T CALL RAREF(XI, RHOL, PL, UL, CSL, VELL, 'L', & RHOA(I), PA(I), UA(I), VELA(I)) ELSE IF (RAD(I).LE.X3) THEN PA(I) = PS RHOA(I) = RHOLS VELA(I) = VELS UA(I) = ULS ELSE IF (RAD(I).LE.X4) THEN PA(I) = PS RHOA(I) = RHORS VELA(I) = VELS UA(I) = URS ELSE IF (RAD(I).LE.X5) THEN XI = (RAD(I) - 0.5D0)/T CALL RAREF(XI, RHOR, PR, UR, CSR, VELR, 'R', & RHOA(I), PA(I), UA(I), VELA(I)) ELSE PA(I) = PR RHOA(I) = RHOR VELA(I) = VELR UA(I) = UR END IF 120 CONTINUE OPEN (3,FILE='solution.dat',FORM='FORMATTED',STATUS='NEW') WRITE(3,150) N, T 150 FORMAT(I5,1X,F10.5) DO 60 I=1,N WRITE(3,200) RAD(I),PA(I),RHOA(I),VELA(I),UA(I) 60 CONTINUE 200 FORMAT(5(E15.8,1X)) CLOSE(3) STOP END C ---------- CN NAME: G E T D V E L C ---------- CP PURPOSE: CP COMPUTE THE DIFFERENCE IN FLOW SPEED BETWEEN LEFT AND RIGHT INTERMEDIATE CP STATES FOR GIVEN LEFT AND RIGHT STATES AND PRESSURE C CU USAGE: CU CALL GETDVEL( P, DVEL ) C CA ARGUMENTS: CA NAME / I,IO,O / TYPE(DIMENSION) CA DESCRIPTION CA P / I / DOUBLEPRECISION CA PRESSURE IN THE INTERMEDIATE STATE CA DVEL / O / DOUBLEPRECISION CA DIFFERENCE IN FLOW SPEED BETWEEN LEFT AND RIGHT INTERMEDIATE SATES C CG GLOBAL DATA (COMMON BLOCKS) REFERENCED: CG /GAMMA/, /STATES/ C CF FILE ACCESS: CF NONE C CM MODULES CALLED: CM GETVEL C SUBROUTINE GETDVEL( P, DVEL ) IMPLICIT NONE C ----- C ARGUMENTS C ----- DOUBLEPRECISION P, DVEL C ------- C COMMON BLOCKS C ------- DOUBLE PRECISION RHOLS,ULS,HLS,CSLS,VELLS,VSHOCKL COMMON /LS/ RHOLS,ULS,HLS,CSLS,VELLS,VSHOCKL DOUBLE PRECISION RHORS,URS,HRS,CSRS,VELRS,VSHOCKR COMMON /RS/ RHORS,URS,HRS,CSRS,VELRS,VSHOCKR DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, WL, & RHOR, PR, UR, HR, CSR, VELR, WR COMMON /STATES/ RHOL, PL, UL, HL, CSL, VELL, WL, & RHOR, PR, UR, HR, CSR, VELR, WR DOUBLE PRECISION GAMMA COMMON /GAMMA/ GAMMA C ----- C LEFT WAVE C ----- CALL GETVEL(P, RHOL, PL, UL, HL, CSL, VELL, WL, 'L', & RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL ) C ----- C RIGHT WAVE C ----- CALL GETVEL(P, RHOR, PR, UR, HR, CSR, VELR, WR, 'R', & RHORS, URS, HRS, CSRS, VELRS, VSHOCKR ) DVEL = VELLS - VELRS RETURN END C ------- CN NAME: G E T P C ------- CP PURPOSE: CP FIND THE PRESSURE IN THE INTERMEDIATE STATE OF A RIEMANN PROBLEM IN CP RELATIVISTIC HYDRODYNAMICS C CD DESCRIPTION: CD THIS ROUTINE USES A COMBINATION OF INTERVAL BISECTION AND INVERSE CD QUADRATIC INTERPOLATION TO FIND THE ROOT IN A SPECIFIED INTERVAL. C CU USAGE: CU CALL GETP( PMIN, PMAX, TOL, PS ) C CA ARGUMENTS: CA NAME / I,IO,O / TYPE(DIMENSION) CA DESCRIPTION CA PMIN / I / DOUBLEPRECISION CA THE LEFT ENDPOINT OF THE INTERVAL CA PMAX / I / DOUBLEPRECISION CA THE RIGHT ENDPOINT OF THE INTERVAL CA TOL / I / DOUBLEPRECISION CA TOLERANCE: DESIRED LENGTH OF THE INTERVAL OF UNCERTAINTY CA OF THE RESULT ( .GE. 0) CA PS / O / DOUBLEPRECISION CA PRESSURE IN THE INTERMEDIATE STATE WITHIN THE SPECIFIED TOLERANCE CA (OR TO MACHINE PRECISION IF THE TOLERANCE IS ZERO). C CG GLOBAL DATA (COMMON BLOCKS) REFERENCED: CG /GAMMA/, /STATES/ C CF FILE ACCESS: CF NONE C CM MODULES CALLED: CM GETDVEL C CC COMMENTS: CC IT IS ASSUMED THAT DVEL(PMIN) AND DVEL(PMAX) HAVE OPPOSITE SIGNS WITHOUT CC A CHECK. CC THIS ROUTINE IS FROM "COMPUTER METHODS FOR MATHEMATICAL COMPUTATION", CC BY G. E. FORSYTHE, M. A. MALCOLM, AND C. B. MOLER, CC PRENTICE-HALL, ENGLEWOOD CLIFFS N.J. C SUBROUTINE GETP( PMIN, PMAX, TOL, PS ) IMPLICIT NONE C ----- C ARGUMENTS C ----- DOUBLEPRECISION PMIN, PMAX, TOL, PS C ------- C COMMON BLOCKS C ------- DOUBLEPRECISION GAMMA COMMON /GAMMA/ GAMMA DOUBLEPRECISION RHOL, PL, UL, HL, CSL, VELL, WL, & RHOR, PR, UR, HR, CSR, VELR, WR COMMON /STATES/ RHOL, PL, UL, HL, CSL, VELL, WL, & RHOR, PR, UR, HR, CSR, VELR, WR C --------- C INTERNAL VARIABLES C --------- DOUBLEPRECISION A, B, C, D, E, EPS, FA, FB, FC, TOL1, & XM, P, Q, R, S C ------------- C COMPUTE MACHINE PRECISION C ------------- EPS = 1.D0 10 EPS = EPS/2.D0 TOL1 = 1.D0 + EPS IF( TOL1 .GT. 1.D0 ) GO TO 10 C ------- C INITIALIZATION C ------- A = PMIN B = PMAX CALL GETDVEL(A,FA) CALL GETDVEL(B,FB) C ----- C BEGIN STEP C ----- 20 C = A FC = FA D = B - A E = D 30 IF( DABS(FC) .GE. DABS(FB) )GO TO 40 A = B B = C C = A FA = FB FB = FC FC = FA C -------- C CONVERGENCE TEST C -------- 40 TOL1 = 2.D0*EPS*DABS(B) + 0.5D0*TOL XM = 0.5D0*(C - B) IF( DABS(XM) .LE. TOL1 ) GO TO 90 IF( FB .EQ. 0.D0 ) GO TO 90 C ------------ C IS BISECTION NECESSARY? C ------------ IF( DABS(E) .LT. TOL1 ) GO TO 70 IF( DABS(FA) .LE. DABS(FB) ) GO TO 70 C ------------------ C IS QUADRATIC INTERPOLATION POSSIBLE? C ------------------ IF( A .NE. C ) GO TO 50 C ---------- C LINEAR INTERPOLATION C ---------- S = FB/FA P = 2.D0*XM*S Q = 1.D0 - S GO TO 60 C ---------------- C INVERSE QUADRATIC INTERPOLATION C ---------------- 50 Q = FA/FC R = FB/FC S = FB/FA P = S*(2.D0*XM*Q*(Q - R) - (B - A)*(R - 1.D0)) Q = (Q - 1.D0)*(R - 1.D0)*(S - 1.D0) C ------ C ADJUST SIGNS C ------ 60 IF( P .GT. 0.D0 ) Q = -Q P = DABS(P) C -------------- C IS INTERPOLATION ACCEPTABLE? C -------------- IF( (2.D0*P) .GE. (3.D0*XM*Q-DABS(TOL1*Q)) ) GO TO 70 IF( P .GE. DABS(0.5D0*E*Q) ) GO TO 70 E = D D = P/Q GO TO 80 C ----- C BISECTION C ----- 70 D = XM E = D C ------- C COMPLETE STEP C ------- 80 A = B FA = FB IF( DABS(D) .GT. TOL1 ) B = B+D IF( DABS(D) .LE. TOL1 ) B = B+DSIGN(TOL1,XM) CALL GETDVEL(B,FB) IF( (FB*(FC/DABS(FC))) .GT. 0.D0) GO TO 20 GO TO 30 C -- C DONE C -- 90 PS = B RETURN END C --------- CN NAME: G E T V E L C --------- CP PURPOSE: CP COMPUTE THE FLOW VELOCITY BEHIND A RAREFACTION OR SHOCK IN TERMS OF THE CP POST-WAVE PRESSURE FOR A GIVEN STATE AHEAD THE WAVE IN A RELATIVISTIC CP FLOW C CD DESCRIPTION: C CU USAGE: CU CALL GETVEL( P, RHOA, PA, UA, HA, CSA, VELA, WA, S, VEL ) C CA ARGUMENTS: CA NAME / I,IO,O / TYPE(DIMENSION) CA DESCRIPTION CA P / I / DOUBLEPRECISION CA THE POST-WAVE PRESSURE CA RHOA / I / DOUBLEPRECISION CA THE DENSITY AHEAD THE WAVE CA PA / I / DOUBLEPRECISION CA THE PRESSURE AHEAD THE WAVE CA UA / I / DOUBLEPRECISION CA THE SPECIFIC INTERNAL ENERGY AHEAD THE WAVE CA HA / I / DOUBLEPRECISION CA THE SPECIFIC ENTHALPY AHEAD THE WAVE CA CSA / I / DOUBLEPRECISION CA THE LOCAL SOUND SPEED AHEAD THE WAVE CA VELA / I / DOUBLEPRECISION CA THE FLOW VELOCITY AHEAD THE WAVE CA WA / I / DOUBLEPRECISION CA THE FLOW LORENTZ FACTOR AHEAD THE WAVE CA S / I / CHARACTER CA THE DIRECTION OF PROPAGATION OF THE WAVE (LEFT OR RIGHT) CA RHO / O / DOUBLEPRECISION CA THE DENSITY IN THE POST-WAVE STATE CA U / O / DOUBLEPRECISION CA THE SPECIFIC INTERNAL ENERGY IN THE POST-WAVE STATE CA H / O / DOUBLEPRECISION CA THE SPECIFIC ENTHALPY IN THE POST-WAVE STATE CA CS / O / DOUBLEPRECISION CA THE LOCAL SOUND SPEED IN THE POST-WAVE STATE CA VEL / O / DOUBLEPRECISION CA THE FLOW VELOCITY IN THE POST-WAVE STATE CA VSHOCK / O / DOUBLEPRECISION CA THE SHOCK VELOCITY C CG GLOBAL DATA (COMMON BLOCKS) REFERENCED: CG /GAMMA/ C CF FILE ACCESS: CF NONE C CM MODULES CALLED: CM NONE C CC COMMENTS: CC THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN MARTI AND MUELLER, CC J. FLUID MECH., (1994) SUBROUTINE GETVEL( P, RHOA, PA, UA, HA, CSA, VELA, WA, S, & RHO, U, H, CS, VEL, VSHOCK ) IMPLICIT NONE C ----- C ARGUMENTS C ----- DOUBLE PRECISION P, RHOA, PA, UA, HA, CSA, VELA, WA CHARACTER*1 S DOUBLE PRECISION RHO, U, H, CS, VEL, VSHOCK C ------- C COMMON BLOCKS C ------- DOUBLE PRECISION GAMMA COMMON /GAMMA/ GAMMA C --------- C INTERNAL VARIABLES C --------- DOUBLE PRECISION A, B, C, SIGN DOUBLE PRECISION J, WSHOCK DOUBLE PRECISION K, SQGL1 C --------------- C LEFT OR RIGHT PROPAGATING WAVE C --------------- IF (S.EQ.'L') SIGN = -1.D0 IF (S.EQ.'R') SIGN = 1.D0 C IF (P.GT.PA) THEN C --- C SHOCK C --- A = 1.D0+(GAMMA-1.D0)*(PA-P)/GAMMA/P B = 1.D0-A C = HA*(PA-P)/RHOA-HA**2 C ---------------- C CHECK FOR UNPHYSICAL ENTHALPIES C ---------------- IF (C.GT.(B**2/4.D0/A)) STOP & 'GETVEL: UNPHYSICAL SPECIFIC ENTHALPY IN INTERMEDIATE STATE' C ----------------------------- C SPECIFIC ENTHALPY IN THE POST-WAVE STATE C (FROM THE EQUATION OF STATE AND THE TAUB ADIABAT, C EQ.(74), MM94) C ----------------------------- H = (-B+DSQRT(B**2-4.D0*A*C))/2.D0/A C --------------- C DENSITY IN THE POST-WAVE STATE C (FROM EQ.(73), MM94) C --------------- RHO = GAMMA*P/(GAMMA-1.D0)/(H-1.D0) C ------------------------ C SPECIFIC INTERNAL ENERGY IN THE POST-WAVE STATE C (FROM THE EQUATION OF STATE) C ------------------------ U = P/(GAMMA-1.D0)/RHO C -------------------------- C MASS FLUX ACROSS THE WAVE C (FROM THE RANKINE-HUGONIOT RELATIONS, EQ.(71), MM94) C -------------------------- J = SIGN*DSQRT((P-PA)/(HA/RHOA-H/RHO)) C ---------- C SHOCK VELOCITY C (FROM EQ.(86), MM94 C ---------- A = J**2+(RHOA*WA)**2 B = -VELA*RHOA**2*WA**2 VSHOCK = (-B+SIGN*J**2*DSQRT(1.D0+RHOA**2/J**2))/A WSHOCK = 1.D0/DSQRT(1.D0-VSHOCK**2) C ------------------- C FLOW VELOCITY IN THE POST-SHOCK STATE C (FROM EQ.(67), MM94) C ------------------- A = WSHOCK*(P-PA)/J+HA*WA*VELA B = HA*WA+(P-PA)*(WSHOCK*VELA/J+1.D0/RHOA/WA) VEL = A/B C --------------------- C LOCAL SOUND SPEED IN THE POST-SHOCK STATE C (FROM THE EQUATION OF STATE) C --------------------- CS = DSQRT(GAMMA*P/RHO/H) ELSE C ------ C RAREFACTION C ------ C --------------------------- C POLITROPIC CONSTANT OF THE GAS ACROSS THE RAREFACTION C --------------------------- K = PA/RHOA**GAMMA C --------------- C DENSITY BEHIND THE RAREFACTION C --------------- RHO = (P/K)**(1.D0/GAMMA) C ------------------------ C SPECIFIC INTERNAL ENERGY BEHIND THE RAREFACTION C (FROM THE EQUATION OF STATE) C ------------------------ U = P/(GAMMA-1.D0)/RHO C -------------------- C LOCAL SOUND SPEED BEHIND THE RAREFACTION C (FROM THE EQUATION OF STATE) C -------------------- CS = DSQRT(GAMMA*P/(RHO+GAMMA*P/(GAMMA-1.D0))) C ------------------ C FLOW VELOCITY BEHIND THE RAREFACTION C ------------------ SQGL1 = DSQRT(GAMMA-1.D0) A = (1.D0+VELA)/(1.D0-VELA)* & ((SQGL1+CSA)/(SQGL1-CSA)* & (SQGL1-CS )/(SQGL1+CS ))**(-SIGN*2.D0/SQGL1) VEL = (A-1.D0)/(A+1.D0) END IF END C -------- CN NAME: R A R E F C -------- CP PURPOSE: CP COMPUTE THE FLOW STATE IN A RAREFACTION FOR GIVEN PRE-WAVE STATE C CD DESCRIPTION: C CU USAGE: CU CALL RAREF( XI, RHOA, PA, UA, CSA, VELA, S, RHO, P, U, VEL ) C CA ARGUMENTS: CA NAME / I,IO,O / TYPE(DIMENSION) CA DESCRIPTION CA XI / I / DOUBLEPRECISION CA SIMILARITY VARIABLE CA RHOA / I / DOUBLEPRECISION CA THE DENSITY AHEAD THE RAREFACTION CA PA / I / DOUBLEPRECISION CA THE PRESSURE AHEAD THE RAREFACTION CA UA / I / DOUBLEPRECISION CA THE SPECIFIC INTERNAL ENERGY AHEAD THE RAREFACTION CA CSA / I / DOUBLEPRECISION CA THE LOCAL SOUND SPEED AHEAD THE RAREFACTION CA VELA / I / DOUBLEPRECISION CA THE FLOW VELOCITY AHEAD THE RAREFACTION CA S / I / CHARACTER CA THE DIRECTION OF PROPAGATION OF THE WAVE (LEFT OR RIGHT) CA RHO / O / DOUBLEPRECISION CA THE DENSITY IN A POINT WITHIN THE RAREFACTION CA P / O / DOUBLEPRECISION CA THE PRESSURE IN A POINT WITHIN THE RAREFACTION CA U / O / DOUBLEPRECISION CA THE SPECIFIC INTERNAL ENERGY IN A POINT WITHIN THE RAREFACTION CA CS / O / DOUBLEPRECISION CA THE LOCAL SOUND SPEED IN IN A POINT WITHIN THE RAREFACTION CA VEL / O / DOUBLEPRECISION CA THE FLOW VELOCITY IN IN A POINT WITHIN THE RAREFACTION C CG GLOBAL DATA (COMMON BLOCKS) REFERENCED: CG /GAMMA/ C CF FILE ACCESS: CF NONE C CM MODULES CALLED: CM NONE C CC COMMENTS: CC THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN MARTI AND MUELLER, CC J. FLUID MECH., (1994) SUBROUTINE RAREF( XI, RHOA, PA, UA, CSA, VELA, S, RHO, P, U, VEL ) IMPLICIT NONE C ----- C ARGUMENTS C ----- DOUBLE PRECISION XI DOUBLE PRECISION RHOA, PA, UA, CSA, VELA CHARACTER S DOUBLE PRECISION RHO, P, U, VEL C ------- C COMMON BLOCKS C ------- DOUBLE PRECISION GAMMA COMMON /GAMMA/ GAMMA C --------- C INTERNAL VARIABLES C --------- DOUBLE PRECISION B, C, D, K, L, V, OCS2, FCS2, DFDCS2, CS2, SIGN C --------------- C LEFT OR RIGHT PROPAGATING WAVE C --------------- IF (S.EQ.'L') SIGN = 1.D0 IF (S.EQ.'R') SIGN = -1.D0 B = DSQRT(GAMMA - 1.D0) C = (B + CSA)/(B - CSA) D = -SIGN*B/2.D0 K = (1.D0 + XI)/(1.D0 - XI) L = C*K**D V = ((1.D0 - VELA)/(1.D0 + VELA))**D OCS2 = CSA 25 FCS2 = L*V*(1.D0 + SIGN*OCS2)**D*(OCS2 - B) + & (1.D0 - SIGN*OCS2)**D*(OCS2 + B) DFDCS2 = L*V*(1.D0 + SIGN*OCS2)**D* & (1.D0 + SIGN*D*(OCS2 - B)/(1.D0 + SIGN*OCS2)) + & (1.D0 - SIGN*OCS2)**D* & (1.D0 - SIGN*D*(OCS2 + B)/(1.D0 - SIGN*OCS2)) CS2 = OCS2 - FCS2/DFDCS2 IF (ABS(CS2 - OCS2)/OCS2.GT.5.E-7)THEN OCS2 = CS2 GOTO 25 END IF VEL = (XI + SIGN*CS2)/(1.D0 + SIGN*XI*CS2) RHO = RHOA*((CS2**2*(GAMMA - 1.D0 - CSA**2))/ & (CSA**2*(GAMMA - 1.D0 - CS2**2))) & **(1.D0/(GAMMA - 1.D0)) P = CS2**2*(GAMMA - 1.D0)*RHO/(GAMMA - 1.D0 - CS2**2)/GAMMA U = P/(GAMMA - 1.D0)/RHO RETURN END .