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 P P M
3
C     --------
4

        
5
CP    PURPOSE: 
6
CP    SOLVE THE 1D EQUATIONS OF RELATIVISTIC HYDRODYNAMICS IN PLANAR SYMMETRY
7
CP    FOR AN IDEAL GAS EQUATION OF STATE.
8
C
9

        
10
CC    COMMENTS: 
11
CC    THIS PROGRAM IS DESCRIBED IN THE PAPER BY MARTI & MUELLER, JCP, 1996.
12
CC    IT USES AN EXACT RIEMANN SOLVER (MARTI & MUELLER, JFM, 1994), PPM SPATIAL
13
CC    RECONSTRUCTION AND AVERAGING IN THE DOMAIN OF DEPENDENCE OF THE 
14
CC    CELL INTERFACES FOR TIME ADVANCE.
15
CC    LIGHT SPEED IN CODE UNITS IS EQUAL TO 1.
16
CC 
17
CC    WRITTEN BY:     Jose-Maria Marti
18
CC                    Departamento de Astronomia y Astrofisica 
19
CC                    Universidad de Valencia 
20
CC                    46100 Burjassot (Valencia), Spain
21
CC                    jose-maria.marti@uv.es
22
CC    AND
23
CC                    Ewald Mueller
24
CC                    Max-Planck-Institut fuer Astrophysik
25
CC                    Karl-Schwarzschild-Str. 1
26
CC                    85741 Garching, Germany
27
CC                    emueller@mpa-garching.mpg.de
28
C
29

        
30
      PROGRAM RPPM
31

        
32
      IMPLICIT NONE
33

        
34
      INCLUDE 'size'
35

        
36
C     ----------
37
C     COMMON BLOCKS
38
C     ----------
39

        
40
      INTEGER         BNDMNX,BNDMXX
41
      COMMON /BOUN/   BNDMNX,BNDMXX
42

        
43
      INTEGER         NEND,NOUT,ITSTP,NX
44
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
45

        
46
      INTEGER         NSTEP
47
      COMMON /NSTEP/  NSTEP
48

        
49
      INTEGER         NOUT1
50
      COMMON /OUTI/   NOUT1
51

        
52
      DOUBLEPRECISION GAMMA
53
      COMMON /ADIND/  GAMMA
54

        
55
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)              
56
      COMMON /GRD/    X,XL,XR,DX
57

        
58
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),W(-4:MN5),
59
     &                U(-4:MN5),CS(-4:MN5),H(-4:MN5),DPDRH(-4:MN5),
60
     &                DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),E(-4:MN5)
61
      COMMON /HYDRO/  P,RHO,VEL,W,U,CS,H,DPDRH,DPDU,R,M,E
62

        
63
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
64
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
65
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
66
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
67

        
68
      DOUBLEPRECISION TOUT1
69
      COMMON /OUTF/   TOUT1
70

        
71
      DOUBLEPRECISION TIME,DT
72
      COMMON /ZEIT/   TIME,DT
73

        
74
      CHARACTER*7     OUTFIL
75
      CHARACTER*8     LABEL
76
      CHARACTER*4     BASENM
77
      CHARACTER*2     SUFFIX
78
      COMMON /CHRC/   LABEL,OUTFIL,BASENM,SUFFIX
79

        
80
C     --------------------
81
C     READ INITIAL PARAMETERS
82
C     --------------------
83

        
84
      CALL INPUT
85

        
86
C     --------------------
87
C     CONSTRUCT A NEW MODEL
88
C     --------------------
89

        
90
      IF (SUFFIX(2:2).NE.'A') THEN
91

        
92
        WRITE(6,2000)
93
 2000   FORMAT('RPPM: CHECK INPUT FILE SUFFIX')
94
        STOP
95

        
96
      END IF
97

        
98
      WRITE(6,2100)
99
 2100 FORMAT('RPPM: CONSTRUCTING NEW INITIAL MODEL')
100

        
101
      CALL GRID
102

        
103
      CALL INIT
104

        
105
      CALL TSTEP
106

        
107
      DT = MIN(DT, DTINI)
108

        
109
      OUTFIL = BASENM//'O'//SUFFIX
110

        
111
      NOUT1 = 0
112
      TOUT1 = 0.D0
113

        
114
C     -------------
115
C     START TIME LOOP
116
C     -------------
117

        
118
      DO 100 NSTEP = 1, NEND
119

        
120
         TIME  = TIME  + DT
121
         NOUT1 = NOUT1 + 1
122
         TOUT1 = TOUT1 + DT
123

        
124
         IF (TIME.GT.TMAX) GOTO 200
125

        
126
         CALL BNDRY
127

        
128
         CALL HYDROW
129

        
130
         IF ( (NOUT1.GE.NOUT) .OR. (TOUT1.GE.TOUT) ) CALL PLTOUT
131

        
132
         CALL TSTEP
133

        
134
 100     CONTINUE
135

        
136
 200  CONTINUE
137

        
138
      STOP 'RPPM: NORMAL TERMINATION'
139
      END
140

        
141
C     -------- 
142
CN    NAME: I N P U T
143
C     --------
144

        
145
CP    PURPOSE: 
146
CP    READS THE INITIAL PARAMETERS FROM FILE inpt.dat
147
C
148

        
149
CC    COMMENTS: 
150
CC    SEE INSERTED COMMENTS
151

        
152
      SUBROUTINE INPUT
153

        
154
      IMPLICIT NONE
155

        
156
      INCLUDE 'size'
157

        
158
C     ---------
159
C     COMMON BLOCKS
160
C     ---------
161

        
162
      INTEGER         NEND,NOUT,ITSTP,NX
163
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
164

        
165
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
166
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
167
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
168
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
169

        
170
      CHARACTER*7     OUTFIL
171
      CHARACTER*8     LABEL,LABEL1
172
      CHARACTER*4     BASENM
173
      CHARACTER*2     SUFFIX
174
      COMMON /CHRC/   LABEL,OUTFIL,BASENM,SUFFIX
175

        
176
C     ---------
177
C     INTERNAL VARIABLES
178
C     ---------
179

        
180
      CHARACTER*72 TEXT
181
      CHARACTER*2  TXT
182
      CHARACTER*15 TXTXT
183
      DATA TXTXT /'.............  '/
184

        
185
      OPEN(1,FILE='inpt.dat',FORM='FORMATTED',STATUS='OLD')
186

        
187
      READ (1,*) TEXT
188
      WRITE(6,*) TEXT
189
      READ (1,*) TEXT
190
      WRITE(6,*) TEXT
191
      READ (1,*) TEXT
192
      WRITE(6,*) TEXT
193
      PRINT*, '-------------------------------------------------------'
194

        
195
C     ----------------------------
196
C     BASENM IS THE ROOT FOR THE OUTPUT, PLOT AND RESTART FILE NAMES (ROOTS
197
C     RST_, RBW_, RSR_, RBWI STAND FOR SPECIAL TESTS)
198
C     ----------------------------
199

        
200
      READ (1,*) TXT,LABEL,BASENM
201
      WRITE(6,*) TXT,LABEL,TXTXT,BASENM
202
      LABEL1 = 'basenm'
203
      IF (LABEL.NE.LABEL1) GOTO 10
204

        
205
C     ----------------------------
206
C     NEND IS THE TOTAL NUMBER OF TIMESTEPS
207
C     ----------------------------
208

        
209
      READ (1,*) TXT,LABEL,NEND
210
      WRITE(6,*) TXT,LABEL,TXTXT,NEND
211
      LABEL1 = 'nend'
212
      IF (LABEL.NE.LABEL1) GOTO 10
213

        
214
C     ----------------------------
215
C     THE PROGRAM STOPS WHEN TIME IS .GE. THAN TMAX
216
C     ----------------------------
217
 
218
      READ (1,*) TXT,LABEL,TMAX
219
      WRITE(6,*) TXT,LABEL,TXTXT,TMAX
220
      LABEL1 = 'tmax'
221
      IF (LABEL.NE.LABEL1) GOTO 10
222

        
223
C     ----------------------------
224
C     SUFFIX IS THE SUFFIX FOR THE OUTPUT AND RESTART FILE NAMES
225
C     ----------------------------
226

        
227
      READ (1,*) TXT,LABEL,SUFFIX
228
      WRITE(6,*) TXT,LABEL,TXTXT,SUFFIX
229
      LABEL1 = 'suffix'
230
      IF (LABEL.NE.LABEL1) GOTO 10
231

        
232
C     ----------------------------
233
C     AN OUTPUT FILE IS WRITEN EVERY NOUT TIMESTEPS
234
C     ----------------------------
235

        
236
      READ (1,*) TXT,LABEL,NOUT
237
      WRITE(6,*) TXT,LABEL,TXTXT,NOUT
238
      LABEL1 = 'nout'
239
      IF (LABEL.NE.LABEL1) GOTO 10
240

        
241
C     ----------------------------
242
C     AN OUTPUT FILE IS WRITEN EVERY TOUT TIME UNITS
243
C     ----------------------------
244

        
245
      READ (1,*) TXT,LABEL,TOUT
246
      WRITE(6,*) TXT,LABEL,TXTXT,TOUT
247
      LABEL1 = 'tout'
248
      IF (LABEL.NE.LABEL1) GOTO 10
249

        
250
C     -----------------------------
251
C     DT IS WRITEN ON THE SCREEN EVERY ITSTP TIMESTEPS
252
C     -----------------------------
253

        
254
      READ (1,*) TXT,LABEL,ITSTP
255
      WRITE(6,*) TXT,LABEL,TXTXT,ITSTP
256
      LABEL1 = 'itstp'
257
      IF (LABEL.NE.LABEL1) GOTO 10
258

        
259
C     -----------------------------
260
C     CFL IS THE TIME-LIMITING FACTOR (.LT. 1)
261
C     -----------------------------
262

        
263
      READ (1,*) TXT,LABEL,CFL
264
      WRITE(6,*) TXT,LABEL,TXTXT,CFL
265
      LABEL1 = 'cfl'
266
      IF (LABEL.NE.LABEL1) GOTO 10
267

        
268
C     -------------------
269
C     DTINI IS THE INITIAL DT
270
C     -------------------
271

        
272
      READ (1,*) TXT,LABEL,DTINI
273
      WRITE(6,*) TXT,LABEL,TXTXT,DTINI
274
      LABEL1 = 'dtini'
275
      IF (LABEL.NE.LABEL1) GOTO 10
276

        
277
C     -------------------------------------------
278
C     SMALL IS THE THRESHOLD FOR NONDIMENSIONAL NUMBERS (I.E. VELOCITY)
279
C     -------------------------------------------
280

        
281
      READ (1,*) TXT,LABEL,SMALL
282
      WRITE(6,*) TXT,LABEL,TXTXT,SMALL
283
      LABEL1 = 'small'
284
      IF (LABEL.NE.LABEL1) GOTO 10 
285

        
286
C     -----------------------------------------
287
C     SMLRHO IS THE THRESHOLD FOR DENSITIES (RHO, R)
288
C     -----------------------------------------
289

        
290
      READ (1,*) TXT,LABEL,SMLRHO
291
      WRITE(6,*) TXT,LABEL,TXTXT,SMLRHO
292
      LABEL1 = 'smlrho'
293
      IF (LABEL.NE.LABEL1) GOTO 10
294

        
295
C     --------------------------------
296
C     SMALLP IS THE THRESHOLD FOR PRESSURE
297
C     --------------------------------
298

        
299
      READ (1,*) TXT,LABEL,SMALLP
300
      WRITE(6,*) TXT,LABEL,TXTXT,SMALLP
301
      LABEL1 = 'smallp'
302
      IF (LABEL.NE.LABEL1) GOTO 10
303

        
304
C     ---------------------------------
305
C     SMALLU IS THE THRESHOLD FOR INTERNAL ENERGY
306
C     ---------------------------------
307

        
308
      READ (1,*) TXT,LABEL,SMALLU
309
      WRITE(6,*) TXT,LABEL,TXTXT,SMALLU
310
      LABEL1 = 'smallu'
311
      IF (LABEL.NE.LABEL1) GOTO 10
312

        
313
C     ------------------------
314
C     GRIDLX IS THE LENGTH OF THE GRID
315
C     ------------------------
316

        
317
      READ (1,*) TXT,LABEL,GRIDLX
318
      WRITE(6,*) TXT,LABEL,TXTXT,GRIDLX
319
      LABEL1 = 'gridlx'
320
      IF (LABEL.NE.LABEL1) GOTO 10
321

        
322
      IF ((BASENM.EQ.'RST_'.OR.BASENM.EQ.'RBW_'.OR.BASENM.EQ.'RSR_'.OR.
323
     &     BASENM.EQ.'RBWI').AND.
324
     &     GRIDLX.NE.1.D0) THEN
325

        
326
         GRIDLX = 1.D0
327
         WRITE(6,1200)
328
 1200    FORMAT('INPUT: GRIDLX RESET TO 1.D0')
329

        
330
      END IF
331

        
332
C     ---------------------------
333
C     NX IS THE NUMBER OF GRID POINTS
334
C     ---------------------------
335

        
336
      READ (1,*) TXT,LABEL,NX
337
      WRITE(6,*) TXT,LABEL,TXTXT,NX
338
      LABEL1 = 'nx'
339
      IF (LABEL.NE.LABEL1) GOTO 10
340

        
341
      IF (NX.LT.4.OR.NX.GT.MNX) STOP 'INPUT: UNSUITABLE NX'
342

        
343
C     ----------------------------------------
344
C     ETA1 IS USED IN SUBROUTINE DETECT (PPM RECONSTRUCTION)
345
C     TYPICAL VALUE: 5.D0
346
C     ----------------------------------------
347

        
348
      READ (1,*) TXT,LABEL,ETA1
349
      WRITE(6,*) TXT,LABEL,TXTXT,ETA1
350
      LABEL1 = 'eta1'
351
      IF (LABEL.NE.LABEL1) GOTO 10
352

        
353
C     -------------------------------------------
354
C     ETA2 IS USED IN SUBROUTINE DETECT (PPM RECONSTRUCTION)
355
C     TYPICAL VALUE: 5.D-2
356
C     -------------------------------------------
357

        
358
      READ (1,*) TXT,LABEL,ETA2
359
      WRITE(6,*) TXT,LABEL,TXTXT,ETA2
360
      LABEL1 = 'eta2'
361
      IF (LABEL.NE.LABEL1) GOTO 10
362

        
363
C     ---------------------------------------
364
C     EPSLN IS USED IN SUBROUTINE DETECT (PPM RECONSTRUCTION)
365
C     TYPICAL VALUE: 1.D-1
366
C     ---------------------------------------
367

        
368
      READ (1,*) TXT,LABEL,EPSLN
369
      WRITE(6,*) TXT,LABEL,TXTXT,EPSLN
370
      LABEL1 = 'epsln'
371
      IF (LABEL.NE.LABEL1) GOTO 10
372

        
373
C     -----------------------------------------
374
C     AK0 IS USED IN SUBROUTINE DETECT (PPM RECONSTRUCTION)
375
C     TYPICAL VALUE: 1.D0
376
C     -----------------------------------------
377
      READ (1,*) TXT,LABEL,AK0
378
      WRITE(6,*) TXT,LABEL,TXTXT,AK0
379
      LABEL1 = 'ak0'
380
      IF (LABEL.NE.LABEL1) GOTO 10
381

        
382
C     ----------------------------------------
383
C     EPSILN IS USED IN SUBROUTINE FLATEN (PPM RECONSTRUCTION)
384
C     TYPICAL VALUE: 1.D0
385
C     ----------------------------------------
386

        
387
      READ (1,*) TXT,LABEL,EPSILN
388
      WRITE(6,*) TXT,LABEL,TXTXT,EPSILN
389
      LABEL1 = 'epsiln'
390
      IF (LABEL.NE.LABEL1) GOTO 10
391

        
392
C     ---------------------------------------
393
C     OMG1 IS USED IN SUBROUTINE FLATEN (PPM RECONSTRUCTION)
394
C     TYPICAL VALUE: 5.2D-1
395
C     ---------------------------------------
396

        
397
      READ (1,*) TXT,LABEL,OMG1
398
      WRITE(6,*) TXT,LABEL,TXTXT,OMG1
399
      LABEL1 = 'omg1'
400
      IF (LABEL.NE.LABEL1) GOTO 10
401

        
402
C     ----------------------------------------
403
C     OMG2 IS USED IN SUBROUTINE FLATEN (PPM RECONSTRUCTION)
404
C     TYPICAL VALUE: 1.D1
405
C     ----------------------------------------
406

        
407
      READ (1,*) TXT,LABEL,OMG2
408
      WRITE(6,*) TXT,LABEL,TXTXT,OMG2
409
      LABEL1 = 'omg2'
410
      IF (LABEL.NE.LABEL1) GOTO 10
411

        
412
      PRINT*, '-------------------------------------------------------'
413

        
414
      RETURN
415

        
416
 10   CONTINUE
417
      PRINT*, ' '
418
      WRITE(6,1020)
419
 1020 FORMAT('INPUT: INCORRECT INPUT DECK')
420
      WRITE(6,1001) LABEL, LABEL1
421
 1001 FORMAT('       LABEL = ',A6,'  EXPECTED LABEL = ',A6)
422
      STOP
423

        
424
      END
425

        
426
C     -------- 
427
CN    NAME: G R I D
428
C     --------
429

        
430
CP    PURPOSE: 
431
CP    ESTABLISHES THE BOUNDARY CONDITIONS AND SETS UP THE NUMERICAL GRID
432
C
433
 
434
CC    COMMENTS: 
435
CC    BOUNDARY CONDITIONS ARE SPECIFIED FOR A SERIES OF 1D TESTS. BOUNDARIES.
436
CC    APPROPRIATE BOUNDARIES MUST BE CHOSEN BY THE USER FOR SPECIFIC PROBLEMS.
437
CC    AN EQUIDISTANT X-GRID IS GENERATED BY DEFAULT, HOWEVER THE CODE IS
438
CC    SUITED FOR NON-EQUIDISTANT GRIDS
439

        
440
      SUBROUTINE GRID
441

        
442
      IMPLICIT NONE
443

        
444
      INCLUDE 'size'
445

        
446
C     ----------
447
C     COMMON BLOCKS
448
C     ----------
449

        
450
      INTEGER         BNDMNX,BNDMXX
451
      COMMON /BOUN/   BNDMNX,BNDMXX
452

        
453
      INTEGER         NEND,NOUT,ITSTP,NX
454
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
455

        
456
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
457
      COMMON /GRD/    X,XL,XR,DX
458

        
459
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
460
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
461
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
462
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
463

        
464
      CHARACTER*7     OUTFIL
465
      CHARACTER*8     LABEL
466
      CHARACTER*4     BASENM
467
      CHARACTER*2     SUFFIX
468
      COMMON /CHRC/   LABEL,OUTFIL,BASENM,SUFFIX
469

        
470
C     -------------
471
C     INTERNAL VARIABLES
472
C     -------------
473

        
474
      INTEGER         I
475

        
476
      DOUBLEPRECISION DELX
477

        
478
C     ---------------------------------
479
C     BOUNDARY CONDITIONS
480
C     BNDM.. = 1 ===> REFLECTING BOUNDARY
481
C     BNDM.. = 2 ===> FLOW OUT BOUNDARY
482
C     BNDM.. = 3 ===> FLOW IN BOUNDARY
483
C     BNDM.. = 4 ===> PERIODIC BOUNDARY
484
C     BNDM.. = 5 ===> ANY OTHER BOUNDARY
485
C     ---------------------------------
486

        
487
      IF (BASENM.EQ.'RST_'.OR.BASENM.EQ.'RBW_'.OR.BASENM.EQ.'RBWI') THEN
488
         BNDMNX = 2
489
         BNDMXX = 2
490
      ELSE IF (BASENM.EQ.'RSR_') THEN
491
         BNDMNX = 2
492
         BNDMXX = 1
493
      ELSE
494
         BNDMNX = 2
495
         BNDMXX = 2
496
      END IF
497

        
498
      IF (BNDMNX.EQ.4.AND.BNDMXX.NE.4) STOP 'GRID: INCORRECT BOUNDARIES'
499

        
500
C     ----------
501
C     SET UP X-GRID
502
C     ----------
503

        
504
      X(1) = 0.D0
505

        
506
      DELX = GRIDLX/DFLOAT(NX)
507

        
508
      DO 10 I=2,NX+1
509
         XL(I) = XL(I-1) + DELX
510
 10      CONTINUE
511

        
512

        
513
      DO 20 I=1,NX
514
         XR(I) = XL(I+1)
515
 20      CONTINUE
516

        
517
      DO 30 I=1,NX
518
         X(I) = 0.5D0*(XL(I) + XR(I))
519
 30      CONTINUE
520

        
521
      DO 40 I=1,NX
522
         DX(I) = XR(I) - XL(I)
523
 40      CONTINUE
524

        
525
      RETURN
526
      END
527

        
528
C     -------- 
529
CN    NAME: I N I T
530
C     --------
531

        
532
CP    PURPOSE: 
533
CP    DEFINES THE INITIAL MODEL
534
C
535
 
536
CC    COMMENTS: 
537
CC    DEFINES INITIAL DATA FOR A SERIES OF STANDARD 1D TEST PROBLEMS.
538
CC    APPROPRIATE INITIAL DATA MUST BE DEFINED BY THE USER FOR SPECIFIC 
539
CC    PROBLEMS.
540
   
541
      SUBROUTINE INIT
542

        
543
      IMPLICIT NONE
544

        
545
      INCLUDE 'size'
546

        
547
C     -----------
548
C     COMMON BLOCKS
549
C     -----------
550

        
551
      INTEGER         NEND,NOUT,ITSTP,NX
552
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
553

        
554
      DOUBLEPRECISION GAMMA
555
      COMMON /ADIND/  GAMMA
556

        
557
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
558
      COMMON /GRD/    X,XL,XR,DX
559

        
560
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),W(-4:MN5),
561
     &                U(-4:MN5),CS(-4:MN5),H(-4:MN5),DPDRH(-4:MN5),
562
     &                DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),E(-4:MN5)
563
      COMMON /HYDRO/  P,RHO,VEL,W,U,CS,H,DPDRH,DPDU,R,M,E
564

        
565
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
566
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
567
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
568
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
569

        
570
      DOUBLEPRECISION TIME,DT
571
      COMMON /ZEIT/   TIME,DT
572

        
573
      CHARACTER*7     OUTFIL
574
      CHARACTER*8     LABEL
575
      CHARACTER*4     BASENM
576
      CHARACTER*2     SUFFIX
577
      COMMON /CHRC/   LABEL,OUTFIL,BASENM,SUFFIX
578

        
579
C     -----------
580
C     INTERNAL VARIABLES
581
C     -----------
582

        
583
      INTEGER         I
584

        
585
C     ---------
586
C     INITIAL TIME
587
C     ---------
588

        
589
      TIME = 0.D0
590

        
591
C     ---------------------
592
C     RELATIVISTIC SOD'S TUBE
593
C     ---------------------
594

        
595
      IF (BASENM.EQ.'RST_') THEN
596

        
597
         GAMMA = 1.4D0
598
         
599
         DO 100 I=1,NX
600

        
601
            IF (X(I).LE.GRIDLX/2.D0) THEN
602
               VEL(I) = 0.D0
603
               RHO(I) = 1.D0
604
               U(I)   = 2.5D0
605
            ELSE
606
               VEL(I) = 0.D0
607
               RHO(I) = 0.125D0
608
               U(I)   = 2.D0
609
            END IF
610

        
611
 100        CONTINUE
612

        
613
         GOTO 450
614

        
615
      END IF
616

        
617
C     ----------------------
618
C     SCHNEIDER ET AL.'S TEST
619
C     ----------------------
620

        
621
      IF (BASENM.EQ.'SCHN') THEN
622

        
623
         GAMMA = 5.D0/3.D0
624
         
625
         DO 125 I=1,NX
626

        
627
            IF (X(I).LE.GRIDLX/2.D0) THEN
628
               VEL(I) = 0.D0
629
               RHO(I) = 10.D0
630
               U(I)   = 2.D0
631
            ELSE
632
               VEL(I) = 0.D0
633
               RHO(I) = 1.D0
634
               U(I)   = 1.D-6
635
            END IF
636

        
637
 125        CONTINUE
638

        
639
         GOTO 450
640

        
641
      END IF
642

        
643
C     --------------------
644
C     RELATIVISTIC BLAST WAVE
645
C     --------------------
646

        
647
      IF (BASENM.EQ.'RBW_') THEN
648

        
649
         GAMMA = 5.D0/3.D0
650

        
651
         DO 200 I=1,NX
652

        
653
            IF (X(I).LE.GRIDLX/2.D0) THEN
654
               VEL(I) = 0.D0
655
               RHO(I) = 1.D0
656
               U(I)   = 1.5D3
657
            ELSE
658
               VEL(I) = 0.D0
659
               RHO(I) = 1.D0
660
               U(I)   = 1.5D-2
661
            END IF
662

        
663
 200        CONTINUE
664

        
665
         GOTO 450
666

        
667
      END IF
668

        
669
C     -----------------------
670
C     RELATIVISTIC SHOCK REFLECTION
671
C     -----------------------
672

        
673
      IF (BASENM.EQ.'RSR_') THEN
674

        
675
         GAMMA = 4.D0/3.D0
676
         
677
         DO 300 I=1,NX
678

        
679
            VEL(I) = 0.99999D0
680
            RHO(I) = 1.D0
681
            U(I)   = 1.D-7/DSQRT(1.D0 - VEL(I)*VEL(I))
682

        
683
 300        CONTINUE
684

        
685
         GOTO 450
686

        
687
      END IF
688

        
689
C     ------------------------
690
C     RELATIVISTIC BLAST WAVE INTERACTION
691
C     ------------------------
692

        
693
      IF (BASENM.EQ.'RBWI') THEN
694

        
695
         GAMMA = 1.4D0
696
         
697
         DO 400 I=1,NX
698

        
699
            IF (X(I).LE.0.1D0*GRIDLX) THEN
700
               VEL(I) = 0.D0
701
               RHO(I) = 1.D0
702
               U(I)   = 2.5D3
703
            ELSE IF (X(I).LE.0.9D0*GRIDLX) THEN
704
               VEL(I) = 0.D0
705
               RHO(I) = 1.D0
706
               U(I)   = 2.5D-2
707
            ELSE
708
               VEL(I) = 0.D0
709
               RHO(I) = 1.D0
710
               U(I)   = 2.5D2
711
            END IF
712

        
713
 400        CONTINUE
714

        
715
        GOTO 450
716
      END IF
717
      
718
      STOP 'INIT: NO INITIAL DATA SPECIFIED'
719

        
720
 450  CONTINUE
721

        
722
      CALL EOS (NX, RHO, U, GAMMA, P, H, CS, DPDRH, DPDU)
723

        
724
      DO 500 I=1,NX
725

        
726
            W(I)     = 1.D0/DSQRT(1.D0 - VEL(I)*VEL(I))
727
      
728
            R(I)     = RHO(I)*W(I)
729
            M(I)     = R(I)*H(I)*W(I)*VEL(I)
730
            E(I)     = R(I)*H(I)*W(I) - P(I) - R(I)
731

        
732
 500        CONTINUE
733

        
734
      RETURN
735
      END
736

        
737
C     -------- 
738
CN    NAME: T S T E P 
739
C     --------
740

        
741
CP    PURPOSE: 
742
CP    COMPUTES THE NEW TIMESTEP VALUE FROM COURANT CONDITION
743
C
744
 
745
CC    COMMENTS: 
746
CC    NONE
747

        
748
      SUBROUTINE TSTEP
749

        
750
      IMPLICIT NONE
751

        
752
      INCLUDE 'size'
753

        
754
C     -----------
755
C     COMMON BLOCKS
756
C     -----------
757

        
758
      INTEGER         NEND,NOUT,ITSTP,NX
759
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
760

        
761
      INTEGER         NSTEP
762
      COMMON /NSTEP/  NSTEP
763

        
764
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
765
      COMMON /GRD/    X,XL,XR,DX
766

        
767
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),W(-4:MN5),
768
     &                U(-4:MN5),CS(-4:MN5),H(-4:MN5),DPDRH(-4:MN5),
769
     &                DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),E(-4:MN5)
770
      COMMON /HYDRO/  P,RHO,VEL,W,U,CS,H,DPDRH,DPDU,R,M,E
771

        
772
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
773
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
774
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
775
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
776

        
777
      DOUBLEPRECISION TIME,DT
778
      COMMON /ZEIT/   TIME,DT
779

        
780
C     ------------
781
C     INTERNAL VARIABLES
782
C     ------------
783

        
784
      INTEGER         I,IC
785

        
786
      DOUBLEPRECISION DTCC,DTEST(MN)
787

        
788
      DOUBLEPRECISION LAMBD1,LAMBD4,LAMBDA
789

        
790
      DOUBLEPRECISION V
791

        
792
      IC   = 0
793
      DTCC = 0.D0
794

        
795
      DO 10 I=1,NX
796

        
797
         LAMBD1   = (VEL(I) - CS(I))/(1.D0 - VEL(I)*CS(I))
798
         LAMBD4   = (VEL(I) + CS(I))/(1.D0 + VEL(I)*CS(I))
799
         LAMBDA   = DMAX1(DABS(LAMBD1),DABS(LAMBD4))
800
         DTEST(I) = LAMBDA/(XR(I) - XL(I))
801

        
802
 10      CONTINUE
803

        
804
      DO 13 I=1,NX
805

        
806
         IF (DTEST(I).GT.DTCC) THEN
807
            IC   = I
808
            DTCC = DTEST(I)
809
         END IF
810

        
811
 13      CONTINUE
812

        
813
      DT = CFL/DTCC
814

        
815
      V  = DABS(VEL(IC))
816

        
817
      IF (MOD(NSTEP,ITSTP).NE.0) RETURN
818

        
819
      WRITE(6,1001) NSTEP,DT,IC,CS(IC),V
820

        
821
 1001 FORMAT(I5,2X,1PE8.1,2X,I5,2X,1P2E11.2,2X,1P2E11.2)
822

        
823
      RETURN
824
      END
825

        
826
C     -------- 
827
CN    NAME: B N D R Y
828
C     --------
829

        
830
CP    PURPOSE: 
831
CP    PROVIDES DIFFERENT TYPES OF BOUNDARY CONDITIONS
832
C
833

        
834
CC    COMMENTS: 
835
CC    NEW BOUNDARY CONDITIONS CAN BE SPECIFIED AT THE USER'S WILL
836

        
837
      SUBROUTINE BNDRY
838

        
839
      IMPLICIT NONE
840

        
841
      INCLUDE 'size'
842

        
843
C     --------
844
C     COMMON BLOCKS
845
C     --------
846

        
847
      INTEGER         BNDMNX,BNDMXX
848
      COMMON /BOUN/   BNDMNX,BNDMXX
849

        
850
      INTEGER         NEND,NOUT,ITSTP,NX
851
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
852

        
853
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
854
      COMMON /GRD/    X,XL,XR,DX
855

        
856
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),W(-4:MN5),
857
     &                U(-4:MN5),CS(-4:MN5),H(-4:MN5),DPDRH(-4:MN5),
858
     &                DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),E(-4:MN5)
859
      COMMON /HYDRO/  P,RHO,VEL,W,U,CS,H,DPDRH,DPDU,R,M,E
860

        
861
      DOUBLEPRECISION TIME,DT
862
      COMMON /ZEIT/   TIME,DT
863
      
864
C     ---------
865
C     INTERNAL VARIABLES
866
C     ---------
867

        
868
      INTEGER         I
869

        
870
C     ------------
871
C     LEFT BOUNDARY
872
C     ------------
873

        
874
      GOTO (410, 420, 430, 440, 450), BNDMNX
875

        
876
C     -----------------
877
C     REFLECTING BOUNDARY
878
C     -----------------
879

        
880
 410  CONTINUE
881

        
882
      DO 415 I=-4,0
883

        
884
         RHO(I)  = RHO(1-I)
885
         VEL(I)  = -VEL(1-I)
886
         U(I)    = U(1-I)
887
         P(I)    = P(1-I)
888
         CS(I)   = CS(1-I)
889
         H(I)    = H(1-I)
890
         W(I)    = W(1-I)
891
         DX(I)   = DX(1-I)
892

        
893
 415     CONTINUE
894

        
895
      GOTO 500
896

        
897
C     ---------------
898
C     FLOW OUT BOUNDARY
899
C     ---------------
900

        
901
 420  CONTINUE
902

        
903
      DO 425 I=-4,0
904

        
905
         RHO(I)  = RHO(1)
906
         VEL(I)  = VEL(1)
907
         U(I)    = U(1)
908
         P(I)    = P(1)
909
         CS(I)   = CS(1)
910
         H(I)    = H(1)
911
         W(I)    = W(1)
912
         DX(I)   = DX(1)
913

        
914
 425     CONTINUE
915

        
916
      GOTO 500
917

        
918
C     ---------------
919
C     FLOW IN BOUNDARY
920
C     ---------------
921

        
922
 430  CONTINUE
923

        
924
      STOP 'BNDRY: INFLOW BOUNDARY MUST BE SUPPLIED BY THE USER'
925

        
926
      GOTO 500
927

        
928
C     --------------
929
C     PERIODIC BOUNDARY
930
C     --------------
931

        
932
 440  CONTINUE
933

        
934
      DO 445 I=-4,0
935

        
936
         RHO(I)  = RHO(NX+I)
937
         VEL(I)  = VEL(NX+I)
938
         U(I)    = U(NX+I)
939
         P(I)    = P(NX+I)
940
         CS(I)   = CS(NX+I)
941
         H(I)    = H(NX+I)
942
         W(I)    = W(NX+I)
943
         DX(I)   = DX(NX+I)
944

        
945
 445     CONTINUE
946

        
947
      GOTO 500
948

        
949
C     ---------------------------------------------------------
950
C     SPECIAL BOUNDARY (ADD ANY NONSTANDARD BOUNDARY CONDITION HERE)
951
C     INFLOW JET BOUNDARY
952
C     ---------------------------------------------------------
953

        
954
 450  CONTINUE
955

        
956
      STOP 'BNDRY: NON STANDARD BOUNDARY. TO BE SUPPLIED BY THE USER'
957
            
958
 500  CONTINUE
959

        
960
      DO 505 I=0,-4,-1
961

        
962
         XL(I) = XL(I+1)-DX(I)
963
         XR(I) = XR(I+1)-DX(I+1)
964
         X(I)  = 0.5*(XL(I)+XR(I))
965

        
966
 505     CONTINUE
967

        
968
C     ------------
969
C     RIGHT BOUNDARY
970
C     ------------
971

        
972
      GOTO (510, 520, 530, 540, 550), BNDMXX
973

        
974
C     -----------------
975
C     REFLECTING BOUNDARY
976
C     -----------------
977

        
978
 510  CONTINUE
979

        
980
      DO 515 I=1,5
981

        
982
         RHO(NX+I)  = RHO(NX+1-I)
983
         VEL(NX+I)  = -VEL(NX+1-I)
984
         U(NX+I)    = U(NX+1-I)
985
         P(NX+I)    = P(NX+1-I)
986
         CS(NX+I)   = CS(NX+1-I)
987
         H(NX+I)    = H(NX+1-I)
988
         W(NX+I)    = W(NX+1-I)
989
         DX(NX+I)   = DX(NX+1-I)
990

        
991
 515     CONTINUE
992

        
993
      GOTO 600
994

        
995
C     --------------
996
C     FLOW OUT BOUNDARY
997
C     --------------
998

        
999
 520  CONTINUE
1000

        
1001
      DO 525 I=NX+1,NX+5
1002

        
1003
         RHO(I)  = RHO(NX)
1004
         VEL(I)  = VEL(NX)
1005
         U(I)    = U(NX)
1006
         P(I)    = P(NX)
1007
         CS(I)   = CS(NX)
1008
         H(I)    = H(NX)
1009
         W(I)    = W(NX)
1010
         DX(I)   = DX(NX)
1011

        
1012
 525     CONTINUE
1013

        
1014
      GOTO 600
1015

        
1016
C     --------------
1017
C     FLOW IN BOUNDARY
1018
C     --------------
1019

        
1020
 530  CONTINUE
1021

        
1022
      STOP 'BNDRY: INFLOW BOUNDARY MUST BE SUPPLIED BY THE USER'
1023

        
1024
      GOTO 600
1025

        
1026
C     ---------------
1027
C     PERIODIC BOUNDARY
1028
C     ---------------
1029

        
1030
 540  CONTINUE
1031

        
1032
      DO 545 I=1,5
1033

        
1034
         RHO(NX+I)  = RHO(I)
1035
         VEL(NX+I)  = VEL(I)
1036
         U(NX+I)    = U(I)
1037
         P(NX+I)    = P(I)
1038
         CS(NX+I)   = CS(I)
1039
         H(NX+I)    = H(I)
1040
         W(NX+I)    = W(I)
1041
         DX(NX+I)   = DX(I)
1042

        
1043
 545     CONTINUE
1044

        
1045
      GOTO 600
1046

        
1047
C     ---------------------------------------
1048
C     SPECIAL BOUNDARY
1049
C     ADD ANY NONSTANDARD BOUNDARY CONDITION HERE
1050
C     ---------------------------------------
1051

        
1052
 550  CONTINUE
1053

        
1054
      DO 555 I=NX+1,NX+5
1055

        
1056
         DX(I)   = DX(NX)
1057
         XL(I)   = XL(I-1) + DX(I-1)
1058
         XR(I)   = XR(I-1) + DX(I)
1059
         X(I)    = 0.5*(XL(I)+XR(I))
1060
         RHO(I)  = 1.D0+DABS(VEL(NX))*TIME/X(I)
1061
         VEL(I)  = VEL(NX)
1062
         U(I)    = U(NX)
1063
         P(I)    = P(NX)
1064
         CS(I)   = CS(NX)
1065
         H(I)    = H(NX)
1066
         W(I)    = W(NX)
1067

        
1068
 555     CONTINUE
1069

        
1070
 600  CONTINUE
1071

        
1072
      DO 650 I=NX+1,NX+5
1073

        
1074
         XL(I) = XL(I-1) + DX(I-1)
1075
         XR(I) = XR(I-1) + DX(I)
1076
         X(I)  = 0.5*(XL(I)+XR(I))
1077

        
1078
 650     CONTINUE
1079

        
1080
      RETURN
1081
      END
1082

        
1083
C     -------- 
1084
CN    NAME: H Y D R O W
1085
C     --------
1086

        
1087
CP    PURPOSE: 
1088
CP    ADVANCE IN TIME THE  1D EQUATIONS OF RELATIVISTIC HYDRODYNAMICS 
1089
CP    (IN CONSERVATION FORM) IN PLANAR COORDINATES.
1090
C
1091

        
1092
CC    COMMENTS: 
1093
CC    NONE
1094

        
1095
      SUBROUTINE HYDROW
1096

        
1097
      IMPLICIT NONE
1098

        
1099
      INCLUDE 'size'
1100

        
1101
C     ---------
1102
C     COMMON BLOCKS
1103
C     ---------
1104

        
1105
      INTEGER         NEND,NOUT,ITSTP,NX
1106
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
1107

        
1108
      DOUBLEPRECISION COEFF1(-4:MN5),COEFF2(-4:MN5),COEFF3(-4:MN5),
1109
     &                COEFF4(-4:MN5),COEFF5(-4:MN5)
1110
      COMMON /COEFF/  COEFF1,COEFF2,COEFF3,COEFF4,COEFF5
1111

        
1112
      DOUBLEPRECISION DELP(-4:MN5),DELRHO(-4:MN5),DELVEL(-4:MN5),
1113
     &                DELU(-4:MN5)
1114
      COMMON /DELU/   DELP,DELRHO,DELVEL,DELU
1115

        
1116
      DOUBLEPRECISION FICT(-4:MN5)
1117
      COMMON /FICT/   FICT
1118

        
1119
      DOUBLEPRECISION FLATN(-4:MN5),FLATN1(-4:MN5)
1120
      COMMON /FLAT/   FLATN,FLATN1
1121

        
1122
      DOUBLEPRECISION GAMMA
1123
      COMMON /ADIND/  GAMMA
1124
      
1125
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
1126
      COMMON /GRD/    X,XL,XR,DX
1127

        
1128
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),
1129
     &                W(-4:MN5),U(-4:MN5),CS(-4:MN5),H(-4:MN5),
1130
     &                DPDRH(-4:MN5),DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),
1131
     &                E(-4:MN5)
1132
      COMMON /HYDRO/  P,RHO,VEL,W,U,CS,H,DPDRH,DPDU,R,M,E
1133

        
1134
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,
1135
     &                SMALLU,GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
1136
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,
1137
     &                SMALLU,GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
1138

        
1139
      DOUBLEPRECISION PL(-4:MN6),PR(-4:MN6),RHOL(-4:MN6),RHOR(-4:MN6),
1140
     &                VELL(-4:MN6),VELR(-4:MN6),
1141
     &                UL(-4:MN6),UR(-4:MN6),CSL(-4:MN6),
1142
     &                CSR(-4:MN6),RL(-4:MN6),RR(-4:MN6),ML(-4:MN6),
1143
     &                MR(-4:MN6),EL(-4:MN6),
1144
     &                ER(-4:MN6)
1145
      COMMON /INTERF/ PL,PR,RHOL,RHOR,VELL,VELR,UL,UR,CSL,
1146
     &                CSR,RL,RR,ML,MR,EL,ER
1147

        
1148
      DOUBLEPRECISION TIME,DT
1149
      COMMON /ZEIT/   TIME,DT
1150

        
1151
      DOUBLEPRECISION DP(-4:MN5),P6(-4:MN5),DRHO(-4:MN5),RHO6(-4:MN5),
1152
     &                DVEL(-4:MN5),VEL6(-4:MN5),
1153
     &                DU(-4:MN5),U6(-4:MN5)
1154
      COMMON /U6/     DP,P6,DRHO,RHO6,DVEL,VEL6,DU,U6
1155

        
1156
      DOUBLEPRECISION PM(-4:MN5),PP(-4:MN5),RHOM(-4:MN5),RHOP(-4:MN5),
1157
     &                VELM(-4:MN5),VELP(-4:MN5),
1158
     &                UM(-4:MN5),UP(-4:MN5)
1159
      COMMON /UMP/    PM,PP,RHOM,RHOP,VELM,VELP,UM,UP
1160

        
1161
C     ------------
1162
C     INTERNAL VARIABLES
1163
C     ------------
1164

        
1165
      INTEGER         I
1166

        
1167
      DOUBLEPRECISION AUX1
1168

        
1169
      DOUBLEPRECISION DTDX(-4:MN5)
1170

        
1171
      DOUBLEPRECISION RFLX(-4:MN6),MFLX(-4:MN6),EFLX(-4:MN6)
1172

        
1173
C     -------------
1174
C     SPATIAL RECONSTRUCTION (PPM)
1175
C     -------------
1176

        
1177
      CALL COEF(NX,COEFF1,COEFF2,COEFF3,COEFF4,COEFF5)
1178

        
1179
      CALL INTERP(NX,PM   ,P   ,PP   ,DELP  )
1180

        
1181
      CALL INTERP(NX,RHOM ,RHO ,RHOP ,DELRHO)
1182

        
1183
      CALL DETECT(NX,RHOM ,RHO ,RHOP ,DELRHO)
1184

        
1185
      CALL INTERP(NX,VELM ,VEL ,VELP ,DELVEL)
1186

        
1187
      CALL FLATEN
1188

        
1189

        
1190
      DO 17 I=0,NX+1
1191
         RHOM(I)  = FLATN(I)*RHO(I) + FLATN1(I)*RHOM(I)
1192
         RHOP(I)  = FLATN(I)*RHO(I) + FLATN1(I)*RHOP(I)
1193
         VELM(I)  = FLATN(I)*VEL(I) + FLATN1(I)*VELM(I)
1194
         VELP(I)  = FLATN(I)*VEL(I) + FLATN1(I)*VELP(I)
1195
         PM(I)    = FLATN(I)*P(I)   + FLATN1(I)*PM(I)
1196
         PP(I)    = FLATN(I)*P(I)   + FLATN1(I)*PP(I)
1197
 17      CONTINUE
1198

        
1199
      CALL MONOT(NX,PM   ,P   ,PP   ,DP   ,P6   )
1200

        
1201
      CALL MONOT(NX,RHOM ,RHO ,RHOP ,DRHO ,RHO6 )
1202

        
1203
      CALL MONOT(NX,VELM ,VEL ,VELP ,DVEL ,VEL6 )
1204

        
1205
C     -------------------
1206
C     AVERAGED STATES FOR TIME ADVANCE
1207
C     -------------------
1208

        
1209
      CALL STAT1D
1210

        
1211
C     ----------------
1212
C     COMPUTATION OF NUMERICAL FLUXES (WITH AN EXACT RIEMANN SOLVER)
1213
C     ----------------
1214

        
1215
      DO 23 I=1,NX+1 
1216

        
1217
         CALL NFLUX(RHOL(I),RHOR(I),PL(I),PR(I),VELL(I),VELR(I),
1218
     &              UL(I),UR(I),CSL(I),CSR(I),RFLX(I),MFLX(I),EFLX(I))
1219
 23      CONTINUE
1220
        
1221
C     ----------
1222
C     TIME ADVANCE
1223
C     ----------
1224

        
1225
      DO 50 I=1,NX
1226
         DTDX(I) = DT/DX(I)
1227
         R(I)    = R(I) - DTDX(I)*(RFLX(I+1) - RFLX(I))
1228
         R(I)    = DMAX1(SMLRHO,R(I))
1229
 50      CONTINUE
1230

        
1231
      DO 60 I=1,NX
1232
         AUX1  = -DTDX(I)*(MFLX(I+1) - MFLX(I) )
1233
         M(I)  = M(I) + AUX1 + DT*FICT(I)
1234
         AUX1  = -DTDX(I)*(EFLX(I+1) - EFLX(I) )
1235
         E(I)  = E(I) + AUX1
1236
 60      CONTINUE
1237

        
1238
C     --------------
1239
C     PRIMITIVE VARIABLES
1240
C     --------------
1241

        
1242
      CALL GETPRFQ(NX,R,M,E,VEL,W,RHO,U,P,H,CS,DPDRH,DPDU)
1243

        
1244
      RETURN
1245
      END
1246

        
1247
C     -------- 
1248
CN    NAME: P L T O U T
1249
C     --------
1250

        
1251
CP    PURPOSE: 
1252
CP    COMPUTES THE ANALYTICAL SOLUTION OF STANDARD 1D PROBLEMS AND WRITES
1253
CP    THE RESULTS IN THE OUTPUT FILES
1254
C
1255
 
1256
CC    COMMENTS: 
1257
CC    NONE
1258

        
1259
      SUBROUTINE PLTOUT
1260

        
1261
      IMPLICIT NONE  
1262

        
1263
      INCLUDE 'size'  
1264

        
1265
C     -----------
1266
C     COMMON BLOCKS
1267
C     -----------
1268

        
1269
      INTEGER         NSTEP  
1270
      COMMON /NSTEP/  NSTEP  
1271

        
1272
      INTEGER         NOUT1
1273
      COMMON /OUTI/   NOUT1  
1274

        
1275
      INTEGER         NEND,NOUT,ITSTP,NX
1276
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
1277
 
1278
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),W(-4:MN5),
1279
     &                U(-4:MN5),CS(-4:MN5),H(-4:MN5),DPDRH(-4:MN5),
1280
     &                DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),E(-4:MN5)
1281
      COMMON /HYDRO/  P,RHO,VEL,W,U,CS,H,DPDRH,DPDU,R,M,E
1282

        
1283
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
1284
      COMMON /GRD/    X,XL,XR,DX
1285

        
1286
      DOUBLEPRECISION TOUT1
1287
      COMMON /OUTF/   TOUT1  
1288

        
1289
      DOUBLEPRECISION TIME,DT  
1290
      COMMON /ZEIT/   TIME,DT  
1291

        
1292
      DOUBLEPRECISION G  
1293
      COMMON /ADIND/  G  
1294

        
1295
      CHARACTER*7     OUTFIL  
1296
      CHARACTER*8     LABEL  
1297
      CHARACTER*4     BASENM  
1298
      CHARACTER*2     SUFFIX  
1299
      COMMON /CHRC/   LABEL,OUTFIL,BASENM,SUFFIX  
1300

        
1301
C     ------------
1302
C     INTERNAL VARIABLES
1303
C     ------------
1304

        
1305
      INTEGER         I  
1306

        
1307
      DOUBLEPRECISION P1,RHO1,VEL1,U1,CS1  
1308

        
1309
      DOUBLEPRECISION OCS2,CS2,FCS2,DFDCS2,CS9  
1310

        
1311
      DOUBLEPRECISION P2,RHO2,VEL2,U2  
1312

        
1313
      DOUBLEPRECISION P3,RHO3,VEL3,U3,CS3  
1314

        
1315
      DOUBLEPRECISION P4,RHO4,VEL4,U4,CS4  
1316

        
1317
      DOUBLEPRECISION P5,RHO5,VEL5,U5,CS5  
1318

        
1319
      DOUBLEPRECISION P6,RHO6,VEL6,U6,CS6  
1320

        
1321
      DOUBLEPRECISION P7,RHO7,VEL7,U7,CS7  
1322

        
1323
      DOUBLEPRECISION P8,RHO8,VEL8,U8,CS8  
1324

        
1325
      DOUBLEPRECISION P10,RHO10,VEL10,U10,CS10  
1326

        
1327
      DOUBLEPRECISION X1,X2,X3,X4,X5,X6,X7,X8,X9  
1328

        
1329
      DOUBLEPRECISION A,B,C,D,K,L  
1330

        
1331
      DOUBLEPRECISION PA(MN),RHOA(MN),VELA(MN),UA(MN)  
1332

        
1333
      DOUBLEPRECISION WW,VS,XS  
1334

        
1335
C     ----------------------
1336
C     ANALYTICAL SOLUTION FOR STANDARD TESTS
1337
C     ----------------------
1338

        
1339
      IF (BASENM.EQ.'RSR_') THEN  
1340

        
1341
        WW   = 1.D0/SQRT(1.D0 - VEL(1)*VEL(1))  
1342
        U1   = 0.D0  
1343
        U2   = WW - 1.D0  
1344
        RHO1 = 1.D0  
1345
        RHO2 = ((G + 1.D0)/(G - 1.D0) + G*U2/(G - 1.D0))*RHO1  
1346
        VEL1 = VEL(1)  
1347
        VEL2 = 0.D0  
1348
        P1   = (G - 1.D0)*RHO1*U1  
1349
        P2   = (G - 1.D0)*RHO2*U2  
1350
        VS   = -VEL1/(RHO2/RHO1/WW - 1.D0)  
1351
        XS   = XR(NX) + VS*TIME
1352
  
1353
        DO 30 I = 1,NX  
1354

        
1355
          IF (X(I).LT.XS) THEN  
1356
            UA(I)   = U1  
1357
            RHOA(I) = RHO1  
1358
            VELA(I) = VEL1  
1359
            PA(I)   = P1  
1360
          ELSE  
1361
            UA(I)   = U2  
1362
            RHOA(I) = RHO2  
1363
            VELA(I) = VEL2  
1364
            PA(I)   = P2  
1365
          END IF  
1366

        
1367
 30       CONTINUE  
1368

        
1369
      END IF  
1370
 
1371
      IF (BASENM.EQ.'RST_') THEN
1372
  
1373
        P1   = 1.D0
1374
        RHO1 = 1.D0  
1375
        VEL1 = 0.D0  
1376
        U1   = P1/(G - 1.D0)/RHO1  
1377
        CS1  = DSQRT(G*P1*(G - 1.D0)/(RHO1*(G - 1.D0) + G*P1))  
1378
 
1379
        P3   = 0.3115D0  
1380
        RHO3 = 0.4345D0  
1381
        VEL3 = 0.4262D0  
1382
        U3   = P3/(G - 1.D0)/RHO3  
1383
        CS3  = DSQRT(G*P3*(G - 1.D0)/(RHO3*(G - 1.D0) + G*P3))  
1384

        
1385
        P4   = P3  
1386
        RHO4 = 0.273D0  
1387
        VEL4 = VEL3  
1388
        U4   = P4/(G - 1.D0)/RHO4  
1389
        CS4  = DSQRT(G*P4*(G - 1.D0)/(RHO4*(G - 1.D0) + G*P4))  
1390
 
1391
        P5   = 0.1D0  
1392
        RHO5 = 0.125D0  
1393
        VEL5 = 0.D0  
1394
        U5   = P5/(G - 1.D0)/RHO5  
1395
        CS5  = DSQRT(G*P5*(G - 1.D0)/(RHO5*(G - 1.D0) + G*P5))  
1396
 
1397
        X1 = 0.5D0 + (VEL1 - CS1)*TIME/(1.D0 - VEL1*CS1)  
1398
 
1399
        X2 = 0.5D0 + (VEL3 - CS3)*TIME/(1.D0 - VEL3*CS3)  
1400
 
1401
        X3 = 0.5D0 + VEL3*TIME  
1402
 
1403
        X4 = 0.5D0 + VEL4*TIME/(1.D0 - RHO5*DSQRT(1.D0 - VEL4**2)/RHO4)  
1404
 
1405
      ELSE IF(BASENM.EQ.'RBW_') THEN  
1406
 
1407
        P1   = 1000.D0  
1408
        RHO1 = 1.D0  
1409
        VEL1 = 0.D0  
1410
        U1   = P1/(G-1.D0)/RHO1  
1411
        CS1  = DSQRT(G*P1*(G - 1.D0)/(RHO1*(G - 1.D0) + G*P1))  
1412
 
1413
        P3   = 18.6D0  
1414
        RHO3 = 9.15D-2  
1415
        VEL3 = 0.960D0  
1416
        U3   = P3/(G-1.D0)/RHO3  
1417
        CS3  = DSQRT(G*P3*(G - 1.D0)/(RHO3*(G - 1.D0) + G*P3))  
1418
 
1419
        P4   = P3  
1420
        RHO4 = 10.75D0  
1421
        VEL4 = VEL3  
1422
        U4   = P4/(G - 1.D0)/RHO4  
1423
        CS4  = DSQRT(G*P4*(G - 1.D0)/(RHO4*(G - 1.D0) + G*P4))  
1424
 
1425
        P5   = 1.D-2  
1426
        RHO5 = 1.D0  
1427
        VEL5 = 0.D0  
1428
        U5   = P5/(G - 1.D0)/RHO5  
1429
        CS5  = DSQRT(G*P5*(G - 1.D0)/(RHO5*(G - 1.D0) + G*P5))  
1430
 
1431
        X1 = 0.5D0 + (VEL1 - CS1)*TIME/(1.D0 - VEL1*CS1)  
1432
 
1433
        X2 = 0.5D0 + (VEL3 - CS3)*TIME/(1.D0 - VEL3*CS3)  
1434
 
1435
        X3 = 0.5D0 + VEL3*TIME  
1436
 
1437
        X4 = 0.5D0 + VEL4*TIME/(1.D0 - RHO5*DSQRT(1.D0 - VEL4**2)/RHO4)  
1438

        
1439
      ELSE IF(BASENM.EQ.'SCHN') THEN  
1440

        
1441
        P1   = 13.33333333D0  
1442
        RHO1 = 10.D0  
1443
        VEL1 = 0.D0  
1444
        U1   = P1/(G - 1.D0)/RHO1  
1445
        CS1  = DSQRT(G*P1*(G - 1.D0)/(RHO1*(G - 1.D0) + G*P1))  
1446
 
1447
        P3   = 1.448D0  
1448
        RHO3 = 2.639D0  
1449
        VEL3 = 0.714D0  
1450
        U3   = P3/(G - 1.D0)/RHO3  
1451
        CS3  = DSQRT(G*P3*(G - 1.D0)/(RHO3*(G - 1.D0) + G*P3))  
1452
 
1453
        P4   = P3  
1454
        RHO4 = 5.071D0  
1455
        VEL4 = VEL3  
1456
        U4   = P4/(G - 1.D0)/RHO4  
1457
        CS4  = DSQRT(G*P4*(G - 1.D0)/(RHO4*(G - 1.D0)+G*P4))  
1458
 
1459
        P5   = 0.666666666666D-6  
1460
        RHO5 = 1.D0  
1461
        VEL5 = 0.D0  
1462
        U5   = P5/(G - 1.D0)/RHO5  
1463
        CS5  = DSQRT(G*P5*(G - 1.D0)/(RHO5*(G - 1.D0) + G*P5))  
1464
 
1465
        X1 = 0.5D0 + (VEL1 - CS1)*TIME/(1.D0 - VEL1*CS1)  
1466
 
1467
        X2 = 0.5D0 + (VEL3 - CS3)*TIME/(1.D0 - VEL3*CS3)  
1468
 
1469
        X3 = 0.5D0 + VEL3*TIME  
1470
 
1471
        X4 = 0.5D0 + VEL4*TIME/(1.D0 - RHO5*DSQRT(1.D0 - VEL4**2)/RHO4)  
1472

        
1473
      ELSE IF (BASENM.EQ.'RBWI') THEN  
1474

        
1475
         U1   = 2.5D3  
1476
         RHO1 = 1.D0  
1477
         VEL1 = 0.D0  
1478
         P1   = (G - 1.D0)*U1*RHO1  
1479
         CS1  = DSQRT(G*P1*(G - 1.D0)/(RHO1*(G - 1.D0) + G*P1))  
1480
           
1481
         RHO3 = 0.0491D0  
1482
         VEL3 = 0.957D0  
1483
         P3   = 14.71D0  
1484
         U3   = P3/(G - 1.D0)/RHO3  
1485
         CS3  = DSQRT(G*P3*(G - 1.D0)/(RHO3*(G - 1.D0) + G*P3))  
1486
            
1487
         RHO4 = 14.39D0
1488
         VEL4 = VEL3  
1489
         P4   = P3  
1490
         U4   = P4/(G - 1.D0)/RHO4  
1491
         CS4  = DSQRT(G*P4*(G - 1.D0)/(RHO4*(G - 1.D0) + G*P4))  
1492
           
1493
         IF (TIME.LT.0.4200) THEN  
1494

        
1495
            U5   = 2.5D-2  
1496
            RHO5 = 1.D0  
1497
            VEL5 = 0.D0  
1498
            P5   = (G - 1.D0)*RHO5*U5  
1499
            CS5  = DSQRT(G*P5*(G - 1.D0)/(RHO5*(G - 1.D0) + G*P5))  
1500
              
1501
            U6   = U5  
1502
            RHO6 = RHO5  
1503
            VEL6 = VEL5  
1504
            P6   = P5  
1505
            CS6  = CS5  
1506

        
1507
         ELSE  
1508

        
1509
            RHO5 = 104.41D0  
1510
            VEL5 = 0.456D0  
1511
            P5   = 369.84D0  
1512
            U5   = P5/(G - 1.D0)/RHO5  
1513
            CS5  = DSQRT(G*P5*(G - 1.D0)/(RHO5*(G - 1.D0) + G*P5))  
1514
              
1515
            RHO6 = 117.25D0  
1516
            VEL6 = VEL5  
1517
            P6   = P5  
1518
            U6   = P6/(G - 1.D0)/RHO6  
1519
            CS6  = DSQRT(G*P6*(G - 1.D0)/(RHO6*(G - 1.D0) + G*P6))  
1520

        
1521
         END IF  
1522

        
1523
         RHO7 = 9.72D0
1524
         VEL7 = -0.882D0  
1525
         P7   = 4.639D0  
1526
         U7   = P7/(G - 1.D0)/RHO7  
1527
         CS7  = DSQRT(G*P7*(G - 1.D0)/(RHO7*(G - 1.D0) + G*P7))  
1528
           
1529
         RHO8 = 0.112D0  
1530
         VEL8 = VEL7  
1531
         P8   = P7  
1532
         U8   = P8/(G - 1.D0)/RHO8  
1533
         CS8  = DSQRT(G*P8*(G - 1.D0)/(RHO8*(G - 1.D0) + G*P8))  
1534
           
1535
         U10  = 2.5D2  
1536
         RHO10= 1.D0  
1537
         VEL10= 0.D0  
1538
         P10  = (G - 1.D0)*U10*RHO10  
1539
         CS10 = DSQRT(G*P10*(G - 1.D0)/(RHO10*(G - 1.D0) + G*P10))  
1540

        
1541
         X1 = 0.1D0 - TIME*CS1  
1542
           
1543
         X2 = 0.1D0 + TIME*(VEL3 - CS3)/(1.D0 - VEL3*CS3)  
1544
           
1545
         X3 = 0.1D0 + TIME*VEL3  
1546
           
1547
         IF (TIME.LT.0.4200D0) THEN  
1548
            X4 = 0.1D0 + TIME*0.9776D0  
1549
              
1550
            X5 = 0.9D0 - TIME*0.9274D0  
1551
              
1552
            X6 = X5  
1553
         ELSE  
1554
            X4 = 0.5106D0 + (TIME - 0.4200D0)*0.088D0  
1555
              
1556
            X5 = 0.5106D0 + (TIME - 0.4200D0)*VEL5  
1557
              
1558
            X6 = 0.5106D0 + (TIME - 0.4200D0)*0.703D0  
1559
         END IF   
1560
            
1561
         X7 = 0.9D0 + TIME*VEL7  
1562
           
1563
         X8 = 0.9D0 + TIME*(VEL8 + CS8)/(1.D0 + VEL8*CS8)  
1564
           
1565
         X9 = 0.9D0 + TIME*CS10  
1566

        
1567
      END IF  
1568
 
1569
      IF (BASENM.EQ.'RST_'.OR.BASENM.EQ.'RBW_'.OR.  
1570
     &    BASENM.EQ.'SCHN') THEN  
1571

        
1572
        DO 70 I=1,NX  
1573

        
1574
          IF (X(I).LT.X1) THEN  
1575

        
1576
            PA(I)   = P1  
1577
            RHOA(I) = RHO1  
1578
            VELA(I) = VEL1  
1579
            UA(I)   = U1  
1580
      
1581
          ELSE IF (X(I).LT.X2) THEN  
1582
         
1583
            A = (X(I) - 0.5D0)/TIME  
1584
            B = DSQRT(G - 1.D0)  
1585
            C = (B + CS1)/(B - CS1)  
1586
            D = -B/2.D0  
1587
            K = (1.D0 + A)/(1.D0 - A)  
1588
            L = C*K**D  
1589
            OCS2 = CS1  
1590
 50         FCS2 = L*(1.D0 + OCS2)**D*(OCS2 - B) +
1591
     &             (1.D0 - OCS2)**D*(OCS2+B)  
1592
            DFDCS2 = L*(1.D0 + OCS2)**D*(1.D0 + D*(OCS2 - B)/
1593
     &               (1.D0 + OCS2)) +  
1594
     &                 (1.D0 - OCS2)**D*(1.D0 - D*(OCS2 + B)/
1595
     &               (1.D0 - OCS2))  
1596
            CS2 = OCS2 - FCS2/DFDCS2  
1597
            
1598
            IF (DABS(CS2 - OCS2)/OCS2.GT.5.D-10) THEN  
1599
               OCS2 = CS2  
1600
               GOTO 50  
1601
            END IF  
1602
         
1603
            VELA(I) = (A + CS2)/(1.D0 + A*CS2)  
1604
            RHOA(I) = RHO1*((CS2**2*(G - 1.D0 - CS1**2))/  
1605
     &                (CS1**2*(G - 1.D0 - CS2**2.)))**(1.D0/(G - 1.D0))  
1606
            PA(I)   = CS2**2*(G - 1.D0)*RHOA(I)/(G - 1.D0 - CS2**2)/G  
1607
            UA(I)   = PA(I)/(G - 1.D0)/RHOA(I)
1608
  
1609
          ELSE IF (X(I).LT.X3) THEN  
1610

        
1611
            PA(I)   = P3  
1612
            RHOA(I) = RHO3  
1613
            VELA(I) = VEL3  
1614
            UA(I)   = U3  
1615

        
1616
          ELSE IF (X(I).LT.X4) THEN  
1617
  
1618
            PA(I)   = P4  
1619
            RHOA(I) = RHO4  
1620
            VELA(I) = VEL4  
1621
            UA(I)   = U4  
1622
        
1623
          ELSE
1624
  
1625
            PA(I)   = P5  
1626
            RHOA(I) = RHO5  
1627
            VELA(I) = VEL5  
1628
            UA(I)   = U5  
1629
          END IF  
1630
 
1631
70        CONTINUE  
1632
 
1633
      END IF  
1634

        
1635
      IF (BASENM.EQ.'RBWI') THEN  
1636

        
1637
        DO 80 I=1,NX  
1638
      
1639
          IF (X(I).LT.X1) THEN  
1640

        
1641
            RHOA(I) = RHO1  
1642
            VELA(I) = VEL1  
1643
            PA(I)   = P1  
1644
      
1645
          ELSE IF (X(I).LT.X2) THEN  
1646
 
1647
            A = (X(I) - 0.1D0)/TIME  
1648
            B = DSQRT(G - 1.D0)  
1649
            C = (B + CS1)/(B - CS1)  
1650
            D = -B/2.D0  
1651
            K = (1.D0 + A)/(1.D0 - A)  
1652
            L = C*K**D  
1653
            OCS2 = CS1  
1654
 52         FCS2 = L*(1.D0 + OCS2)**D*(OCS2 - B) +
1655
     &               (1.D0 - OCS2)**D*(OCS2 + B)  
1656
            DFDCS2 = L*(1.D0 + OCS2)**D*(1.D0 + D*(OCS2 - B)/
1657
     &               (1.D0 + OCS2)) +  
1658
     &                 (1.D0 - OCS2)**D*(1.D0 - D*(OCS2 + B)/
1659
     &               (1.D0 - OCS2))  
1660
            CS2 = OCS2 - FCS2/DFDCS2  
1661
         
1662
            IF (ABS(CS2-OCS2)/OCS2.GT.5.E-10)THEN  
1663
               OCS2 = CS2  
1664
               GOTO 52  
1665
            END IF 
1666
 
1667
            VELA(I) = (A+CS2)/(1.D0+A*CS2)  
1668
            RHOA(I) = RHO1*((CS2**2.*(G-1.D0-CS1**2))/  
1669
     &                (CS1**2*(G-1.D0-CS2**2)))**(1.D0/(G-1.D0))  
1670
            PA(I)   = CS2**2*(G-1.D0)*RHOA(I)/(G-1.D0-CS2**2)/G  
1671
            UA(I)   = PA(I)/(G-1.D0)/RHOA(I)  
1672

        
1673
          ELSE IF (X(I).LT.X3) THEN  
1674

        
1675
            RHOA(I) = RHO3  
1676
            VELA(I) = VEL3  
1677
            PA(I)   = P3  
1678

        
1679
          ELSE IF (X(I).LT.X4) THEN  
1680

        
1681
            RHOA(I) = RHO4  
1682
            VELA(I) = VEL4  
1683
            PA(I)   = P4  
1684

        
1685
          ELSE IF (X(I).LT.X5) THEN  
1686

        
1687
            RHOA(I) = RHO5  
1688
            VELA(I) = VEL5  
1689
            PA(I)   = P5  
1690

        
1691
          ELSE IF (X(I).LT.X6) THEN  
1692

        
1693
            RHOA(I) = RHO6  
1694
            VELA(I) = VEL6  
1695
            PA(I)   = P6  
1696

        
1697
          ELSE IF (X(I).LT.X7) THEN  
1698

        
1699
            RHOA(I) = RHO7  
1700
            VELA(I) = VEL7  
1701
            PA(I)   = P7  
1702

        
1703
          ELSE IF (X(I).LT.X8) THEN  
1704

        
1705
            RHOA(I) = RHO8  
1706
            VELA(I) = VEL8  
1707
            PA(I)   = P8  
1708
          
1709
          ELSE IF (X(I).LT.X9) THEN  
1710
         
1711
            A = (X(I) - 0.9D0)/TIME  
1712
            B = SQRT(G - 1.D0)  
1713
            C = (B + CS10)/(B - CS10)  
1714
            D = B/2.D0  
1715
            K = (1.D0 + A)/(1.D0 - A)  
1716
            L = C*K**D  
1717
            OCS2 = CS10  
1718
 54         FCS2 = L*(1.D0 - OCS2)**D*(OCS2 - B) +
1719
     &               (1.D0 + OCS2)**D*(OCS2 + B)  
1720
            DFDCS2 = L*(1.D0 - OCS2)**D*(1.D0 + D*(OCS2 - B)/
1721
     &               (1.D0 - OCS2))+  
1722
     &                 (1.D0 + OCS2)**D*(1.D0 - D*(OCS2 + B)/
1723
     &               (1.D0 + OCS2))  
1724
            CS9 = OCS2 - FCS2/DFDCS2  
1725
         
1726
            IF (DABS(CS9-OCS2)/OCS2.GT.5.D-10)THEN  
1727
               OCS2 = CS9  
1728
               GOTO 54  
1729
            END IF
1730
  
1731
            VELA(I) = (A - CS9)/(1.D0 - A*CS9)  
1732
            RHOA(I) = RHO10*((CS9 **2*(G - 1.D0 - CS10**2))/  
1733
     &               (CS10**2*(G - 1.D0 - CS9 **2)))**(1.D0/(G - 1.D0))  
1734
            PA(I)   = CS9**2*(G - 1.D0)*RHOA(I)/(G - 1.D0 - CS9**2)/G  
1735
            UA(I)   = PA(I)/(G - 1.D0)/RHOA(I)  
1736
 
1737
          ELSE   
1738
    
1739
            RHOA(I) = RHO10  
1740
            VELA(I) = VEL10  
1741
            PA(I)   = P10  
1742

        
1743
         END IF  
1744

        
1745
 80      CONTINUE
1746
  
1747
      END IF  
1748

        
1749
C     ----
1750
C     OUTPUT 
1751
C     ----
1752

        
1753
      OPEN(10,FILE='DATA/'//OUTFIL,FORM='FORMATTED',STATUS='NEW')  
1754
      WRITE(10,111) NSTEP,TIME
1755
 111  FORMAT('N =', I6, 3X, 'TIME = ', 1PE10.3)   
1756

        
1757
      DO 85 I=1,NX  
1758
        WRITE(10,200) X(I),P(I),PA(I),RHO(I),RHOA(I),  
1759
     &                VEL(I),VELA(I)  
1760
 85     CONTINUE  
1761

        
1762
 200  FORMAT(F6.4,1X,2(F11.4,1X),2(F9.4,1X),2(F7.5,1X))  
1763
      CLOSE(10)  
1764

        
1765
      NOUT1 = 0  
1766
      TOUT1 = 0.D0  
1767

        
1768
      CALL FILNAM  
1769

        
1770
      RETURN  
1771
      END  
1772

        
1773
C     -------- 
1774
CN    NAME: E O S 
1775
C     --------
1776

        
1777
CP    PURPOSE: 
1778
CP    COMPUTES DERIVED THERMODYNAMICAL QUANTITIES
1779
C
1780

        
1781
CC    COMMENTS: 
1782
CC    GAMMA-LAW EOS
1783

        
1784
      SUBROUTINE EOS(N,RHO,U,G,P,H,CS,DPDRH,DPDU)
1785

        
1786
      IMPLICIT NONE
1787

        
1788
      INCLUDE 'size'
1789

        
1790
C     -------
1791
C     ARGUMENTS
1792
C     -------
1793

        
1794
      INTEGER         N
1795

        
1796
      DOUBLEPRECISION U(-4:MN5),RHO(-4:MN5),P(-4:MN5),H(-4:MN5),
1797
     &                CS(-4:MN5),DPDRH(-4:MN5),DPDU(-4:MN5)
1798
      
1799
      DOUBLEPRECISION G
1800

        
1801
C     -----------
1802
C     INTERNAL VARIABLES
1803
C     -----------
1804

        
1805
      INTEGER         I
1806

        
1807
      DOUBLEPRECISION GAM1
1808

        
1809
      DO 10 I=1,N
1810
         GAM1     = G -1.D0
1811
         P(I)     = GAM1*RHO(I)*U(I)
1812
         DPDRH(I) = GAM1*U(I)
1813
         DPDU(I)  = GAM1*RHO(I)
1814
         H(I)     = 1.D0 + U(I) + P(I)/RHO(I)
1815
         CS(I)    = DSQRT((DPDRH(I) + P(I)*DPDU(I)/RHO(I)/RHO(I))/H(I))
1816
 10      CONTINUE
1817

        
1818
      RETURN
1819
      END
1820

        
1821
C     -------- 
1822
CN    NAME: C O E F
1823
C     --------
1824

        
1825
CP    PURPOSE: 
1826
CP    COMPUTES THE COEFFICIENTS FOR INTERPOLATED VALUES 
1827
C
1828

        
1829
CC    COMMENTS: 
1830
CC    COMPUTES DE GRID-DEPENDENT COEFFICIENTS APPEARING IN EQS.59-61 OF MARTI
1831
CC    AND MUELLER (1996), JCP, VOL. 123, 1-14
1832

        
1833
      SUBROUTINE COEF(N,COEFF1,COEFF2,COEFF3,COEFF4,COEFF5)
1834
 
1835
      IMPLICIT NONE
1836
 
1837
      INCLUDE 'size'
1838
 
1839
C     --------
1840
C     ARGUMENTS
1841
C     --------
1842

        
1843
      INTEGER N
1844
 
1845
      DOUBLEPRECISION COEFF1(-4:MN5),COEFF2(-4:MN5),COEFF3(-4:MN5)
1846
      DOUBLEPRECISION COEFF4(-4:MN5),COEFF5(-4:MN5)
1847

        
1848
C     ----------
1849
C     COMMON BLOCKS
1850
C     ----------
1851

        
1852
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
1853
      COMMON /GRD/    X,XL,XR,DX
1854

        
1855
C     --------------
1856
C     INTERNAL VARIABLES
1857
C     --------------
1858
      
1859
      INTEGER         I
1860

        
1861
      DOUBLEPRECISION SCRCH1(-4:MN6),SCRCH2(-4:MN6),SCRCH3(-4:MN6),
1862
     &                SCRCH4(-4:MN6)
1863
 
1864
 
1865
      DO 10 I=0,N+2
1866
 
1867
        SCRCH1(I) = DX(I) + DX(I-1)
1868
        SCRCH2(I) = SCRCH1(I) + DX(I)
1869
        SCRCH3(I) = SCRCH1(I) + DX(I-1)
1870

        
1871
10      CONTINUE
1872
 
1873
      DO 20 I=0,N+1
1874

        
1875
        SCRCH4(I) = DX(I)/(SCRCH1(I) + DX(I+1))
1876
        COEFF1(I) = SCRCH4(I)*SCRCH3(I)/SCRCH1(I+1)
1877
        COEFF2(I) = SCRCH4(I)*SCRCH2(I+1)/SCRCH1(I)
1878

        
1879
20      CONTINUE
1880
 
1881
      DO 30 I=0,N
1882
        
1883
        SCRCH4(I) = 1.D0/(SCRCH1(I) + SCRCH1(I+2))
1884
        COEFF3(I) = -SCRCH4(I)*DX(I)*SCRCH1(I)/SCRCH3(I+1)
1885
        COEFF4(I) = SCRCH4(I)*DX(I+1)*SCRCH1(I+2)/SCRCH2(I+1)
1886
        COEFF5(I) = DX(I) - 2.D0*(DX(I+1)*COEFF3(I) + DX(I)*COEFF4(I))
1887
        COEFF5(I) = COEFF5(I)/SCRCH1(I+1)
1888

        
1889
30      CONTINUE
1890
  
1891
      RETURN
1892
      END
1893
 
1894
C     -------- 
1895
CN    NAME: I N T E R P
1896
C     --------
1897

        
1898
CP    PURPOSE: 
1899
CP    COMPUTES INTERPOLATED VALUES AT INTERFACES
1900
C
1901

        
1902
CC    COMMENTS: 
1903
CC    STEP 1 IN THE RECONSTRUCTION PROCEDURE (SEE APPENDIX I IN MARTI 
1904
CC    & MUELLER 1996)
1905

        
1906
      SUBROUTINE INTERP(N,UM,U,UP,DELU)
1907
 
1908
      IMPLICIT NONE
1909
 
1910
      INCLUDE 'size'
1911

        
1912
C     --------
1913
C     ARGUMENTS
1914
C     --------
1915

        
1916
      INTEGER N
1917
 
1918
      DOUBLEPRECISION UM(-4:MN5),U(-4:MN5),UP(-4:MN5),DELU(-4:MN5)
1919

        
1920
C     ---------
1921
C     COMMON BLOCKS
1922
C     ---------
1923
 
1924
      DOUBLEPRECISION COEFF1(-4:MN5),COEFF2(-4:MN5),COEFF3(-4:MN5)
1925
      DOUBLEPRECISION COEFF4(-4:MN5),COEFF5(-4:MN5)
1926
      COMMON /COEFF/  COEFF1,COEFF2,COEFF3,COEFF4,COEFF5
1927
 
1928
C     -----------
1929
C     INTERNAL VARIABLES
1930
C     -----------
1931

        
1932
      INTEGER         I
1933

        
1934
      DOUBLEPRECISION SCRCH1(-4:MN6),SCRCH2(-4:MN6)
1935
 
1936
      DOUBLEPRECISION SDELU
1937
 
1938
      DO 10 I=-2,N+3
1939

        
1940
        SCRCH1(I) = U(I) - U(I-1)
1941

        
1942
10    CONTINUE
1943
 
1944
C     ------------------------------------------------------
1945
C     DELU(I) AS IN EQ.61 OF MARTI AND MUELLER (1996), JCP, VOL. 123, 1-14
1946
C     ------------------------------------------------------
1947

        
1948
      DO 20 I=0,N+1
1949

        
1950
        DELU(I) = COEFF1(I)*SCRCH1(I+1) + COEFF2(I)*SCRCH1(I)
1951

        
1952
20    CONTINUE
1953
 
1954
C     -------------------------------------------------------
1955
C     DELU(I) AS IN EQ.60 OF MARTI AND MUELLER (1996), JCP, VOL. 123, 1-14
1956
C     -------------------------------------------------------
1957

        
1958
      DO 30 I=0,N+1
1959

        
1960
        IF (SCRCH1(I+1)*SCRCH1(I).GT.0.D0) THEN
1961
          SDELU     = DELU(I)/DABS(DELU(I))
1962
          SCRCH2(I) = MIN(DABS(SCRCH1(I)),DABS(SCRCH1(I+1)))
1963
          DELU(I)   = MIN(DABS(DELU(I)),2.D0*SCRCH2(I))*SDELU
1964
        ELSE
1965
          DELU(I)   = 0.D0
1966
        END IF
1967

        
1968
30    CONTINUE
1969

        
1970
C     ----------------------------------------------
1971
C     INTERFACE VALUES AS IN  EQ.59 OF MARTI AND MUELLER (1996), JCP, 
1972
C     VOL. 123, 1-14
1973
C     ----------------------------------------------
1974

        
1975
      DO 40 I=0,N
1976

        
1977
        UP(I)   = U(I)  + COEFF5(I)*SCRCH1(I+1) + COEFF3(I)*DELU(I+1)
1978
        UP(I)   = UP(I) + COEFF4(I)*DELU(I)
1979
        UM(I+1) = UP(I)
1980

        
1981
40    CONTINUE
1982
 
1983
      RETURN
1984
      END
1985
 
1986
C     -------- 
1987
CN    NAME: D E T E C T
1988
C     --------
1989

        
1990
CP    PURPOSE: 
1991
CP    DETECTS CONTACT DISCONTINUITIES AND STEEPENS THE CORRESPONDING
1992
CP    RECONSTRUCTED VALUES AT INTERFACES
1993
C
1994

        
1995
CC    COMMENTS: 
1996
CC    STEP 2 IN THE RECONSTRUCTION PROCEDURE (SEE APPENDIX I IN MARTI 
1997
CC    & MUELLER 1996)
1998

        
1999
      SUBROUTINE DETECT(N,UM,U,UP,DELU)
2000
 
2001
      IMPLICIT NONE
2002
 
2003
      INCLUDE 'size'
2004
 
2005
C     -------
2006
C     ARGUMENTS
2007
C     -------
2008

        
2009
      INTEGER         N
2010
 
2011
      DOUBLEPRECISION UM(-4:MN5),U(-4:MN5),UP(-4:MN5),DELU(-4:MN5)
2012

        
2013
C     -----------
2014
C     COMMON BLOCKS
2015
C     -----------
2016
 
2017
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
2018
      COMMON /GRD/    X,XL,XR,DX
2019

        
2020
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),W(-4:MN5),
2021
     &                UU(-4:MN5),CS(-4:MN5),H(-4:MN5),DPDRH(-4:MN5),
2022
     &                DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),E(-4:MN5)
2023
      COMMON /HYDRO/  P,RHO,VEL,W,UU,CS,H,DPDRH,DPDU,R,M,E
2024
 
2025
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
2026
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
2027
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
2028
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
2029
 
2030
      DOUBLEPRECISION GB
2031
      COMMON /ADIND/  GB
2032
 
2033
C     ---------
2034
C     INTERNAL VARIABLES
2035
C     ---------
2036

        
2037
      INTEGER         I
2038

        
2039
      DOUBLEPRECISION D2U(-4:MN6),ETA(-4:MN6),ETATIL(-4:MN6)
2040
      
2041
      DOUBLEPRECISION SCRCH1(-4:MN6),SCRCH2(-4:MN6),SCRCH3(-4:MN6),
2042
     &                SCRCH4(-4:MN6)
2043

        
2044
 
2045
      DO 10 I=-1,N+3
2046

        
2047
        SCRCH1(I) = DX(I) + DX(I-1)
2048
        SCRCH2(I) = SCRCH1(I) + DX(I+1)
2049
        SCRCH3(I) = U(I) - U(I-1)
2050
        SCRCH1(I) = SCRCH3(I)/SCRCH1(I)
2051

        
2052
10      CONTINUE
2053
 
2054
      DO 20 I=-1,N+2
2055
        
2056
        D2U(I)    = (SCRCH1(I+1) - SCRCH1(I))/SCRCH2(I)
2057
        SCRCH4(I) = X(I) - X(I-1)
2058
        SCRCH4(I) = SCRCH4(I)*SCRCH4(I)*SCRCH4(I)
2059

        
2060
20      CONTINUE
2061
 
2062
      DO 30 I=0,N+1
2063

        
2064
        SCRCH1(I) = D2U(I+1)*D2U(I-1)
2065
        SCRCH3(I) = DABS(U(I+1) - U(I-1))
2066
        SCRCH3(I) = SCRCH3(I) - EPSLN*MIN(DABS(U(I+1)),DABS(U(I-1)))
2067

        
2068
30      CONTINUE
2069
 
2070
C     ------------------------------------------------------
2071
C     ETATIL(I) AS IN EQ.67 OF MARTI AND MUELLER (1996), JCP, VOL. 123, 1-14
2072
C     -------------------------------------------------------
2073

        
2074
      DO 40 I=0,N+1
2075
      
2076
        IF ((U(I+1) - U(I-1)).EQ.0.D0) THEN
2077
           SCRCH2(I) = SMALL*SMLRHO
2078
        ELSE
2079
           SCRCH2(I) = U(I+1) - U(I-1)
2080
        END IF
2081

        
2082
        IF ((SCRCH1(I).GT.0.D0).OR.(SCRCH3(I).LT.0.D0)) THEN
2083
           ETATIL(I) = 0.D0
2084
        ELSE
2085
           ETATIL(I) = (D2U(I-1) - D2U(I+1))*(SCRCH4(I) + SCRCH4(I+1))
2086
           ETATIL(I) = ETATIL(I)/(X(I+1) - X(I-1))/SCRCH2(I)
2087
        END IF
2088

        
2089
40      CONTINUE
2090
 
2091
C     -------------------------------------------------------
2092
C     ETA(I) AS IN EQ.66 OF MARTI AND MUELLER (1996), JCP, VOL. 123, 1-14
2093
C     ONLY FOR ZONES VERIFYING EQ.63 (OTHERWISE, ZERO)
2094
C     -------------------------------------------------------
2095

        
2096
      DO 50 I=0,N+1
2097

        
2098
        ETA(I)    = MAX(0.D0,MIN(ETA1*(ETATIL(I) - ETA2),1.D0))
2099
        SCRCH1(I) = DABS(P  (I+1) - P  (I-1))/MIN(P  (I+1),P  (I-1))
2100
        SCRCH2(I) = DABS(RHO(I+1) - RHO(I-1))/MIN(RHO(I+1),RHO(I-1))
2101

        
2102
50      CONTINUE
2103
 
2104
      DO 60 I=0,N+1
2105

        
2106
        IF (GB*AK0*SCRCH2(I).LT.SCRCH1(I)) THEN
2107
          ETA(I) = 0.D0
2108
        END IF
2109

        
2110
60      CONTINUE
2111
 
2112
C     ----------------------------
2113
C     NEW RECONSTRUCTED VALUES (EQ.65)
2114
C     ----------------------------
2115

        
2116
      DO 70 I=0,N+1
2117

        
2118
        SCRCH1(I) = U(I-1) + 0.5D0*DELU(I-1)
2119
        SCRCH2(I) = U(I+1) - 0.5D0*DELU(I+1)
2120
        UM(I)     = UM(I)*(1.D0-ETA(I)) + SCRCH1(I)*ETA(I)
2121
        UP(I)     = UP(I)*(1.D0-ETA(I)) + SCRCH2(I)*ETA(I)
2122

        
2123
70    CONTINUE
2124
 
2125
      RETURN
2126
      END
2127
 
2128
C     -------- 
2129
CN    NAME: M O N O T 
2130
C     --------
2131

        
2132
CP    PURPOSE: 
2133
CP    RECONSTRUCTED VALUES OBTAINED IN STEPS 1 TO 3 ARE MODIFIED SUCH THAT 
2134
CP    THE INTERPOLATION PARABOLA IN EACH ZONE IS MONOTONE
2135
C
2136

        
2137
CC    COMMENTS: 
2138
CC    STEP 4 IN THE RECONSTRUCTION PROCEDURE (SEE APPENDIX I IN MARTI 
2139
CC    & MUELLER 1996)
2140

        
2141
      SUBROUTINE MONOT(N,UM,U,UP,DU,U6)
2142
 
2143
      IMPLICIT NONE
2144
 
2145
      INCLUDE 'size'
2146

        
2147
C     --------
2148
C     ARGUMENTS
2149
C     --------
2150

        
2151
      INTEGER         N
2152
 
2153
      DOUBLEPRECISION UM(-4:MN5),U(-4:MN5),UP(-4:MN5),DU(-4:MN5)
2154
      DOUBLEPRECISION U6(-4:MN5)
2155
 
2156
C     --------
2157
C     INTERNAL VARIABLES
2158
C     --------
2159

        
2160
      INTEGER         I
2161

        
2162
      DOUBLEPRECISION SCRCH1(-4:MN6),SCRCH2(-4:MN6),SCRCH3(-4:MN6)
2163
  
2164
C     -----------------------------------------------------
2165
C     NEW RECONSRUCTED VALUES IF CONDITION IN EQ.73 OF MARTI & MUELLER
2166
C     1996 HOLDS
2167
C     -----------------------------------------------------
2168

        
2169
      DO 10 I=0,N+1
2170

        
2171
        DU(I)     = UP(I) - UM(I)
2172
        SCRCH1(I) = UP(I) - U(I)
2173
        SCRCH1(I) = SCRCH1(I)*(UM(I)-U(I))
2174

        
2175
10      CONTINUE
2176
 
2177
      DO 20 I=0,N+1
2178

        
2179
        IF (SCRCH1(I).GE.0.D0) THEN
2180
          UM(I) = U(I)
2181
          UP(I) = U(I)
2182
        END IF
2183

        
2184
20      CONTINUE
2185

        
2186
C     --------------------------------------------------------
2187
C     NEW RECONSTRUCTED VALUES IF CONDITION IN EQ.74 OR EQ.75 OF MARTI 
2188
C     & MUELLER 1996 HOLDS
2189
C     --------------------------------------------------------
2190
 
2191
      DO 30 I=0,N+1
2192

        
2193
        DU(I)     = UP(I) - UM(I)
2194
        SCRCH1(I) = (UP(I) - U(I))*(UM(I) - U(I))
2195

        
2196
        IF (SCRCH1(I).EQ.0.D0) THEN
2197
          SCRCH2(I) = UM(I)
2198
          SCRCH3(I) = UP(I)
2199
        ELSE
2200
          SCRCH2(I) = 3.D0*U(I) - 2.D0*UP(I)
2201
          SCRCH3(I) = 3.D0*U(I) - 2.D0*UM(I)
2202
        END IF
2203

        
2204
30      CONTINUE
2205
 
2206
      DO 40 I=0,N+1
2207

        
2208
        IF (DU(I)*(UM(I) - SCRCH2(I)).LT.0.D0) THEN
2209
          UM(I) = SCRCH2(I)
2210
        END IF
2211

        
2212
        IF (DU(I)*(SCRCH3(I) - UP(I)).LT.0.D0) THEN
2213
          UP(I) = SCRCH3(I)
2214
        END IF
2215

        
2216
40      CONTINUE
2217
 
2218
C     -------------------------------------------
2219
C     COMPUTATION OF AUXILIAR VARIABLES FOR TIME ADVANCE
2220
C     -------------------------------------------
2221
 
2222
      DO 50 I=0,N+1
2223

        
2224
        DU(I) = UP(I) - UM(I)
2225
        U6(I) = 6.D0*U(I) - 3.D0*(UM(I) + UP(I))
2226

        
2227
50      CONTINUE
2228
 
2229
      RETURN
2230
      END
2231

        
2232
C     -------- 
2233
CN    NAME: F L A T E N  
2234
C     --------
2235

        
2236
CP    PURPOSE: 
2237
CP    THIS SUBROUTINE FLATTENS ZONE STRUCTURE IN REGIONS WHERE SHOCKS
2238
CP    ARE TOO THIN    
2239
C
2240

        
2241
CC    COMMENTS: 
2242
CC    STEP 3 IN THE RECONSTRUCTION PROCEDURE (SEE APPENDIX I IN MARTI 
2243
CC    & MUELLER 1996)
2244

        
2245
      SUBROUTINE FLATEN
2246

        
2247
      IMPLICIT NONE
2248
 
2249
      INCLUDE 'size'
2250

        
2251
C     -------------
2252
C     COMMON BLOCKS
2253
C     -------------
2254

        
2255
      INTEGER         NEND,NOUT,ITSTP,NX
2256
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
2257

        
2258
      DOUBLEPRECISION FLATN(-4:MN5),FLATN1(-4:MN5)
2259
      COMMON /FLAT/   FLATN,FLATN1
2260

        
2261
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),W(-4:MN5),
2262
     &                U(-4:MN5),CS(-4:MN5),H(-4:MN5),DPDRH(-4:MN5),
2263
     &                DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),E(-4:MN5)
2264
      COMMON /HYDRO/  P,RHO,VEL,W,U,CS,H,DPDRH,DPDU,R,M,E
2265

        
2266
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
2267
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
2268
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
2269
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
2270
 
2271
C     ------------
2272
C     INTERNAL VARIABLES
2273
C     ------------
2274

        
2275
      INTEGER         I
2276

        
2277
      DOUBLEPRECISION SCRCH1(-4:MN6),SCRCH2(-4:MN6),SCRCH3(-4:MN6)
2278
 
2279
      DOUBLEPRECISION DP(-4:MN5),DVEL(-4:MN5),DP2
2280
 
2281
      DO 20 I=-2,NX+3
2282

        
2283
        DP(I)     = P(I+1)   - P(I-1)
2284
        DVEL(I)   = VEL(I+1) - VEL(I-1)
2285
        SCRCH1(I) = EPSILN*DMIN1(P(I+1),P(I-1)) - DABS(DP(I))
2286

        
2287
        IF (SCRCH1(I).LT.0.D0.AND.DVEL(I).LT.0.D0) THEN
2288
          SCRCH1(I) = 1.D0
2289
        ELSE
2290
          SCRCH1(I) = 0.D0
2291
        END IF
2292

        
2293
20      CONTINUE
2294
 
2295
      DO 30 I=-1,NX+2
2296

        
2297
        DP2 = P(I+2) - P(I-2)
2298

        
2299
          IF (DP2.EQ.0.D0) THEN
2300

        
2301
            IF (DP(I).EQ.0.D0) THEN
2302
              SCRCH2(I) = -OMG1
2303
            ELSE
2304
              SCRCH2(I) = 1.D0 - OMG1
2305
            END IF
2306

        
2307
          ELSE
2308
            SCRCH2(I) = DP(I)/DP2 - OMG1
2309
        END IF
2310

        
2311
        SCRCH3(I) = SCRCH1(I)*DMAX1(0.D0,SCRCH2(I)*OMG2)
2312

        
2313
30    CONTINUE
2314
 
2315
      DO 40 I=0,NX+1
2316

        
2317
        IF (DP(I).LT.0.D0) THEN
2318
          SCRCH2(I) = SCRCH3(I+1)
2319
        ELSE
2320
          SCRCH2(I) = SCRCH3(I-1)
2321
        END IF
2322

        
2323
40    CONTINUE
2324
 
2325
      DO 45 I=0,NX+1
2326

        
2327
        FLATN(I)  = DMAX1(SCRCH2(I),SCRCH3(I))
2328
        FLATN(I)  = DMAX1(0.D0,DMIN1(1.D0,FLATN(I)))
2329
        FLATN1(I) = 1.D0 - FLATN(I)
2330

        
2331
45    CONTINUE
2332
 
2333
      RETURN
2334
      END
2335

        
2336

        
2337
C     -------- 
2338
CN    NAME: S T A T 1 D 
2339
C     --------
2340

        
2341
CP    PURPOSE: 
2342
CP    THIS SUBROUTINE CALCULATES EFFECTIVE SECOND-ORDER-ACCURATE LEFT 
2343
CP    AND RIGHT STATES FOR RIEMANN PROBLEMS IN ONE DIMENSIONAL 
2344
CP    CALCULATIONS.
2345
C
2346

        
2347
CC    COMMENTS: 
2348
CC    THIS ROUTINE CLOSELY FOLLOWS THE ANALYTICAL DEVELOPMENTS DESCRIBED IN 
2349
CC    MARTI & MUELLER, JCP, 1996
2350

        
2351
      SUBROUTINE STAT1D
2352

        
2353
      IMPLICIT NONE
2354
 
2355
      INCLUDE 'size'
2356

        
2357
C     ---------
2358
C     COMMON BLOCKS
2359
C     ---------
2360

        
2361
      INTEGER         NEND,NOUT,ITSTP,NX
2362
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
2363

        
2364
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
2365
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
2366
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,SMALLU,
2367
     &                GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
2368

        
2369
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),W(-4:MN5),
2370
     &                U(-4:MN5),CS(-4:MN5),H(-4:MN5),DPDRH(-4:MN5),
2371
     &                DPDU(-4:MN5),R(-4:MN5),M(-4:MN5),E(-4:MN5)
2372
      COMMON /HYDRO/  P,RHO,VEL,W,U,CS,H,DPDRH,DPDU,R,M,E
2373

        
2374
      DOUBLEPRECISION PM(-4:MN5),PP(-4:MN5),RHOM(-4:MN5),RHOP(-4:MN5),
2375
     &                VELM(-4:MN5),VELP(-4:MN5),UM(-4:MN5),UP(-4:MN5)
2376
      COMMON /UMP/    PM,PP,RHOM,RHOP,VELM,VELP,UM,UP
2377

        
2378
      DOUBLEPRECISION DP(-4:MN5),P6(-4:MN5),DRHO(-4:MN5),RHO6(-4:MN5),
2379
     &                DVEL(-4:MN5),VEL6(-4:MN5),DU(-4:MN5),U6(-4:MN5)
2380
      COMMON /U6/     DP,P6,DRHO,RHO6,DVEL,VEL6,DU,U6
2381

        
2382
      DOUBLEPRECISION PL(-4:MN6),PR(-4:MN6),RHOL(-4:MN6),RHOR(-4:MN6),
2383
     &                VELL(-4:MN6),VELR(-4:MN6),UL(-4:MN6),UR(-4:MN6),
2384
     &                CSL(-4:MN6),CSR(-4:MN6),RL(-4:MN6),RR(-4:MN6),
2385
     &                ML(-4:MN6),MR(-4:MN6),EL(-4:MN6),ER(-4:MN6)
2386
      COMMON /INTERF/ PL,PR,RHOL,RHOR,VELL,VELR,UL,UR,CSL,CSR,RL,RR,ML,
2387
     &                MR,EL,ER
2388

        
2389
      DOUBLEPRECISION FICT(-4:MN5)
2390
      COMMON /FICT/   FICT
2391

        
2392
      DOUBLEPRECISION TIME,DT
2393
      COMMON /ZEIT/   TIME,DT
2394

        
2395
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
2396
      COMMON /GRD/    X,XL,XR,DX
2397

        
2398
      DOUBLEPRECISION GB
2399
      COMMON /ADIND/  GB
2400
 
2401
C     ---------
2402
C     INTERNAL VARIABLES
2403
C     ---------
2404

        
2405
      INTEGER I
2406
 
2407
      DOUBLEPRECISION LAMB1(-4:MN5),LAMB2(-4:MN5),LAMB3(-4:MN5)
2408
 
2409
      DOUBLEPRECISION P1L(-4:MN6),P1R(-4:MN6),P2L(-4:MN6),P2R(-4:MN6),
2410
     &                P3L(-4:MN6),P3R(-4:MN6),RHO1L(-4:MN6),
2411
     &                RHO1R(-4:MN6),RHO2L(-4:MN6),RHO2R(-4:MN6),
2412
     &                RHO3L(-4:MN6),RHO3R(-4:MN6),VEL1L(-4:MN6),
2413
     &                VEL1R(-4:MN6),VEL2L(-4:MN6),VEL2R(-4:MN6),
2414
     &                VEL3L(-4:MN6),VEL3R(-4:MN6)
2415
 
2416
      DOUBLEPRECISION BETA1L(-4:MN6),BETA2L(-4:MN6),BETA3L(-4:MN6),
2417
     &                BETA1R(-4:MN6),BETA2R(-4:MN6),BETA3R(-4:MN6)
2418

        
2419
 
2420
      DOUBLEPRECISION SCRCH1(-4:MN6),SCRCH2(-4:MN6),SCRCH3(-4:MN6),
2421
     &                SCRCH4(-4:MN6),SCRCH5(-4:MN6)
2422

        
2423
      DO 10 I=0,NX+1
2424

        
2425
        LAMB1(I) = (VEL(I) - CS(I))/(1.D0 - VEL(I)*CS(I))
2426
        LAMB2(I) =  VEL(I)
2427
        LAMB3(I) = (VEL(I) + CS(I))/(1.D0 + VEL(I)*CS(I))
2428

        
2429
10      CONTINUE
2430

        
2431
      CALL AVRG1D(NX,DT,LAMB1,PM,PP,DP,P6,RHOM,RHOP,DRHO,RHO6,
2432
     &            VELM,VELP,DVEL,VEL6,
2433
     &            P1L,P1R,RHO1L,RHO1R,VEL1L,VEL1R)
2434
 
2435
      CALL AVRG1D(NX,DT,LAMB2,PM,PP,DP,P6,RHOM,RHOP,DRHO,RHO6,
2436
     &            VELM,VELP,DVEL,VEL6,
2437
     &            P2L,P2R,RHO2L,RHO2R,VEL2L,VEL2R)
2438
 
2439
      CALL AVRG1D(NX,DT,LAMB3,PM,PP,DP,P6,RHOM,RHOP,DRHO,RHO6,
2440
     &            VELM,VELP,DVEL,VEL6,
2441
     &            P3L,P3R,RHO3L,RHO3R,VEL3L,VEL3R)
2442
 
2443
C     -------------
2444
C     EFFECTIVE LEFT STATES
2445
C     -------------
2446
 
2447
      DO 20 I=1,NX+1
2448

        
2449
        SCRCH1(I) = P3L(I)/(GB - 1.D0)/RHO3L(I)
2450
        SCRCH1(I) = 1.D0 + SCRCH1(I) + P3L(I)/RHO3L(I)
2451
        SCRCH2(I) = DSQRT(GB*P3L(I)/RHO3L(I)/SCRCH1(I))
2452
        SCRCH1(I) = SCRCH1(I)*SCRCH2(I)
2453
        SCRCH3(I) = 1.D0/(1.D0 - VEL3L(I)**2)
2454
        SCRCH3(I) = RHO3L(I)*SCRCH3(I)
2455
        SCRCH4(I) = 0.D0
2456
        SCRCH5(I) = 0.D0
2457
 
2458
        BETA1L(I) = 0.5D0*(VEL3L(I) - VEL1L(I) - (P3L(I) - P1L(I))/
2459
     &              SCRCH3(I)/SCRCH1(I) - DT*SCRCH4(I))
2460
        BETA2L(I) = RHO3L(I) - RHO2L(I) - (P3L(I) - P2L(I))/
2461
     &              SCRCH1(I)/SCRCH2(I)
2462
        BETA3L(I) = -0.5D0*DT*SCRCH5(I)
2463

        
2464
20      CONTINUE
2465
 
2466
      DO 25 I=1,NX+1
2467

        
2468
        IF (LAMB3(I-1).LE.0.D0) THEN
2469
          BETA3L(I) = 0.D0
2470
        END IF
2471

        
2472
        IF (LAMB2(I-1).LE.0.D0) THEN
2473
          BETA2L(I) = 0.D0
2474
        END IF
2475

        
2476
        IF (LAMB1(I-1).LE.0.D0) THEN
2477
          BETA1L(I) = 0.D0
2478
        END IF
2479

        
2480
25      CONTINUE
2481
 
2482
      DO 30 I=1,NX+1
2483
        
2484
        PL(I)     = P3L(I) + SCRCH3(I)*SCRCH1(I)*(BETA1L(I) + BETA3L(I))
2485
        PL(I)     = DMAX1(SMALLP,PL(I))
2486
 
2487
        RHOL(I)   = RHO3L(I) + SCRCH3(I)/SCRCH2(I)*(BETA1L(I) + 
2488
     &              BETA3L(I)) - BETA2L(I)
2489
        RHOL(I)   = DMAX1(SMLRHO,RHOL(I))
2490
 
2491
        VELL(I)   = VEL3L(I) + BETA3L(I) - BETA1L(I)
2492
 
2493
        UL(I)     = PL(I)/(GB - 1.D0)/RHOL(I)
2494
 
2495
        SCRCH1(I) = 1.D0 + UL(I) + PL(I)/RHOL(I)
2496
 
2497
        CSL(I)    = DSQRT(GB*PL(I)/RHOL(I)/SCRCH1(I))
2498
 
2499
30      CONTINUE
2500
 
2501
C     ----------
2502
C     EFFECTIVE RIGHT STATES
2503
C     ----------
2504

        
2505
      DO 35 I=1,NX+1
2506

        
2507
        SCRCH1(I) = P1R(I)/(GB - 1.D0)/RHO1R(I)
2508
        SCRCH1(I) = 1.D0 + SCRCH1(I) + P1R(I)/RHO1R(I)
2509
        SCRCH2(I) = DSQRT(GB*P1R(I)/RHO1R(I)/SCRCH1(I))
2510
        SCRCH1(I) = SCRCH1(I)*SCRCH2(I)
2511
        SCRCH3(I) = 1.D0/(1.D0 - VEL1R(I)**2)
2512
        SCRCH3(I) = RHO1R(I)*SCRCH3(I)
2513
        SCRCH4(I) = 0.D0
2514
        SCRCH5(I) = 0.D0
2515
 
2516
        BETA1R(I) = -0.5D0*DT*SCRCH4(I)
2517
        BETA2R(I) = RHO1R(I) - RHO2R(I) - (P1R(I)-P2R(I))/
2518
     &              SCRCH1(I)/SCRCH2(I)
2519
        BETA3R(I) = -0.5D0*(VEL1R(I) - VEL3R(I) + (P1R(I) - P3R(I))/
2520
     &              SCRCH3(I)/SCRCH1(I) + DT*SCRCH5(I))
2521
35      CONTINUE
2522
 
2523
      DO 40 I=1,NX+1
2524

        
2525
        IF (LAMB3(I).GE.0.D0) THEN
2526
          BETA3R(I) = 0.D0
2527
        END IF
2528

        
2529
        IF (LAMB2(I).GE.0.D0) THEN
2530
          BETA2R(I) = 0.D0
2531
        END IF
2532

        
2533
        IF (LAMB1(I).GE.0.D0) THEN
2534
          BETA1R(I) = 0.D0
2535
        END IF
2536

        
2537
40    CONTINUE
2538
 
2539
      DO 45 I=1,NX+1
2540

        
2541
        PR(I)     = P1R(I) + SCRCH3(I)*SCRCH1(I)*(BETA1R(I) + BETA3R(I))
2542
        PR(I)     = DMAX1(SMALLP,PR(I))
2543
 
2544
        RHOR(I)   = RHO1R(I) + SCRCH3(I)/SCRCH2(I)*(BETA1R(I) + 
2545
     &              BETA3R(I)) - BETA2R(I)
2546
        RHOR(I)   = DMAX1(SMLRHO,RHOR(I))
2547
 
2548
        VELR(I)   = VEL1R(I) + BETA3R(I) - BETA1R(I)
2549
 
2550
        UR(I)     = PR(I)/(GB - 1.D0)/RHOR(I)
2551
 
2552
        SCRCH1(I) = 1.D0 + UR(I) + PR(I)/RHOR(I)
2553
 
2554
        CSR(I)    = DSQRT(GB*PR(I)/RHOR(I)/SCRCH1(I))
2555
 
2556
45      CONTINUE
2557
 
2558
C     --------------
2559
C     CONSERVED VARIABLES
2560
C     --------------
2561

        
2562
      DO 60 I=1,NX+1
2563

        
2564
         SCRCH1(I) = 1.D0/DSQRT(1.D0 - VELL(I)*VELL(I))
2565
         SCRCH2(I) = 1.D0 + UL(I) + PL(I)/RHOL(I)
2566
         RL(I)     = RHOL(I)*SCRCH1(I)
2567
         ML(I)     = RL(I)*SCRCH2(I)*SCRCH1(I)*VELL(I)
2568
         EL(I)     = RL(I)*SCRCH2(I)*SCRCH1(I) - PL(I) - RL(I)
2569

        
2570
60       CONTINUE
2571
         
2572
      DO 70 I=1,NX+1
2573

        
2574
         SCRCH1(I) = 1.D0/DSQRT(1.D0 - VELR(I)*VELR(I))
2575
         SCRCH2(I) = 1.D0 + UR(I) + PR(I)/RHOR(I)
2576
         RR(I)     = RHOR(I)*SCRCH1(I)
2577
         MR(I)     = RR(I)*SCRCH2(I)*SCRCH1(I)*VELR(I)
2578
         ER(I)     = RR(I)*SCRCH2(I)*SCRCH1(I) - PR(I) - RR(I)
2579

        
2580
70       CONTINUE
2581
 
2582
      RETURN
2583
      END
2584

        
2585
C     -------- 
2586
CN    NAME: A V R G 1 D 
2587
C     --------
2588

        
2589
CP    PURPOSE: 
2590
CP    THIS SUBROUTINE CALCULATES AVERAGES OF QUANTITIES P,RHO,VEL, OVER
2591
CP    THE PART OF THE DOMAIN OF DEPENDENCE FOR THE LAMBDA 
2592
CP    CHARACTERISTIC OF RADM(I) FOR THE TIME INTERVAL (T(N),T(N+1)).
2593
C
2594

        
2595
CC    COMMENTS: 
2596
CC    THIS ROUTINE CLOSELY FOLLOWS THE ANALYTICAL DEVELOPMENTS DESCRIBED IN 
2597
CC    MARTI & MUELLER, JCP, 1996
2598

        
2599
      SUBROUTINE AVRG1D(N,DT,LAMB,PM,PP,DP,P6,RHOM,RHOP,DRHO,
2600
     &                  RHO6,VELM,VELP,DVEL,VEL6,PL,PR,RHOL,RHOR,
2601
     &                  VELL,VELR)
2602
 
2603
      IMPLICIT NONE
2604
 
2605
      INCLUDE 'size'
2606
 
2607
C     --------
2608
C     ARGUMENTS
2609
C     --------
2610

        
2611
      INTEGER         N
2612
 
2613
      DOUBLEPRECISION DT
2614
 
2615
      DOUBLEPRECISION LAMB(-4:MN5)
2616

        
2617
      DOUBLEPRECISION PM(-4:MN5),PP(-4:MN5),DP(-4:MN5),P6(-4:MN5),
2618
     &                RHOM(-4:MN5),RHOP(-4:MN5),DRHO(-4:MN5),
2619
     &                RHO6(-4:MN5),VELM(-4:MN5),VELP(-4:MN5),
2620
     &                DVEL(-4:MN5),VEL6(-4:MN5)
2621

        
2622
      DOUBLEPRECISION PL(-4:MN6),PR(-4:MN6),RHOL(-4:MN6),RHOR(-4:MN6),
2623
     &                VELL(-4:MN6),VELR(-4:MN6)
2624

        
2625
C     ------
2626
C     COMMON BLOCKS
2627
C     ------
2628

        
2629
      DOUBLEPRECISION X(-4:MN5),XL(-4:MN5),XR(-4:MN5),DX(-4:MN5)
2630
      COMMON /GRD/    X,XL,XR,DX
2631

        
2632
C     ----------
2633
C     INTERNAL VARIABLES
2634
C     ----------
2635

        
2636
      INTEGER         I
2637
 
2638
      DOUBLEPRECISION SCRCH1(-4:MN6)
2639
 
2640

        
2641
      DO 10 I=0,N+1
2642

        
2643
        SCRCH1(I) = DMAX1(0.D0,DT*LAMB(I)/DX(I))
2644

        
2645
10      CONTINUE
2646
 
2647
      DO 20 I=1,N+1
2648
        
2649
        PL(I)   = PP(I-1)   - SCRCH1(I-1)/2.D0*
2650
     &           (DP(I-1)   - (1.D0 - 2.D0*SCRCH1(I-1)/3.D0)*P6(I-1)  )
2651
 
2652
        RHOL(I) = RHOP(I-1) - SCRCH1(I-1)/2.D0*
2653
     &           (DRHO(I-1) - (1.D0 - 2.D0*SCRCH1(I-1)/3.D0)*RHO6(I-1))
2654
 
2655
        VELL(I) = VELP(I-1) - SCRCH1(I-1)/2.D0*
2656
     &           (DVEL(I-1) - (1.D0 - 2.D0*SCRCH1(I-1)/3.D0)*VEL6(I-1))
2657

        
2658
20      CONTINUE
2659
 
2660
      DO 30 I=0,N+1
2661

        
2662
        SCRCH1(I) = DMAX1(0.D0,-DT*LAMB(I)/DX(I))
2663

        
2664
30      CONTINUE
2665
 
2666
      DO 40 I=1,N+1
2667
        
2668
        PR(I)   = PM(I)   + SCRCH1(I)/2.D0*
2669
     &           (DP(I)   + (1.D0 - 2.D0*SCRCH1(I)/3.D0)*P6(I)  )
2670
 
2671
        RHOR(I) = RHOM(I) + SCRCH1(I)/2.D0*
2672
     &           (DRHO(I) + (1.D0 - 2.D0*SCRCH1(I)/3.D0)*RHO6(I))
2673
 
2674
        VELR(I) = VELM(I) + SCRCH1(I)/2.D0*
2675
     &           (DVEL(I) + (1.D0 - 2.D0*SCRCH1(I)/3.D0)*VEL6(I))
2676

        
2677
40      CONTINUE
2678

        
2679
      RETURN
2680
      END
2681

        
2682
C     -------- 
2683
CN    NAME: N F L U X 
2684
C     --------
2685

        
2686
CP    PURPOSE: 
2687
CP    COMPUTES THE NUMERICAL FLUXES 
2688
C
2689

        
2690
CC    COMMENTS: 
2691
CC    COMPUTES THE NUMERICAL FLUXES FROM THE EXACT SOLUTION OF THE 
2692
CC    RELATIVISTIC RIEMANN PROBLEM AS DESCRIBED IN MARTI AND MUELLER, JFM, 
2693
CC    1994
2694
      
2695
      SUBROUTINE NFLUX(RHOL1,RHOR1,PL1,PR1,VELL1,VELR1,UL1,UR1,
2696
     &                 CSL1,CSR1,FR,FM,FE)
2697

        
2698
      IMPLICIT NONE
2699

        
2700
C     -----------
2701
C     ARGUMENTS
2702
C     -----------
2703

        
2704
      DOUBLE PRECISION RHOL1, PL1, UL1, CSL1, VELL1, 
2705
     &                 RHOR1, PR1, UR1, CSR1, VELR1 
2706

        
2707
      DOUBLE PRECISION FR, FM, FE
2708

        
2709
C     ------- 
2710
C     COMMON BLOCKS 
2711
C     -------
2712

        
2713
      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, WL, 
2714
     &                 RHOR, PR, UR, HR, CSR, VELR, WR  
2715
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, WL, 
2716
     &                 RHOR, PR, UR, HR, CSR, VELR, WR
2717

        
2718
      DOUBLE PRECISION RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL 
2719
      COMMON /LS/      RHOLS, ULS, HLS, CSLS, VELLS, VSHOCKL
2720

        
2721
      DOUBLE PRECISION RHORS, URS, HRS, CSRS, VELRS, VSHOCKR 
2722
      COMMON /RS/      RHORS, URS, HRS, CSRS, VELRS, VSHOCKR
2723

        
2724
      DOUBLE PRECISION GAMMA 
2725
      COMMON /ADIND/   GAMMA
2726

        
2727
C     --------- 
2728
C     INTERNAL VARIABLES 
2729
C     ---------
2730

        
2731
      INTEGER          ILOOP
2732

        
2733
      DOUBLE PRECISION TOL, PMIN, PMAX, DVEL1, DVEL2, CHECK
2734

        
2735
      DOUBLE PRECISION PS, VELS
2736

        
2737
      DOUBLE PRECISION RHOA, PA, VELA, UA
2738

        
2739
      DOUBLE PRECISION XI, XI1, XI2, XI3, XI4, XI5
2740

        
2741
C     ------------------------------ 
2742
C     SPECIFIC ENTHALPY AND  
2743
C     FLOW LORENTZ FACTORS IN THE INITIAL STATES 
2744
C     ------------------------------ 
2745

        
2746
      RHOL = RHOL1
2747
      RHOR = RHOR1
2748

        
2749
      PL   = PL1
2750
      PR   = PR1
2751

        
2752
      UL   = UL1
2753
      UR   = UR1
2754

        
2755
      VELL = VELL1
2756
      VELR = VELR1
2757

        
2758
      CSL  = CSL1
2759
      CSR  = CSR1
2760

        
2761
      HL  = 1.D0+UL+PL/RHOL 
2762
      HR  = 1.D0+UR+PR/RHOR
2763

        
2764
      WL  = 1.D0/DSQRT(1.D0-VELL**2) 
2765
      WR  = 1.D0/DSQRT(1.D0-VELR**2)
2766

        
2767
C     ------------- 
2768
C     TOLERANCE FOR THE SOLUTION 
2769
C     -------------
2770

        
2771
      TOL = 1.D-8
2772

        
2773
C
2774

        
2775
      ILOOP = 0
2776

        
2777
      PMIN  = (PL + PR)/2.D0 
2778
      PMAX  = PMIN
2779

        
2780
 5    ILOOP = ILOOP + 1
2781

        
2782
      PMIN  = 0.5D0*MAX(PMIN,0.D0) 
2783
      PMAX  = 2.D0*PMAX
2784

        
2785
      CALL GETDVEL(PMIN, DVEL1)
2786

        
2787
      CALL GETDVEL(PMAX, DVEL2)
2788

        
2789
      CHECK = DVEL1*DVEL2 
2790
      IF (CHECK.GT.0.D0) GOTO 5
2791

        
2792
C     --------------------------- 
2793
C     PRESSURE AND FLOW VELOCITY IN THE INTERMEDIATE STATES 
2794
C     ---------------------------
2795

        
2796
      CALL GETP(PMIN, PMAX, TOL, PS)
2797

        
2798
      VELS = 0.5D0*(VELLS + VELRS)
2799

        
2800
C     --------------- 
2801
C     SOLUTION ON THE NUMERICAL INTERFACE 
2802
C     ---------------
2803

        
2804
C     ----------- 
2805
C     POSITIONS OF THE WAVES 
2806
C     -----------
2807

        
2808
      IF (PL.GE.PS) THEN
2809

        
2810
        XI1 = (VELL - CSL )/(1.D0 - VELL*CSL ) 
2811
        XI2 = (VELS - CSLS)/(1.D0 - VELS*CSLS)
2812

        
2813
      ELSE
2814

        
2815
        XI1 = VSHOCKL
2816
        XI2 = XI1
2817

        
2818
      END IF
2819

        
2820
        XI3 = VELS
2821

        
2822
      IF (PR.GE.PS) THEN
2823

        
2824
        XI4 = (VELS + CSRS)/(1.D0 + VELS*CSRS) 
2825
        XI5 = (VELR + CSR )/(1.D0 + VELR*CSR )
2826

        
2827
      ELSE
2828

        
2829
        XI4 = VSHOCKR
2830
        XI5 = XI4
2831

        
2832
      END IF
2833

        
2834
C     ---------- 
2835
C     SOLUTION ON THE INTERFACE AT X = 0 (XI = 0)
2836
C     ----------
2837

        
2838
      XI = 0.D0
2839

        
2840
      IF (XI1.GE.XI) THEN
2841

        
2842
        PA   = PL 
2843
        RHOA = RHOL 
2844
        VELA = VELL 
2845
        UA   = UL
2846

        
2847
      ELSE IF (XI2.GE.XI) THEN
2848

        
2849
        CALL RAREF(XI,RHOL,CSL,VELL,'L',RHOA,PA,UA,VELA)
2850

        
2851
      ELSE IF (XI3.GE.XI) THEN
2852

        
2853
        PA   = PS 
2854
        RHOA = RHOLS 
2855
        VELA = VELS 
2856
        UA   = ULS
2857

        
2858
      ELSE IF (XI4.GE.XI) THEN
2859

        
2860
        PA   = PS 
2861
        RHOA = RHORS 
2862
        VELA = VELS 
2863
        UA   = URS
2864

        
2865
      ELSE IF (XI5.GE.XI) THEN
2866

        
2867
        CALL RAREF(XI,RHOR,CSR,VELR,'R',RHOA,PA,UA,VELA)
2868

        
2869
      ELSE
2870

        
2871
        PA   = PR 
2872
        RHOA = RHOR 
2873
        VELA = VELR 
2874
        UA   = UR
2875

        
2876
      END IF
2877

        
2878
C     -----------
2879
C     NUMERICAL FLUXES 
2880
C     -----------
2881

        
2882
      FR = RHOA*VELA/DSQRT(1.D0 - VELA**2)  
2883
      FM = RHOA*(1.D0 + UA + PA/RHOA)*VELA**2/(1.D0 - VELA**2) + PA  
2884
      FE = RHOA*(1.D0 + UA + PA/RHOA)*VELA/(1.D0 - VELA**2) -
2885
     &     RHOA*VELA/DSQRT(1.D0 - VELA**2)
2886
   
2887
      RETURN
2888
      END
2889

        
2890
C     ---------- 
2891
CN    NAME: G E T D V E L 
2892
C     ----------
2893

        
2894
CP    PURPOSE: 
2895
CP    COMPUTE THE DIFFERENCE IN FLOW SPEED BETWEEN LEFT AND RIGHT INTERMEDIATE 
2896
CP    STATES FOR GIVEN LEFT AND RIGHT STATES AND PRESSURE 
2897
C 
2898

        
2899
CC    COMMENTS:
2900
CC    NONE
2901
C 
2902
      SUBROUTINE GETDVEL( P, DVEL )
2903

        
2904
      IMPLICIT NONE
2905

        
2906
C     ----- 
2907
C     ARGUMENTS 
2908
C     -----
2909

        
2910
      DOUBLEPRECISION P, DVEL
2911

        
2912
C     ------- 
2913
C     COMMON BLOCKS 
2914
C     -------
2915

        
2916
      DOUBLE PRECISION RHOLS,ULS,HLS,CSLS,VELLS,VSHOCKL 
2917
      COMMON /LS/      RHOLS,ULS,HLS,CSLS,VELLS,VSHOCKL
2918

        
2919
      DOUBLE PRECISION RHORS,URS,HRS,CSRS,VELRS,VSHOCKR 
2920
      COMMON /RS/      RHORS,URS,HRS,CSRS,VELRS,VSHOCKR
2921

        
2922
      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, WL, 
2923
     & RHOR, PR, UR, HR, CSR, VELR, WR 
2924
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, WL, 
2925
     & RHOR, PR, UR, HR, CSR, VELR, WR
2926

        
2927
      DOUBLE PRECISION GAMMA 
2928
      COMMON /ADIND/   GAMMA
2929

        
2930
C     ----- 
2931
C     LEFT WAVE 
2932
C     -----
2933

        
2934
      CALL GETVEL(P, RHOL, PL, UL,  HL,  CSL,  VELL,  WL, 'L', 
2935
     & RHOLS,    ULS, HLS, CSLS, VELLS, VSHOCKL )
2936

        
2937
C     ----- 
2938
C     RIGHT WAVE 
2939
C     -----
2940

        
2941
      CALL GETVEL(P, RHOR, PR, UR,  HR,  CSR,  VELR,  WR, 'R', 
2942
     & RHORS,    URS, HRS, CSRS, VELRS, VSHOCKR )
2943

        
2944
      DVEL = VELLS - VELRS
2945

        
2946
      RETURN 
2947
      END
2948

        
2949
C     ------- 
2950
CN    NAME: G E T P 
2951
C     -------
2952

        
2953
CP    PURPOSE: 
2954
CP    FIND THE PRESSURE IN THE INTERMEDIATE STATE OF A RIEMANN PROBLEM IN 
2955
CP    RELATIVISTIC HYDRODYNAMICS 
2956
C 
2957

        
2958
CC    COMMENTS: 
2959
CC    THIS ROUTINE USES A COMBINATION OF INTERVAL BISECTION AND INVERSE 
2960
CC    QUADRATIC INTERPOLATION TO FIND THE ROOT IN A SPECIFIED INTERVAL. 
2961
CC    IT IS ASSUMED THAT DVEL(PMIN) AND DVEL(PMAX) HAVE OPPOSITE SIGNS WITHOUT 
2962
CC    A CHECK. 
2963
CC    ADAPTED FROM "COMPUTER METHODS FOR MATHEMATICAL COMPUTATION", 
2964
CC    BY G. E. FORSYTHE, M. A. MALCOLM, AND C. B. MOLER, 
2965
CC    PRENTICE-HALL, ENGLEWOOD CLIFFS N.J. 
2966
C
2967
      SUBROUTINE GETP( PMIN, PMAX, TOL, PS )
2968

        
2969
      IMPLICIT NONE
2970

        
2971
C     ----- 
2972
C     ARGUMENTS 
2973
C     -----
2974

        
2975
      DOUBLEPRECISION PMIN, PMAX, TOL, PS
2976

        
2977
C     ------- 
2978
C     COMMON BLOCKS 
2979
C     -------
2980

        
2981
      DOUBLEPRECISION GAMMA 
2982
      COMMON /ADIND/  GAMMA
2983

        
2984
      DOUBLEPRECISION RHOL, PL, UL, HL, CSL, VELL, WL, 
2985
     & RHOR, PR, UR, HR, CSR, VELR, WR 
2986
      COMMON /STATES/ RHOL, PL, UL, HL, CSL, VELL, WL, 
2987
     & RHOR, PR, UR, HR, CSR, VELR, WR
2988

        
2989
C     --------- 
2990
C     INTERNAL VARIABLES 
2991
C     ---------
2992

        
2993
      DOUBLEPRECISION A, B, C, D, E, EPS, FA, FB, FC, TOL1, 
2994
     & XM, P, Q, R, S
2995

        
2996
C     ------------- 
2997
C     COMPUTE MACHINE PRECISION 
2998
C     -------------
2999

        
3000
      EPS  = 1.D0 
3001
10    EPS  = EPS/2.D0 
3002
      TOL1 = 1.D0 + EPS 
3003
      IF( TOL1 .GT. 1.D0 ) GO TO 10
3004

        
3005
C     ------- 
3006
C     INITIALIZATION 
3007
C     -------
3008

        
3009
      A  = PMIN 
3010
      B  = PMAX 
3011
      CALL GETDVEL(A,FA) 
3012
      CALL GETDVEL(B,FB)
3013

        
3014
C     ----- 
3015
C     BEGIN STEP 
3016
C     -----
3017

        
3018
20    C  = A 
3019
      FC = FA 
3020
      D  = B - A 
3021
      E  = D 
3022
30    IF( DABS(FC) .GE. DABS(FB) )GO TO 40 
3023
      A  = B 
3024
      B  = C 
3025
      C  = A 
3026
      FA = FB 
3027
      FB = FC 
3028
      FC = FA
3029

        
3030
C     -------- 
3031
C     CONVERGENCE TEST 
3032
C     --------
3033

        
3034
40    TOL1 = 2.D0*EPS*DABS(B) + 0.5D0*TOL 
3035
      XM   = 0.5D0*(C - B) 
3036
      IF( DABS(XM) .LE. TOL1 ) GO TO 90 
3037
      IF( FB .EQ. 0.D0 ) GO TO 90
3038

        
3039
C     ------------ 
3040
C     IS BISECTION NECESSARY? 
3041
C     ------------
3042

        
3043
      IF( DABS(E) .LT. TOL1 ) GO TO 70 
3044
      IF( DABS(FA) .LE. DABS(FB) ) GO TO 70
3045

        
3046
C     ------------------ 
3047
C     IS QUADRATIC INTERPOLATION POSSIBLE? 
3048
C     ------------------
3049

        
3050
      IF( A .NE. C ) GO TO 50
3051

        
3052
C     ---------- 
3053
C     LINEAR INTERPOLATION 
3054
C     ----------
3055

        
3056
      S = FB/FA 
3057
      P = 2.D0*XM*S 
3058
      Q = 1.D0 - S 
3059
      GO TO 60
3060

        
3061
C     ---------------- 
3062
C     INVERSE QUADRATIC INTERPOLATION 
3063
C     ----------------
3064

        
3065
50    Q = FA/FC 
3066
      R = FB/FC 
3067
      S = FB/FA 
3068
      P = S*(2.D0*XM*Q*(Q - R) - (B - A)*(R - 1.D0)) 
3069
      Q = (Q - 1.D0)*(R - 1.D0)*(S - 1.D0)
3070

        
3071
C     ------ 
3072
C     ADJUST SIGNS 
3073
C     ------
3074

        
3075
60    IF( P .GT. 0.D0 ) Q = -Q 
3076
      P = DABS(P)
3077

        
3078
C     -------------- 
3079
C     IS INTERPOLATION ACCEPTABLE? 
3080
C     --------------
3081

        
3082
      IF( (2.D0*P) .GE. (3.D0*XM*Q-DABS(TOL1*Q)) ) GO TO 70 
3083
      IF( P .GE. DABS(0.5D0*E*Q) ) GO TO 70 
3084
      E = D 
3085
      D = P/Q 
3086
      GO TO 80
3087

        
3088
C     ----- 
3089
C     BISECTION 
3090
C     -----
3091

        
3092
70    D = XM 
3093
      E = D
3094

        
3095
C     ------- 
3096
C     COMPLETE STEP 
3097
C     -------
3098

        
3099
80    A  = B 
3100
      FA = FB 
3101
      IF( DABS(D) .GT. TOL1 ) B = B+D 
3102
      IF( DABS(D) .LE. TOL1 ) B = B+DSIGN(TOL1,XM) 
3103
      CALL GETDVEL(B,FB) 
3104
      IF( (FB*(FC/DABS(FC))) .GT. 0.D0) GO TO 20 
3105
      GO TO 30
3106

        
3107
C     -- 
3108
C     DONE 
3109
C     --
3110

        
3111
90    PS = B
3112

        
3113
      RETURN 
3114
      END
3115

        
3116
C     --------- 
3117
CN    NAME: G E T V E L 
3118
C     ---------
3119

        
3120
CP    PURPOSE: 
3121
CP    COMPUTE THE FLOW VELOCITY BEHIND A RAREFACTION OR SHOCK IN TERMS OF THE 
3122
CP    POST-WAVE PRESSURE FOR A GIVEN STATE AHEAD THE WAVE IN A RELATIVISTIC 
3123
CP    FLOW 
3124
C 
3125

        
3126
CC    COMMENTS: 
3127
CC    THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN MARTI AND MUELLER, 
3128
CC    J. FLUID MECH., (1994)
3129

        
3130
      SUBROUTINE GETVEL( P, RHOA, PA, UA, HA, CSA, VELA, WA, S, 
3131
     & RHO,      U,  H,  CS,  VEL,  VSHOCK )
3132

        
3133
      IMPLICIT NONE
3134

        
3135
C     ----- 
3136
C     ARGUMENTS 
3137
C     -----
3138

        
3139
      DOUBLE PRECISION P, RHOA, PA, UA, HA, CSA, VELA, WA  
3140
      CHARACTER*1      S 
3141
      DOUBLE PRECISION RHO, U, H, CS, VEL, VSHOCK
3142

        
3143
C     ------- 
3144
C     COMMON BLOCKS 
3145
C     -------
3146

        
3147
      DOUBLE PRECISION GAMMA 
3148
      COMMON /ADIND/   GAMMA
3149

        
3150
C     --------- 
3151
C     INTERNAL VARIABLES 
3152
C     ---------
3153

        
3154
      DOUBLE PRECISION A, B, C, SIGN 
3155
      DOUBLE PRECISION J, WSHOCK 
3156
      DOUBLE PRECISION K, SQGL1
3157

        
3158
C     --------------- 
3159
C     LEFT OR RIGHT PROPAGATING WAVE 
3160
C     ---------------
3161

        
3162
      IF (S.EQ.'L') SIGN = -1.D0
3163

        
3164
      IF (S.EQ.'R') SIGN =  1.D0
3165

        
3166
C
3167

        
3168
      IF (P/PA - 1.D0.GT.1.D-10) THEN
3169

        
3170
C       --- 
3171
C       SHOCK 
3172
C       ---
3173

        
3174
        A  = 1.D0+(GAMMA-1.D0)*(PA-P)/GAMMA/P 
3175
        B  = 1.D0-A 
3176
        C  = HA*(PA-P)/RHOA-HA**2
3177

        
3178
C       ---------------- 
3179
C       CHECK FOR UNPHYSICAL ENTHALPIES 
3180
C       ----------------
3181

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

        
3185
C       ----------------------------- 
3186
C       SPECIFIC ENTHALPY IN THE POST-WAVE STATE 
3187
C       (FROM THE EQUATION OF STATE AND THE TAUB ADIABAT, 
3188
C       EQ.(74), MM94) 
3189
C       -----------------------------
3190

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

        
3193
C       --------------- 
3194
C       DENSITY IN THE POST-WAVE STATE 
3195
C       (FROM EQ.(73), MM94) 
3196
C       ---------------
3197

        
3198
        RHO = GAMMA*P/(GAMMA-1.D0)/(H-1.D0)
3199

        
3200
C       ------------------------ 
3201
C       SPECIFIC INTERNAL ENERGY IN THE POST-WAVE STATE 
3202
C       (FROM THE EQUATION OF STATE) 
3203
C       ------------------------
3204

        
3205
        U = P/(GAMMA-1.D0)/RHO
3206

        
3207
C       -------------------------- 
3208
C       MASS FLUX ACROSS THE WAVE  
3209
C       (FROM THE RANKINE-HUGONIOT RELATIONS, EQ.(71), MM94) 
3210
C       --------------------------
3211

        
3212
        J = SIGN*DSQRT((P-PA)/(HA/RHOA-H/RHO))
3213

        
3214
C       ---------- 
3215
C       SHOCK VELOCITY 
3216
C       (FROM EQ.(86), MM94 
3217
C       ----------
3218

        
3219
        A      = J**2+(RHOA*WA)**2 
3220
        B      = -VELA*RHOA**2*WA**2 
3221
        VSHOCK = (-B+SIGN*J**2*DSQRT(1.D0+RHOA**2/J**2))/A 
3222
        WSHOCK = 1.D0/DSQRT(1.D0-VSHOCK**2)
3223

        
3224
C       ------------------- 
3225
C       FLOW VELOCITY IN THE POST-SHOCK STATE 
3226
C       (FROM EQ.(67), MM94) 
3227
C       -------------------
3228

        
3229
        A = WSHOCK*(P-PA)/J+HA*WA*VELA 
3230
        B = HA*WA+(P-PA)*(WSHOCK*VELA/J+1.D0/RHOA/WA)
3231

        
3232
        VEL = A/B
3233

        
3234
C       --------------------- 
3235
C       LOCAL SOUND SPEED IN THE POST-SHOCK STATE 
3236
C       (FROM THE EQUATION OF STATE) 
3237
C       ---------------------
3238

        
3239
        CS = DSQRT(GAMMA*P/RHO/H)
3240

        
3241
      ELSE IF (P/PA - 1.D0.GT.0.D0) THEN
3242
       
3243
C       --------------
3244
C       ALMOST CONSTANT INTERMEDIATE STATE
3245
C       --------------
3246

        
3247
        RHO    = RHOA
3248
        U      = UA
3249
        H      = 1.D0 + U + P/RHO
3250
        CS     = DSQRT(GAMMA*P/RHO/H)
3251
        VEL    = VELA
3252
        VSHOCK = VELA
3253

        
3254
      ELSE
3255

        
3256
C       ------ 
3257
C       RAREFACTION 
3258
C       ------
3259

        
3260
C       --------------------------- 
3261
C       POLITROPIC CONSTANT OF THE GAS ACROSS THE RAREFACTION 
3262
C       ---------------------------
3263

        
3264
        K = PA/RHOA**GAMMA
3265

        
3266
C       --------------- 
3267
C       DENSITY BEHIND THE RAREFACTION 
3268
C       ---------------
3269

        
3270
        RHO = (P/K)**(1.D0/GAMMA)
3271

        
3272
C       ------------------------ 
3273
C       SPECIFIC INTERNAL ENERGY BEHIND THE RAREFACTION 
3274
C       (FROM THE EQUATION OF STATE) 
3275
C       ------------------------
3276

        
3277
        U = P/(GAMMA-1.D0)/RHO
3278

        
3279
C       -----------
3280
C       SPECIFIC ENTHALPY
3281
C       -----------
3282

        
3283
        H = 1.D0 + U + P/RHO
3284
C       -------------------- 
3285
C       LOCAL SOUND SPEED BEHIND THE RAREFACTION 
3286
C       (FROM THE EQUATION OF STATE) 
3287
C       --------------------
3288

        
3289
        CS = DSQRT(GAMMA*P/RHO/H)
3290

        
3291
C       ------------------ 
3292
C       FLOW VELOCITY BEHIND THE RAREFACTION 
3293
C       ------------------
3294

        
3295
        SQGL1 = DSQRT(GAMMA-1.D0) 
3296
        A   = (1.D0+VELA)/(1.D0-VELA)* 
3297
     & ((SQGL1+CSA)/(SQGL1-CSA)* 
3298
     & (SQGL1-CS )/(SQGL1+CS ))**(-SIGN*2.D0/SQGL1)
3299

        
3300
        VEL = (A-1.D0)/(A+1.D0)
3301

        
3302
      END IF
3303
      END
3304

        
3305
C     -------- 
3306
CN    NAME: R A R E F 
3307
C     --------
3308

        
3309
CP    PURPOSE: 
3310
CP    COMPUTE THE FLOW STATE IN A RAREFACTION FOR GIVEN PRE-WAVE STATE 
3311
C
3312

        
3313
CC    COMMENTS: 
3314
CC    THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN MARTI AND MUELLER, 
3315
CC    J. FLUID MECH., (1994)
3316

        
3317
      SUBROUTINE RAREF( XI, RHOA, CSA, VELA, S, RHO, P, U, VEL )
3318

        
3319
      IMPLICIT NONE
3320

        
3321
C     ----- 
3322
C     ARGUMENTS 
3323
C     -----
3324

        
3325
      DOUBLE PRECISION XI
3326

        
3327
      DOUBLE PRECISION RHOA, CSA, VELA
3328

        
3329
      CHARACTER        S
3330

        
3331
      DOUBLE PRECISION RHO, P, U, VEL
3332

        
3333
C     ------- 
3334
C     COMMON BLOCKS 
3335
C     -------
3336

        
3337
      DOUBLE PRECISION GAMMA 
3338
      COMMON /ADIND/   GAMMA
3339

        
3340
C     --------- 
3341
C     INTERNAL VARIABLES 
3342
C     ---------
3343

        
3344
      DOUBLE PRECISION B, C, D, K, L, V, OCS2, FCS2, DFDCS2, CS2, SIGN
3345

        
3346
C     --------------- 
3347
C     LEFT OR RIGHT PROPAGATING WAVE 
3348
C     ---------------
3349

        
3350
      IF (S.EQ.'L') SIGN =  1.D0
3351

        
3352
      IF (S.EQ.'R') SIGN = -1.D0
3353

        
3354
      B    = DSQRT(GAMMA - 1.D0) 
3355
      C    = (B + CSA)/(B - CSA) 
3356
      D    = -SIGN*B/2.D0 
3357
      K    = (1.D0 + XI)/(1.D0 - XI) 
3358
      L    = C*K**D 
3359
      V    = ((1.D0 - VELA)/(1.D0 + VELA))**D
3360

        
3361
      OCS2 = CSA
3362

        
3363
25    FCS2   = L*V*(1.D0 + SIGN*OCS2)**D*(OCS2 - B) + 
3364
     & (1.D0 - SIGN*OCS2)**D*(OCS2 + B)
3365

        
3366
      DFDCS2 = L*V*(1.D0 + SIGN*OCS2)**D* 
3367
     & (1.D0 + SIGN*D*(OCS2 - B)/(1.D0 + SIGN*OCS2)) + 
3368
     & (1.D0 - SIGN*OCS2)**D*
3369
     & (1.D0 - SIGN*D*(OCS2 + B)/(1.D0 - SIGN*OCS2))
3370

        
3371
      CS2 = OCS2 - FCS2/DFDCS2
3372

        
3373
      IF (ABS(CS2 - OCS2)/OCS2.GT.5.E-7)THEN 
3374
        OCS2 = CS2 
3375
        GOTO 25 
3376
      END IF
3377

        
3378
      VEL = (XI + SIGN*CS2)/(1.D0 + SIGN*XI*CS2)
3379

        
3380
      RHO = RHOA*((CS2**2*(GAMMA - 1.D0 - CSA**2))/ 
3381
     & (CSA**2*(GAMMA - 1.D0 - CS2**2))) 
3382
     & **(1.D0/(GAMMA - 1.D0))
3383

        
3384
      P   = CS2**2*(GAMMA - 1.D0)*RHO/(GAMMA - 1.D0 - CS2**2)/GAMMA
3385

        
3386
      U   = P/(GAMMA - 1.D0)/RHO
3387

        
3388
      RETURN  
3389
      END 
3390

        
3391
C     -------- 
3392
CN    NAME: F I L N A M
3393
C     --------
3394

        
3395
CP    PURPOSE: 
3396
CP    CONSTRUCTS NEW FILENAMES FOR OUTPUT AND RESTART FILES
3397
C
3398

        
3399
CC    COMMENTS: 
3400
CC    NONE
3401

        
3402
      SUBROUTINE FILNAM
3403

        
3404
      IMPLICIT NONE
3405

        
3406
C     ---------
3407
C     COMMON BLOCKS
3408
C     ---------
3409

        
3410
      CHARACTER*7   OUTFIL
3411
      CHARACTER*8   LABEL
3412
      CHARACTER*4   BASENM
3413
      CHARACTER*2   SUFFIX
3414
      CHARACTER*1   SF1,SF2
3415

        
3416
      COMMON /CHRC/ LABEL,OUTFIL,BASENM,SUFFIX
3417

        
3418
C     ---------
3419
C     INTERNAL VARIABLES
3420
C     ---------
3421

        
3422
      INTEGER       ISF1,ISF2
3423

        
3424
      IF (SUFFIX(2:2).EQ.'z'.OR.SUFFIX(2:2).EQ.'Z') THEN
3425
         SF1  = SUFFIX(1:1)
3426
         ISF1 = ICHAR(SF1)
3427
         SF2  = SUFFIX(2:2)
3428
         ISF2 = ICHAR(SF2)
3429

        
3430
         ISF1 = ISF1 + 1
3431
         ISF2 = ISF2 - 25
3432
         SUFFIX(1:1) = CHAR(ISF1)
3433
         SUFFIX(2:2) = CHAR(ISF2)
3434
      ELSE
3435
         SF2  = SUFFIX(2:2)
3436
         ISF2 = ICHAR(SF2)
3437
         ISF2 = ISF2 + 1
3438
         SUFFIX(2:2) = CHAR(ISF2)
3439
      END IF
3440

        
3441
      OUTFIL = BASENM // 'O' // SUFFIX
3442

        
3443
      RETURN
3444
      END
3445

        
3446
C     -------- 
3447
CN    NAME: G E T P R F Q
3448
C     --------
3449

        
3450
CP    PURPOSE: 
3451
CP    COMPUTE THE PRIMITIVE QUANTITIES FROM THE CONSERVED ONES  
3452
C
3453

        
3454
CC    COMMENTS: 
3455
CC    PRIMITIVE VARIABLES ARE OBTAINED BY SOLVING AN IMPLICIT EQUATION FOR 
3456
CC    THE PRESSURE BY MEANS OF A NEWTON-RAPHSON METHOD. HARDWIRED FOR A 
3457
CC    CONSTANT-GAMMA IDEAL GAS
3458

        
3459
      SUBROUTINE GETPRFQ(N,R,M,E,
3460
     &                   VEL,W,RHO,U,P,H,CS,DPDRH,DPDU)
3461
 
3462
      IMPLICIT NONE
3463

        
3464
      INCLUDE 'size'
3465

        
3466
C     -------
3467
C     ARGUMENTS
3468
C     -------
3469

        
3470
      INTEGER N
3471

        
3472
      DOUBLEPRECISION P(-4:MN5),RHO(-4:MN5),VEL(-4:MN5),
3473
     &                U(-4:MN5),CS(-4:MN5),W(-4:MN5),H(-4:MN5),
3474
     &                DPDRH(-4:MN5),DPDU(-4:MN5)
3475

        
3476
      DOUBLEPRECISION R(-4:MN5),M(-4:MN5),E(-4:MN5)
3477

        
3478
C     ---------
3479
C     COMMON BLOCKS
3480
C     ---------
3481

        
3482
      INTEGER         NEND,NOUT,ITSTP,NX
3483
      COMMON /INPTI/  NEND,NOUT,ITSTP,NX
3484

        
3485
      DOUBLEPRECISION TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,
3486
     &                SMALLU,GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
3487
      COMMON /INPTF/  TMAX,TOUT,CFL,DTINI,SMALL,SMLRHO,SMALLP,
3488
     &                SMALLU,GRIDLX,EPSILN,ETA1,ETA2,EPSLN,AK0,OMG1,OMG2
3489

        
3490
      DOUBLEPRECISION GB
3491
      COMMON /ADIND/  GB
3492

        
3493
      DOUBLEPRECISION TIME,DT
3494
      COMMON /ZEIT/   TIME,DT
3495

        
3496
C     ---------
3497
C     INTERNAL VARIABLES
3498
C     ---------
3499

        
3500
      INTEGER I,COUNT
3501

        
3502
      DOUBLEPRECISION MOMENT,VELCTY
3503

        
3504
      DOUBLEPRECISION PMIN(-4:MN5),PMAX,OP(-4:MN5)
3505

        
3506
      DOUBLEPRECISION FP,DFDP,ERRP
3507

        
3508
      DO 5 I=1,N
3509

        
3510
         R(I) = DMAX1(R(I),SMLRHO)
3511
         E(I) = DMAX1(E(I),SMALLU)
3512

        
3513
 5       CONTINUE
3514

        
3515
      DO 9 I=1,N
3516

        
3517
         COUNT   = 0
3518

        
3519
         MOMENT  = DABS(M(I))
3520
         PMIN(I) = DMAX1(MOMENT - E(I) - R(I) + MOMENT*SMALL,SMALLP)
3521
         PMAX    = (GB-1.D0)*E(I)
3522
         IF (PMIN(I).GT.PMAX) GOTO 990
3523

        
3524
         OP(I)      = 0.5D0*(PMIN(I)+PMAX)
3525

        
3526
 8       COUNT   = COUNT + 1
3527
         OP(I)   = DMAX1(OP(I),PMIN(I))
3528
         VELCTY  = MOMENT/(E(I) + R(I) + OP(I))
3529
         W(I)    = 1.D0/DSQRT(1.D0 - VELCTY*VELCTY)
3530
         FP      = (GB - 1.D0)*(E(I) + R(I)*(1.D0 - W(I))+
3531
     &             OP(I)*(1.D0 - W(I)*W(I)))/W(I)/W(I) - OP(I)
3532
         DFDP    = (GB - 1.D0)*VELCTY*VELCTY*
3533
     &             (E(I) + R(I)*(1.D0 - W(I))+OP(I))/
3534
     &             (E(I) + R(I) + OP(I)) - 1.D0
3535

        
3536
         P(I)    = DMAX1(OP(I) - FP/DFDP,PMIN(I))
3537

        
3538
         ERRP    = DABS(1.D0 - P(I)/OP(I))
3539

        
3540
         OP(I)      = P(I)
3541

        
3542
         IF (COUNT.GE.10000) GOTO 999
3543

        
3544
         IF (ERRP.GT.1.D-8) GOTO 8
3545

        
3546
         VEL(I)  = M(I)/(E(I)+R(I)+OP(I))
3547
         IF (DABS(VEL(I)) .LT.SMALL*SMALL) VEL(I)  = 0.D0
3548
 
3549
         RHO(I)  = R(I)/W(I)
3550

        
3551
 9       CONTINUE
3552

        
3553
      DO 30 I=1,N
3554

        
3555
	 U(I) = P(I)/(GB - 1.D0)/RHO(I)
3556

        
3557
         IF (P(I).EQ.PMIN(I)) THEN
3558
            WRITE(6,*) 'GETPRFQ: MINIMUM PRESSURE REACHED AT POINT:'
3559
            WRITE(6,*) '         I = ', I,' T = ', TIME
3560
         END IF
3561

        
3562
 30      CONTINUE
3563

        
3564
      CALL EOS(N,RHO,U,GB,P,H,CS,DPDRH,DPDU)
3565

        
3566
      GOTO 1000
3567

        
3568
 990  WRITE(6,*) 'GETPRFQ: NO PHYSICAL PRESSURE AVAILABLE'
3569
      WRITE(6,*) '         T          = ', TIME
3570
      WRITE(6,*) '         I          = ', I
3571
      WRITE(6,*) '         R          = ', R(I),   ' MOMENT = ', MOMENT 
3572
      WRITE(6,*) '         E          = ', E(I)
3573
      WRITE(6,*) '         MOMENT-E-D = ', MOMENT - R(I) - E(I)
3574
      WRITE(6,*) '         (GB-1)E    = ', (GB - 1.D0)*E(I) 
3575
      STOP
3576

        
3577
 999  WRITE(6,*) 'GETPRFQ: NON CONVERGENCE IN PRESSURE'
3578
      WRITE(6,*) '         T   = ', TIME
3579
      WRITE(6,*) '         I   = ', I
3580
      WRITE(6,*) '         P   = ', P(I),   ' PMIN = ', PMIN(I)
3581
      WRITE(6,*) '         R   = ', R(I),   ' M    = ', M(I)
3582
      WRITE(6,*) '         E   = ', E(I)
3583
      WRITE(6,*) '         VEL = ', VEL(I)
3584
      STOP
3585

        
3586
 1000 CONTINUE
3587

        
3588
      RETURN
3589
      END