Warning: This page is an HTML page, i.e. contains markup. Follow the download-link to obtain an unaltered copy of the file the content of which is displayed below.

1
C     ------------
2
CN    NAME: R I E M A N N 
3
C     ------------
4

        
5
CP    PURPOSE:
6
CP    THIS PROGRAM COMPUTES THE SOLUTION OF A 1D   
7
CP    RELATIVISTIC RIEMANN PROBLEM (FOR CONSTANT-GAMMA IDEAL GASES) WITH  
8
CP    INITIAL DATA UL IF X<0.5 AND UR IF X>0.5  
9
CP    IN THE WHOLE SPATIAL DOMAIN [0, 1] 
10
C
11

        
12
CC    COMMENTS:
13
CC    SEE MARTI AND MUELLER, JFM, 1994
14
CC
15
CC    WRITTEN BY:     Jose-Maria Marti
16
CC                    Departamento de Astronomia y Astrofisica 
17
CC                    Universidad de Valencia 
18
CC                    46100 Burjassot (Valencia), Spain
19
CC                    jose-maria.marti@uv.es
20
CC    AND
21
CC                    Ewald Mueller
22
CC                    Max-Planck-Institut fuer Astrophysik
23
CC                    Karl-Schwarzschild-Str. 1
24
CC                    85741 Garching, Germany
25
CC                    emueller@mpa-garching.mpg.de
26
C
27

        
28
      PROGRAM RIEMANN
29

        
30
      IMPLICIT NONE
31

        
32
C     ------- 
33
C     COMMON BLOCKS 
34
C     -------
35

        
36
      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, WL, 
37
     &                 RHOR, PR, UR, HR, CSR, VELR, WR 
38
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, WL, 
39
     &                 RHOR, PR, UR, HR, CSR, VELR, WR
40

        
41
      DOUBLE PRECISION RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL 
42
      COMMON /LS/      RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL
43

        
44
      DOUBLE PRECISION RHORS, URS, HRS, CSRS, VELRS, VSHOCKR 
45
      COMMON /RS/      RHORS, URS, HRS, CSRS, VELRS, VSHOCKR
46

        
47
      DOUBLE PRECISION GAMMA 
48
      COMMON /ADIND/   GAMMA
49

        
50
C     --------- 
51
C     INTERNAL VARIABLES 
52
C     ---------
53

        
54
      INTEGER         MN, N, I, ILOOP 
55
      PARAMETER      (MN = 400)
56

        
57
      DOUBLE PRECISION TOL, PMIN, PMAX, DVEL1, DVEL2, CHECK
58

        
59
      DOUBLE PRECISION PS, VELS
60

        
61
      DOUBLE PRECISION RHOA(MN), PA(MN), VELA(MN), UA(MN)
62

        
63
      DOUBLE PRECISION XI
64

        
65
      DOUBLE PRECISION RAD(MN), X1, X2, X3, X4, X5, T
66

        
67
C     ------- 
68
C     INITIAL STATES 
69
C     -------
70

        
71
      WRITE(*,*) ' ADIABATIC INDEX OF THE GAS: ' 
72
      READ (*,*)   GAMMA
73

        
74
      WRITE(*,*) ' TIME FOR THE SOLUTION: ' 
75
      READ (*,*)   T
76

        
77
C     -----
78
C     LEFT STATE
79
C     -----
80

        
81
      WRITE(*,*) ' -- LEFT STATE -- ' 
82
      WRITE(*,*) '      PRESSURE     : ' 
83
      READ (*,*) PL 
84
      WRITE(*,*) '      DENSITY      : ' 
85
      READ (*,*) RHOL 
86
      WRITE(*,*) '      FLOW VELOCITY: ' 
87
      READ (*,*) VELL
88

        
89
C     ------ 
90
C     RIGHT STATE 
91
C     ------
92

        
93
      WRITE(*,*) ' -- RIGHT STATE -- ' 
94
      WRITE(*,*) '      PRESSURE     : ' 
95
      READ (*,*) PR 
96
      WRITE(*,*) '      DENSITY      : ' 
97
      READ (*,*) RHOR 
98
      WRITE(*,*) '      FLOW VELOCITY: ' 
99
      READ (*,*) VELR
100

        
101
C     ------------------------------ 
102
C     SPECIFIC INTERNAL ENERGY, SPECIFIC ENTHALPY, SOUND SPEED AND  
103
C     FLOW LORENTZ FACTORS IN THE INITIAL STATES 
104
C     ------------------------------ 
105

        
106
      UL  = PL/(GAMMA-1.D0)/RHOL 
107
      UR  = PR/(GAMMA-1.D0)/RHOR
108

        
109
      HL  = 1.D0+UL+PL/RHOL 
110
      HR  = 1.D0+UR+PR/RHOR
111

        
112
      CSL = DSQRT(GAMMA*PL/RHOL/HL) 
113
      CSR = DSQRT(GAMMA*PR/RHOR/HR)
114

        
115
      WL  = 1.D0/DSQRT(1.D0-VELL**2) 
116
      WR  = 1.D0/DSQRT(1.D0-VELR**2)
117

        
118
C     -------- 
119
C     NUMBER OF POINTS 
120
C     --------
121

        
122
      N   = 400
123

        
124
C     ------------- 
125
C     TOLERANCE FOR THE SOLUTION 
126
C     -------------
127

        
128
      TOL = 0.D0
129

        
130
C
131

        
132
      ILOOP = 0
133

        
134
      PMIN  = (PL + PR)/2.D0 
135
      PMAX  = PMIN
136

        
137
 5    ILOOP = ILOOP + 1
138

        
139
      PMIN  = 0.5D0*MAX(PMIN,0.D0) 
140
      PMAX  = 2.D0*PMAX
141

        
142
      CALL GETDVEL(PMIN, DVEL1)
143

        
144
      CALL GETDVEL(PMAX, DVEL2)
145

        
146
      CHECK = DVEL1*DVEL2 
147
      IF (CHECK.GT.0.D0) GOTO 5
148

        
149
C     --------------------------- 
150
C     PRESSURE AND FLOW VELOCITY IN THE INTERMEDIATE STATES 
151
C     ---------------------------
152

        
153
      CALL GETP(PMIN, PMAX, TOL, PS)
154

        
155
      VELS = 0.5D0*(VELLS + VELRS)
156

        
157
C     --------------- 
158
C     SOLUTION ON THE NUMERICAL MESH 
159
C     ---------------
160

        
161
C     ----------- 
162
C     POSITIONS OF THE WAVES 
163
C     -----------
164

        
165
      IF (PL.GE.PS) THEN
166

        
167
        X1 = 0.5D0 + (VELL - CSL )/(1.D0 - VELL*CSL )*T 
168
        X2 = 0.5D0 + (VELS - CSLS)/(1.D0 - VELS*CSLS)*T
169

        
170
      ELSE
171

        
172
        X1 = 0.5D0 + VSHOCKL*T 
173
        X2 = X1
174

        
175
      END IF
176

        
177
      X3 = 0.5D0 + VELS*T
178

        
179
      IF (PR.GE.PS) THEN
180

        
181
        X4 = 0.5D0 + (VELS + CSRS)/(1.D0 + VELS*CSRS)*T 
182
        X5 = 0.5D0 + (VELR + CSR )/(1.D0 + VELR*CSR )*T
183

        
184
      ELSE
185

        
186
        X4 = 0.5D0 + VSHOCKR*T 
187
        X5 = X4
188

        
189
      END IF
190

        
191
C     ---------- 
192
C     SOLUTION ON THE MESH 
193
C     ----------
194

        
195
      DO 100 I=1,N
196

        
197
        RAD(I) = DFLOAT(I)/DFLOAT(N)
198

        
199
 100  CONTINUE
200

        
201
      DO 120 I=1,N
202

        
203
        IF (RAD(I).LE.X1) THEN
204

        
205
          PA(I)   = PL 
206
          RHOA(I) = RHOL 
207
          VELA(I) = VELL 
208
          UA(I)   = UL
209

        
210
        ELSE IF (RAD(I).LE.X2) THEN
211

        
212
          XI = (RAD(I) - 0.5D0)/T
213

        
214
          CALL RAREF(XI, RHOL,    PL,    UL,    CSL,  VELL,  'L', 
215
     &                   RHOA(I), PA(I), UA(I),       VELA(I))
216

        
217
        ELSE IF (RAD(I).LE.X3) THEN
218

        
219
          PA(I)   = PS 
220
          RHOA(I) = RHOLS 
221
          VELA(I) = VELS 
222
          UA(I)   = ULS
223

        
224
        ELSE IF (RAD(I).LE.X4) THEN
225

        
226
          PA(I)   = PS 
227
          RHOA(I) = RHORS 
228
          VELA(I) = VELS 
229
          UA(I)   = URS
230

        
231
        ELSE IF (RAD(I).LE.X5) THEN
232

        
233
          XI = (RAD(I) - 0.5D0)/T
234

        
235
          CALL RAREF(XI, RHOR,    PR,    UR,    CSR,  VELR,  'R', 
236
     &                   RHOA(I), PA(I), UA(I),       VELA(I))
237

        
238
        ELSE
239

        
240
          PA(I)   = PR 
241
          RHOA(I) = RHOR 
242
          VELA(I) = VELR 
243
          UA(I)   = UR
244

        
245
        END IF
246

        
247
 120  CONTINUE
248

        
249
      OPEN (3,FILE='solution.dat',FORM='FORMATTED',STATUS='NEW')
250

        
251
      WRITE(3,150) N, T 
252
 150  FORMAT(I5,1X,F10.5)
253

        
254
      DO 60 I=1,N 
255
        WRITE(3,200) RAD(I),PA(I),RHOA(I),VELA(I),UA(I) 
256
 60   CONTINUE
257

        
258
 200  FORMAT(5(E15.8,1X))
259

        
260
      CLOSE(3)
261

        
262
      STOP 
263
      END
264

        
265
C     ---------- 
266
CN    NAME: G E T D V E L 
267
C     ----------
268

        
269
CP    PURPOSE: 
270
CP    COMPUTE THE DIFFERENCE IN FLOW SPEED BETWEEN LEFT AND RIGHT INTERMEDIATE 
271
CP    STATES FOR GIVEN LEFT AND RIGHT STATES AND PRESSURE 
272
C 
273

        
274
CC    COMMENTS
275
CC    NONE
276

        
277
      SUBROUTINE GETDVEL( P, DVEL )
278

        
279
      IMPLICIT NONE
280

        
281
C     ----- 
282
C     ARGUMENTS 
283
C     -----
284

        
285
      DOUBLEPRECISION P, DVEL
286

        
287
C     ------- 
288
C     COMMON BLOCKS 
289
C     -------
290

        
291
      DOUBLE PRECISION RHOLS,ULS,HLS,CSLS,VELLS,VSHOCKL 
292
      COMMON /LS/      RHOLS,ULS,HLS,CSLS,VELLS,VSHOCKL
293

        
294
      DOUBLE PRECISION RHORS,URS,HRS,CSRS,VELRS,VSHOCKR 
295
      COMMON /RS/      RHORS,URS,HRS,CSRS,VELRS,VSHOCKR
296

        
297
      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, WL, 
298
     &                 RHOR, PR, UR, HR, CSR, VELR, WR 
299
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, WL, 
300
     &                 RHOR, PR, UR, HR, CSR, VELR, WR
301

        
302
      DOUBLE PRECISION GAMMA 
303
      COMMON /ADIND/   GAMMA
304

        
305
C     ----- 
306
C     LEFT WAVE 
307
C     -----
308

        
309
      CALL GETVEL(P, RHOL, PL, UL,  HL,  CSL,  VELL,  WL, 'L', 
310
     &               RHOLS,    ULS, HLS, CSLS, VELLS, VSHOCKL )
311

        
312
C     ----- 
313
C     RIGHT WAVE 
314
C     -----
315

        
316
      CALL GETVEL(P, RHOR, PR, UR,  HR,  CSR,  VELR,  WR, 'R', 
317
     &               RHORS,    URS, HRS, CSRS, VELRS, VSHOCKR )
318

        
319
      DVEL = VELLS - VELRS
320

        
321
      RETURN 
322
      END
323

        
324
C     ------- 
325
CN    NAME: G E T P 
326
C     -------
327

        
328
CP    PURPOSE: 
329
CP    FIND THE PRESSURE IN THE INTERMEDIATE STATE OF A RIEMANN PROBLEM IN 
330
CP    RELATIVISTIC HYDRODYNAMICS 
331
C 
332

        
333
CC    COMMENTS: 
334
CC    THIS ROUTINE USES A COMBINATION OF INTERVAL BISECTION AND INVERSE 
335
CC    QUADRATIC INTERPOLATION TO FIND THE ROOT IN A SPECIFIED INTERVAL. 
336
CC    IT IS ASSUMED THAT DVEL(PMIN) AND DVEL(PMAX) HAVE OPPOSITE SIGNS WITHOUT 
337
CC    A CHECK. 
338
CC    ADAPTED FROM "COMPUTER METHODS FOR MATHEMATICAL COMPUTATION", 
339
CC    BY G. E. FORSYTHE, M. A. MALCOLM, AND C. B. MOLER, 
340
CC    PRENTICE-HALL, ENGLEWOOD CLIFFS N.J. 
341
C
342
      SUBROUTINE GETP( PMIN, PMAX, TOL, PS )
343

        
344
      IMPLICIT NONE
345

        
346
C     ----- 
347
C     ARGUMENTS 
348
C     -----
349

        
350
      DOUBLEPRECISION PMIN, PMAX, TOL, PS
351

        
352
C     ------- 
353
C     COMMON BLOCKS 
354
C     -------
355

        
356
      DOUBLEPRECISION GAMMA 
357
      COMMON /ADIND/  GAMMA
358

        
359
      DOUBLEPRECISION RHOL, PL, UL, HL, CSL, VELL, WL, 
360
     &                RHOR, PR, UR, HR, CSR, VELR, WR 
361
      COMMON /STATES/ RHOL, PL, UL, HL, CSL, VELL, WL, 
362
     &                RHOR, PR, UR, HR, CSR, VELR, WR
363

        
364
C     --------- 
365
C     INTERNAL VARIABLES 
366
C     ---------
367

        
368
      DOUBLEPRECISION A, B, C, D, E, EPS, FA, FB, FC, TOL1, 
369
     & XM, P, Q, R, S
370

        
371
C     ------------- 
372
C     COMPUTE MACHINE PRECISION 
373
C     -------------
374

        
375
      EPS  = 1.D0 
376
10    EPS  = EPS/2.D0 
377
      TOL1 = 1.D0 + EPS 
378
      IF( TOL1 .GT. 1.D0 ) GO TO 10
379

        
380
C     ------- 
381
C     INITIALIZATION 
382
C     -------
383

        
384
      A  = PMIN 
385
      B  = PMAX 
386
      CALL GETDVEL(A,FA) 
387
      CALL GETDVEL(B,FB)
388

        
389
C     ----- 
390
C     BEGIN STEP 
391
C     -----
392

        
393
20    C  = A 
394
      FC = FA 
395
      D  = B - A 
396
      E  = D 
397
30    IF( DABS(FC) .GE. DABS(FB) )GO TO 40 
398
      A  = B 
399
      B  = C 
400
      C  = A 
401
      FA = FB 
402
      FB = FC 
403
      FC = FA
404

        
405
C     -------- 
406
C     CONVERGENCE TEST 
407
C     --------
408

        
409
40    TOL1 = 2.D0*EPS*DABS(B) + 0.5D0*TOL 
410
      XM   = 0.5D0*(C - B) 
411
      IF( DABS(XM) .LE. TOL1 ) GO TO 90 
412
      IF( FB .EQ. 0.D0 ) GO TO 90
413

        
414
C     ------------ 
415
C     IS BISECTION NECESSARY? 
416
C     ------------
417

        
418
      IF( DABS(E) .LT. TOL1 ) GO TO 70 
419
      IF( DABS(FA) .LE. DABS(FB) ) GO TO 70
420

        
421
C     ------------------ 
422
C     IS QUADRATIC INTERPOLATION POSSIBLE? 
423
C     ------------------
424

        
425
      IF( A .NE. C ) GO TO 50
426

        
427
C     ---------- 
428
C     LINEAR INTERPOLATION 
429
C     ----------
430

        
431
      S = FB/FA 
432
      P = 2.D0*XM*S 
433
      Q = 1.D0 - S 
434
      GO TO 60
435

        
436
C     ---------------- 
437
C     INVERSE QUADRATIC INTERPOLATION 
438
C     ----------------
439

        
440
50    Q = FA/FC 
441
      R = FB/FC 
442
      S = FB/FA 
443
      P = S*(2.D0*XM*Q*(Q - R) - (B - A)*(R - 1.D0)) 
444
      Q = (Q - 1.D0)*(R - 1.D0)*(S - 1.D0)
445

        
446
C     ------ 
447
C     ADJUST SIGNS 
448
C     ------
449

        
450
60    IF( P .GT. 0.D0 ) Q = -Q 
451
      P = DABS(P)
452

        
453
C     -------------- 
454
C     IS INTERPOLATION ACCEPTABLE? 
455
C     --------------
456

        
457
      IF( (2.D0*P) .GE. (3.D0*XM*Q-DABS(TOL1*Q)) ) GO TO 70 
458
      IF( P .GE. DABS(0.5D0*E*Q) ) GO TO 70 
459
      E = D 
460
      D = P/Q 
461
      GO TO 80
462

        
463
C     ----- 
464
C     BISECTION 
465
C     -----
466

        
467
70    D = XM 
468
      E = D
469

        
470
C     ------- 
471
C     COMPLETE STEP 
472
C     -------
473

        
474
80    A  = B 
475
      FA = FB 
476
      IF( DABS(D) .GT. TOL1 ) B = B+D 
477
      IF( DABS(D) .LE. TOL1 ) B = B+DSIGN(TOL1,XM) 
478
      CALL GETDVEL(B,FB) 
479
      IF( (FB*(FC/DABS(FC))) .GT. 0.D0) GO TO 20 
480
      GO TO 30
481

        
482
C     -- 
483
C     DONE 
484
C     --
485

        
486
90    PS = B
487

        
488
      RETURN 
489
      END
490

        
491
C     --------- 
492
CN    NAME: G E T V E L 
493
C     ---------
494

        
495
CP    PURPOSE: 
496
CP    COMPUTE THE FLOW VELOCITY BEHIND A RAREFACTION OR SHOCK IN TERMS OF THE 
497
CP    POST-WAVE PRESSURE FOR A GIVEN STATE AHEAD THE WAVE IN A RELATIVISTIC 
498
CP    FLOW 
499
C 
500

        
501
CC    COMMENTS: 
502
CC    THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN MARTI AND MUELLER, 
503
CC    J. FLUID MECH., (1994)
504

        
505
      SUBROUTINE GETVEL( P, RHOA, PA, UA, HA, CSA, VELA, WA, S, 
506
     & RHO,      U,  H,  CS,  VEL,  VSHOCK )
507

        
508
      IMPLICIT NONE
509

        
510
C     ----- 
511
C     ARGUMENTS 
512
C     -----
513

        
514
      DOUBLE PRECISION P, RHOA, PA, UA, HA, CSA, VELA, WA  
515
      CHARACTER*1      S 
516
      DOUBLE PRECISION RHO, U, H, CS, VEL, VSHOCK
517

        
518
C     ------- 
519
C     COMMON BLOCKS 
520
C     -------
521

        
522
      DOUBLE PRECISION GAMMA 
523
      COMMON /ADIND/   GAMMA
524

        
525
C     --------- 
526
C     INTERNAL VARIABLES 
527
C     ---------
528

        
529
      DOUBLE PRECISION A, B, C, SIGN 
530
      DOUBLE PRECISION J, WSHOCK 
531
      DOUBLE PRECISION K, SQGL1
532

        
533
C     --------------- 
534
C     LEFT OR RIGHT PROPAGATING WAVE 
535
C     ---------------
536

        
537
      IF (S.EQ.'L') SIGN = -1.D0
538

        
539
      IF (S.EQ.'R') SIGN =  1.D0
540

        
541
C
542

        
543
      IF (P.GT.PA) THEN
544

        
545
C       --- 
546
C       SHOCK 
547
C       ---
548

        
549
        A  = 1.D0+(GAMMA-1.D0)*(PA-P)/GAMMA/P 
550
        B  = 1.D0-A 
551
        C  = HA*(PA-P)/RHOA-HA**2
552

        
553
C       ---------------- 
554
C       CHECK FOR UNPHYSICAL ENTHALPIES 
555
C       ----------------
556

        
557
        IF (C.GT.(B**2/4.D0/A)) STOP 
558
     & 'GETVEL: UNPHYSICAL SPECIFIC ENTHALPY IN INTERMEDIATE STATE'
559

        
560
C       ----------------------------- 
561
C       SPECIFIC ENTHALPY IN THE POST-WAVE STATE 
562
C       (FROM THE EQUATION OF STATE AND THE TAUB ADIABAT, 
563
C       EQ.(74), MM94) 
564
C       -----------------------------
565

        
566
        H = (-B+DSQRT(B**2-4.D0*A*C))/2.D0/A
567

        
568
C       --------------- 
569
C       DENSITY IN THE POST-WAVE STATE 
570
C       (FROM EQ.(73), MM94) 
571
C       ---------------
572

        
573
        RHO = GAMMA*P/(GAMMA-1.D0)/(H-1.D0)
574

        
575
C       ------------------------ 
576
C       SPECIFIC INTERNAL ENERGY IN THE POST-WAVE STATE 
577
C       (FROM THE EQUATION OF STATE) 
578
C       ------------------------
579

        
580
        U = P/(GAMMA-1.D0)/RHO
581

        
582
C       -------------------------- 
583
C       MASS FLUX ACROSS THE WAVE  
584
C       (FROM THE RANKINE-HUGONIOT RELATIONS, EQ.(71), MM94) 
585
C       --------------------------
586

        
587
        J = SIGN*DSQRT((P-PA)/(HA/RHOA-H/RHO))
588

        
589
C       ---------- 
590
C       SHOCK VELOCITY 
591
C       (FROM EQ.(86), MM94 
592
C       ----------
593

        
594
        A      = J**2+(RHOA*WA)**2 
595
        B      = -VELA*RHOA**2*WA**2 
596
        VSHOCK = (-B+SIGN*J**2*DSQRT(1.D0+RHOA**2/J**2))/A 
597
        WSHOCK = 1.D0/DSQRT(1.D0-VSHOCK**2)
598

        
599
C       ------------------- 
600
C       FLOW VELOCITY IN THE POST-SHOCK STATE 
601
C       (FROM EQ.(67), MM94) 
602
C       -------------------
603

        
604
        A = WSHOCK*(P-PA)/J+HA*WA*VELA 
605
        B = HA*WA+(P-PA)*(WSHOCK*VELA/J+1.D0/RHOA/WA)
606

        
607
        VEL = A/B
608

        
609
C       --------------------- 
610
C       LOCAL SOUND SPEED IN THE POST-SHOCK STATE 
611
C       (FROM THE EQUATION OF STATE) 
612
C       ---------------------
613

        
614
        CS = DSQRT(GAMMA*P/RHO/H)
615

        
616
      ELSE
617

        
618
C       ------ 
619
C       RAREFACTION 
620
C       ------
621

        
622
C       --------------------------- 
623
C       POLITROPIC CONSTANT OF THE GAS ACROSS THE RAREFACTION 
624
C       ---------------------------
625

        
626
        K = PA/RHOA**GAMMA
627

        
628
C       --------------- 
629
C       DENSITY BEHIND THE RAREFACTION 
630
C       ---------------
631

        
632
        RHO = (P/K)**(1.D0/GAMMA)
633

        
634
C       ------------------------ 
635
C       SPECIFIC INTERNAL ENERGY BEHIND THE RAREFACTION 
636
C       (FROM THE EQUATION OF STATE) 
637
C       ------------------------
638

        
639
        U = P/(GAMMA-1.D0)/RHO
640

        
641
C       -------------------- 
642
C       LOCAL SOUND SPEED BEHIND THE RAREFACTION 
643
C       (FROM THE EQUATION OF STATE) 
644
C       --------------------
645

        
646
        CS = DSQRT(GAMMA*P/(RHO+GAMMA*P/(GAMMA-1.D0)))
647

        
648
C       ------------------ 
649
C       FLOW VELOCITY BEHIND THE RAREFACTION 
650
C       ------------------
651

        
652
        SQGL1 = DSQRT(GAMMA-1.D0) 
653
        A   = (1.D0+VELA)/(1.D0-VELA)* 
654
     & ((SQGL1+CSA)/(SQGL1-CSA)* 
655
     & (SQGL1-CS )/(SQGL1+CS ))**(-SIGN*2.D0/SQGL1)
656

        
657
        VEL = (A-1.D0)/(A+1.D0)
658

        
659
      END IF
660

        
661
      END
662

        
663
C     -------- 
664
CN    NAME: R A R E F 
665
C     --------
666

        
667
CP    PURPOSE: 
668
CP    COMPUTE THE FLOW STATE IN A RAREFACTION FOR GIVEN PRE-WAVE STATE 
669
C
670
 
671
CC    COMMENTS: 
672
CC    THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN MARTI AND MUELLER, 
673
CC    J. FLUID MECH., (1994)
674

        
675
      SUBROUTINE RAREF( XI, RHOA, PA, UA, CSA, VELA, S, RHO, P, U, VEL )
676

        
677
      IMPLICIT NONE
678

        
679
C     ----- 
680
C     ARGUMENTS 
681
C     -----
682

        
683
      DOUBLE PRECISION XI
684

        
685
      DOUBLE PRECISION RHOA, PA, UA, CSA, VELA
686

        
687
      CHARACTER        S
688

        
689
      DOUBLE PRECISION RHO, P, U, VEL
690

        
691
C     ------- 
692
C     COMMON BLOCKS 
693
C     -------
694

        
695
      DOUBLE PRECISION GAMMA 
696
      COMMON /ADIND/   GAMMA
697

        
698
C     --------- 
699
C     INTERNAL VARIABLES 
700
C     ---------
701

        
702
      DOUBLE PRECISION B, C, D, K, L, V, OCS2, FCS2, DFDCS2, CS2, SIGN
703

        
704
C     --------------- 
705
C     LEFT OR RIGHT PROPAGATING WAVE 
706
C     ---------------
707

        
708
      IF (S.EQ.'L') SIGN =  1.D0
709

        
710
      IF (S.EQ.'R') SIGN = -1.D0
711

        
712
      B    = DSQRT(GAMMA - 1.D0) 
713
      C    = (B + CSA)/(B - CSA) 
714
      D    = -SIGN*B/2.D0 
715
      K    = (1.D0 + XI)/(1.D0 - XI) 
716
      L    = C*K**D 
717
      V    = ((1.D0 - VELA)/(1.D0 + VELA))**D
718

        
719
      OCS2 = CSA
720

        
721
25    FCS2   = L*V*(1.D0 + SIGN*OCS2)**D*(OCS2 - B) + 
722
     & (1.D0 - SIGN*OCS2)**D*(OCS2 + B)
723

        
724
      DFDCS2 = L*V*(1.D0 + SIGN*OCS2)**D* 
725
     & (1.D0 + SIGN*D*(OCS2 - B)/(1.D0 + SIGN*OCS2)) + 
726
     & (1.D0 - SIGN*OCS2)**D*
727
     & (1.D0 - SIGN*D*(OCS2 + B)/(1.D0 - SIGN*OCS2))
728

        
729
      CS2 = OCS2 - FCS2/DFDCS2
730

        
731
      IF (ABS(CS2 - OCS2)/OCS2.GT.5.E-7)THEN 
732
        OCS2 = CS2 
733
        GOTO 25 
734
      END IF
735

        
736
      VEL = (XI + SIGN*CS2)/(1.D0 + SIGN*XI*CS2)
737

        
738
      RHO = RHOA*((CS2**2*(GAMMA - 1.D0 - CSA**2))/ 
739
     & (CSA**2*(GAMMA - 1.D0 - CS2**2))) 
740
     & **(1.D0/(GAMMA - 1.D0))
741

        
742
      P   = CS2**2*(GAMMA - 1.D0)*RHO/(GAMMA - 1.D0 - CS2**2)/GAMMA
743

        
744
      U   = P/(GAMMA - 1.D0)/RHO
745

        
746
      RETURN  
747
      END