text:
      PROGRAM PCNEPDRIVER
*
*  -- ScaLAPACK testing driver (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     March, 2000
*
*  Purpose
*  =======
*
*  PCNEPDRIVER is the main test program for the COMPLEX
*  SCALAPACK NEP routines.  This test driver performs a Schur
*  decomposition followed by residual check of a Hessenberg matrix.
*
*  The program must be driven by a short data file.  An annotated
*  example of a data file can be obtained by deleting the first 3
*  characters from the following 18 lines:
*  'SCALAPACK, Version 1.4, NEP (Nonsymmetric EigenProblem) input file'
*  'Intel iPSC/860 hypercube, gamma model.'
*  'NEP.out'            output file name (if any)
*  6                    device out
*  8                    number of problems sizes
*  1 2 3 4 6 10 100 200 vales of N
*  3                    number of NB's
*  6 20 40              values of NB
*  4                    number of process grids (ordered pairs of P & Q)
*  1 2 1 4              values of P
*  1 2 4 1              values of Q
*  20.0                 threshold
*
*  Internal Parameters
*  ===================
*
*  TOTMEM   INTEGER, default = 2000000
*           TOTMEM is a machine-specific parameter indicating the
*           maximum amount of available memory in bytes.
*           The user should customize TOTMEM to his platform.  Remember
*           to leave room in memory for the operating system, the BLACS
*           buffer, etc.  For example, on a system with 8 MB of memory
*           per process (e.g., one processor on an Intel iPSC/860), the
*           parameters we use are TOTMEM=6200000 (leaving 1.8 MB for OS,
*           code, BLACS buffer, etc).  However, for PVM, we usually set
*           TOTMEM = 2000000.  Some experimenting with the maximum value
*           of TOTMEM may be required.
*
*  CPLXSZ   INTEGER, default = 8 bytes.
*           CPLXSZ indicate the length in bytes on the given platform
*           for a single precision complex.
*  MEM      COMPLEX array, dimension ( TOTMEM / CPLXSZ )
*
*           All arrays used by SCALAPACK routines are allocated from
*           this array and referenced by pointers.  The integer IPA,
*           for example, is a pointer to the starting element of MEM for
*           the matrix A.
*
*  Further Details
*  ===============
*
*  Contributed by Mark Fahey, March 2000.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DT_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DT_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      INTEGER            CPLXSZ, TOTMEM, MEMSIZ, NTESTS
      PARAMETER          ( CPLXSZ = 16, TOTMEM = 200000000,
     $                   MEMSIZ = TOTMEM / CPLXSZ, NTESTS = 20 )
      COMPLEX            PADVAL, ZERO, ONE
      PARAMETER          ( PADVAL = ( -9923.0E+0, -9923.0E+0 ),
     $                   ZERO = ( 0.0E+0, 0.0E+0 ),
     $                   ONE = ( 1.0E+0, 0.0E+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            CHECK
      CHARACTER*6        PASSED
      CHARACTER*80       OUTFILE
      INTEGER            I, IAM, IASEED, ICTXT, III, IMIDPAD, INFO, IPA,
     $                   IPOSTPAD, IPREPAD, IPW, IPWR, IPC, J, K, KFAIL,
     $                   KPASS, KSKIP, KTESTS, LDA, LDWORK, LDZ, LWORK,
     $                   MYCOL, MYROW, N, NB, NGRIDS, NMAT, NNB, NOUT,
     $                   NP, NPCOL, NPROCS, NPROW, NQ, WORKSIZ
      REAL               THRESH
      REAL               ANORM, FRESID, QRESID, ZNORM
      DOUBLE PRECISION   NOPS, TMFLOPS
*     ..
*     .. Local Arrays ..
      INTEGER            DESCA( DLEN_ ), DESCZ( DLEN_ ), IDUM( 1 ),
     $                   IERR( 2 ), NBVAL( NTESTS ), NVAL( NTESTS ),
     $                   PVAL( NTESTS ), QVAL( NTESTS )
      DOUBLE PRECISION   CTIME( 2 ), WTIME( 2 )
      COMPLEX            MEM( MEMSIZ )
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_BARRIER, BLACS_EXIT, BLACS_GET,
     $                   BLACS_GRIDEXIT, BLACS_GRIDINFO, BLACS_GRIDINIT,
     $                   BLACS_PINFO, DESCINIT, IGSUM2D, PCCHEKPAD,
     $                   PCFILLPAD, PCGEMM, PCLAHQR, PCLASET, PCMATGEN,
     $                   PCNEPFCHK, PCNEPINFO, SLBOOT, SLCOMBINE,
     $                   SLTIMER
*     ..
*     .. External Functions ..
      INTEGER            ILCM, NUMROC
      REAL               PSLAMCH, PCLANGE, PCLANHS
      EXTERNAL           ILCM, NUMROC, PSLAMCH, PCLANGE, PCLANHS
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          REAL, MAX, MIN
*     ..
*     .. Data statements ..
      DATA               KFAIL, KPASS, KSKIP, KTESTS / 4*0 /
*     ..
*     .. Executable Statements ..
*
*     Get starting information
*
      CALL BLACS_PINFO( IAM, NPROCS )
      IASEED = 100
      CALL PCNEPINFO( OUTFILE, NOUT, NMAT, NVAL, NTESTS, NNB, NBVAL,
     $                NTESTS, NGRIDS, PVAL, NTESTS, QVAL, NTESTS,
     $                THRESH, MEM, IAM, NPROCS )
      CHECK = ( THRESH.GE.0.0E+0 )
*
*     Print headings
*
      IF( IAM.EQ.0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = 9995 )
         WRITE( NOUT, FMT = 9994 )
         WRITE( NOUT, FMT = * )
      END IF
*
*     Loop over different process grids
*
      DO 30 I = 1, NGRIDS
*
         NPROW = PVAL( I )
         NPCOL = QVAL( I )
*
*        Make sure grid information is correct
*
         IERR( 1 ) = 0
         IF( NPROW.LT.1 ) THEN
            IF( IAM.EQ.0 )
     $         WRITE( NOUT, FMT = 9999 )'GRID', 'nprow', NPROW
            IERR( 1 ) = 1
         ELSE IF( NPCOL.LT.1 ) THEN
            IF( IAM.EQ.0 )
     $         WRITE( NOUT, FMT = 9999 )'GRID', 'npcol', NPCOL
            IERR( 1 ) = 1
         ELSE IF( NPROW*NPCOL.GT.NPROCS ) THEN
            IF( IAM.EQ.0 )
     $         WRITE( NOUT, FMT = 9998 )NPROW*NPCOL, NPROCS
            IERR( 1 ) = 1
         END IF
*
         IF( IERR( 1 ).GT.0 ) THEN
            IF( IAM.EQ.0 )
     $         WRITE( NOUT, FMT = 9997 )'grid'
            KSKIP = KSKIP + 1
            GO TO 30
         END IF
*
*        Define process grid
*
         CALL BLACS_GET( -1, 0, ICTXT )
         CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
         CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*        Go to bottom of process grid loop if this case doesn't use my
*        process
*
         IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL )
     $      GO TO 30
*
         DO 20 J = 1, NMAT
*
            N = NVAL( J )
*
*           Make sure matrix information is correct
*
            IERR( 1 ) = 0
            IF( N.LT.1 ) THEN
               IF( IAM.EQ.0 )
     $            WRITE( NOUT, FMT = 9999 )'MATRIX', 'N', N
               IERR( 1 ) = 1
            END IF
*
*           Check all processes for an error
*
            CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, -1, 0 )
*
            IF( IERR( 1 ).GT.0 ) THEN
               IF( IAM.EQ.0 )
     $            WRITE( NOUT, FMT = 9997 )'matrix'
               KSKIP = KSKIP + 1
               GO TO 20
            END IF
*
            DO 10 K = 1, NNB
*
               NB = NBVAL( K )
*
*              Make sure nb is legal
*
               IERR( 1 ) = 0
               IF( NB.LT.6 ) THEN
                  IERR( 1 ) = 1
                  IF( IAM.EQ.0 )
     $               WRITE( NOUT, FMT = 9999 )'NB', 'NB', NB
               END IF
*
*              Check all processes for an error
*
               CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, -1, 0 )
*
               IF( IERR( 1 ).GT.0 ) THEN
                  IF( IAM.EQ.0 )
     $               WRITE( NOUT, FMT = 9997 )'NB'
                  KSKIP = KSKIP + 1
                  GO TO 10
               END IF
*
*              Padding constants
*
               NP = NUMROC( N, NB, MYROW, 0, NPROW )
               NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )
               IF( CHECK ) THEN
                  IPREPAD = MAX( NB, NP )
                  IMIDPAD = NB
                  IPOSTPAD = MAX( NB, NQ )
                  IPREPAD = IPREPAD + 1000
                  IMIDPAD = IMIDPAD + 1000
                  IPOSTPAD = IPOSTPAD + 1000
               ELSE
                  IPREPAD = 0
                  IMIDPAD = 0
                  IPOSTPAD = 0
               END IF
*
*              Initialize the array descriptor for the matrix A
*
               CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, ICTXT,
     $                        MAX( 1, NP )+IMIDPAD, IERR( 1 ) )
*
*              Initialize the array descriptor for the matrix Z
*
               CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, ICTXT,
     $                        MAX( 1, NP )+IMIDPAD, IERR( 2 ) )
*
               LDA = DESCA( LLD_ )
               LDZ = DESCZ( LLD_ )
               LDWORK = DESCZ( LLD_ )
*
*              Check all processes for an error
*
               CALL IGSUM2D( ICTXT, 'All', ' ', 2, 1, IERR, 2, -1, 0 )
*
               IF( IERR( 1 ).LT.0 .OR. IERR( 2 ).LT.0 ) THEN
                  IF( IAM.EQ.0 )
     $               WRITE( NOUT, FMT = 9997 )'descriptor'
                  KSKIP = KSKIP + 1
                  GO TO 10
               END IF
*
*              Assign pointers into MEM for SCALAPACK arrays, A is
*              allocated starting at position MEM( IPREPAE+1 )
*
               IPA = IPREPAD + 1
               IPC = IPA + DESCA( LLD_ )*NQ + IPOSTPAD + IPREPAD
               IPWR = IPC + DESCZ( LLD_ )*NQ + IPOSTPAD + IPREPAD
               IPW = IPWR + DESCZ( LLD_ )*NQ + IPOSTPAD + IPREPAD
               III = N / NB
               IF( III*NB.LT.N )
     $            III = III + 1
               III = 7*III / ILCM( NPROW, NPCOL )
*
*
               LWORK = 3*N + MAX( 2*MAX( LDA, LDZ )+2*NQ, III )
               LWORK = LWORK + MAX( 2*N, ( 8*ILCM( NPROW, NPCOL )+2 )**
     $                 2 )
*
               IF( CHECK ) THEN
*
*                 Figure the amount of workspace required by the
*                 checking routines PCNEPFCHK and PCLANHS
*
                  WORKSIZ = LWORK + MAX( NP*DESCA( NB_ ),
     $                      DESCA( MB_ )*NQ ) + IPOSTPAD
*
               ELSE
*
                  WORKSIZ = LWORK + IPOSTPAD
*
               END IF
*
*              Check for adequate memory for problem size
*
               IERR( 1 ) = 0
               IF( IPW+WORKSIZ.GT.MEMSIZ ) THEN
                  IF( IAM.EQ.0 )
     $               WRITE( NOUT, FMT = 9996 )'Schur reduction',
     $               ( IPW+WORKSIZ )*CPLXSZ
                  IERR( 1 ) = 1
               END IF
*
*              Check all processes for an error
*
               CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, -1, 0 )
*
               IF( IERR( 1 ).GT.0 ) THEN
                  IF( IAM.EQ.0 )
     $               WRITE( NOUT, FMT = 9997 )'MEMORY'
                  KSKIP = KSKIP + 1
                  GO TO 10
               END IF
*
*              Generate matrix Z = In
*
               CALL PCLASET( 'All', N, N, ZERO, ONE, MEM( IPC ), 1, 1,
     $                       DESCZ )
*
*              Generate matrix A upper Hessenberg
*
               CALL PCMATGEN( ICTXT, 'No transpose', 'No transpose',
     $                        DESCA( M_ ), DESCA( N_ ), DESCA( MB_ ),
     $                        DESCA( NB_ ), MEM( IPA ), DESCA( LLD_ ),
     $                        DESCA( RSRC_ ), DESCA( CSRC_ ), IASEED, 0,
     $                        NP, 0, NQ, MYROW, MYCOL, NPROW, NPCOL )
               CALL PCLASET( 'Lower', MAX( 0, N-2 ), MAX( 0, N-2 ),
     $                       ZERO, ZERO, MEM( IPA ), MIN( N, 3 ), 1,
     $                       DESCA )
*
*              Calculate inf-norm of A for residual error-checking
*
               IF( CHECK ) THEN
                  CALL PCFILLPAD( ICTXT, NP, NQ, MEM( IPA-IPREPAD ),
     $                            DESCA( LLD_ ), IPREPAD, IPOSTPAD,
     $                            PADVAL )
                  CALL PCFILLPAD( ICTXT, NP, NQ, MEM( IPC-IPREPAD ),
     $                            DESCZ( LLD_ ), IPREPAD, IPOSTPAD,
     $                            PADVAL )
                  CALL PCFILLPAD( ICTXT, WORKSIZ-IPOSTPAD, 1,
     $                            MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD,
     $                            IPREPAD, IPOSTPAD, PADVAL )
                  ANORM = PCLANHS( 'I', N, MEM( IPA ), 1, 1, DESCA,
     $                    MEM( IPW ) )
                  CALL PCCHEKPAD( ICTXT, 'PCLANHS', NP, NQ,
     $                            MEM( IPA-IPREPAD ), DESCA( LLD_ ),
     $                            IPREPAD, IPOSTPAD, PADVAL )
                  CALL PCCHEKPAD( ICTXT, 'PCLANHS', WORKSIZ-IPOSTPAD, 1,
     $                            MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD,
     $                            IPREPAD, IPOSTPAD, PADVAL )
*
                  CALL PCFILLPAD( ICTXT, N, 1, MEM( IPWR-IPREPAD ), N,
     $                            IPREPAD, IPOSTPAD, PADVAL )
                  CALL PCFILLPAD( ICTXT, LWORK, 1, MEM( IPW-IPREPAD ),
     $                            LWORK, IPREPAD, IPOSTPAD, PADVAL )
*
               END IF
*
               CALL SLBOOT()
               CALL BLACS_BARRIER( ICTXT, 'All' )
               CALL SLTIMER( 1 )
*
*              Perform NEP factorization
*
               CALL PCLAHQR( .TRUE., .TRUE., N, 1, N, MEM( IPA ), DESCA,
     $                       MEM( IPWR ), 1, N, MEM( IPC ), DESCZ,
     $                       MEM( IPW ), LWORK, IDUM, 0, INFO )
*
               CALL SLTIMER( 1 )
*
               IF( INFO.NE.0 ) THEN
                  IF( IAM.EQ.0 )
     $               WRITE( NOUT, FMT = * )'PCLAHQR INFO=', INFO
                  KFAIL = KFAIL + 1
                  GO TO 10
               END IF
*
               IF( CHECK ) THEN
*
*                 Check for memory overwrite in NEP factorization
*
                  CALL PCCHEKPAD( ICTXT, 'PCLAHQR (A)', NP, NQ,
     $                            MEM( IPA-IPREPAD ), DESCA( LLD_ ),
     $                            IPREPAD, IPOSTPAD, PADVAL )
                  CALL PCCHEKPAD( ICTXT, 'PCLAHQR (Z)', NP, NQ,
     $                            MEM( IPC-IPREPAD ), DESCZ( LLD_ ),
     $                            IPREPAD, IPOSTPAD, PADVAL )
                  CALL PCCHEKPAD( ICTXT, 'PCLAHQR (WR)', N, 1,
     $                            MEM( IPWR-IPREPAD ), N, IPREPAD,
     $                            IPOSTPAD, PADVAL )
                  CALL PCCHEKPAD( ICTXT, 'PCLAHQR (WORK)', LWORK, 1,
     $                            MEM( IPW-IPREPAD ), LWORK, IPREPAD,
     $                            IPOSTPAD, PADVAL )
*
                  CALL PCFILLPAD( ICTXT, WORKSIZ-IPOSTPAD, 1,
     $                            MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD,
     $                            IPREPAD, IPOSTPAD, PADVAL )
*
*                 Compute || Z * H * Z**T - H0 || / ( N*|| H0 ||*EPS )
*
                  CALL PCNEPFCHK( N, MEM( IPA ), 1, 1, DESCA, IASEED,
     $                            MEM( IPC ), 1, 1, DESCZ, ANORM,
     $                            FRESID, MEM( IPW ) )
*
                  CALL PCCHEKPAD( ICTXT, 'PCNEPFCHK (A)', NP, NQ,
     $                            MEM( IPA-IPREPAD ), DESCA( LLD_ ),
     $                            IPREPAD, IPOSTPAD, PADVAL )
                  CALL PCCHEKPAD( ICTXT, 'PCNEPFCHK (Z)', NP, NQ,
     $                            MEM( IPC-IPREPAD ), DESCZ( LLD_ ),
     $                            IPREPAD, IPOSTPAD, PADVAL )
                  CALL PCCHEKPAD( ICTXT, 'PCNEPFCHK (WORK)',
     $                            WORKSIZ-IPOSTPAD, 1,
     $                            MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD,
     $                            IPREPAD, IPOSTPAD, PADVAL )
*
*                 Compute || (Z**T)*Z - In ||_1
*
                  CALL PCLASET( 'All', N, N, ZERO, ONE, MEM( IPA ), 1,
     $                          1, DESCA )
                  CALL PCGEMM( 'Cong Tran', 'No transpose', N, N, N,
     $                         -ONE, MEM( IPC ), 1, 1, DESCZ,
     $                         MEM( IPC ), 1, 1, DESCZ, ONE, MEM( IPA ),
     $                         1, 1, DESCA )
                  ZNORM = PCLANGE( '1', N, N, MEM( IPA ), 1, 1, DESCA,
     $                    MEM( IPW ) )
                  QRESID = ZNORM / ( REAL( N )*PSLAMCH( ICTXT, 'P' ) )
*
*                 Test residual and detect NaN result
*
                  IF( ( FRESID.LE.THRESH ) .AND.
     $                ( ( FRESID-FRESID ).EQ.0.0E+0 ) .AND.
     $                ( QRESID.LE.THRESH ) .AND.
     $                ( ( QRESID-QRESID ).EQ.0.0E+0 ) ) THEN
                     KPASS = KPASS + 1
                     PASSED = 'PASSED'
                  ELSE
                     KFAIL = KFAIL + 1
                     PASSED = 'FAILED'
                     IF( IAM.EQ.0 ) THEN
                        WRITE( NOUT, FMT = 9986 )FRESID
                        WRITE( NOUT, FMT = 9985 )QRESID
                     END IF
                  END IF
*
               ELSE
*
*                 Don't perform the checking, only timing
*
                  KPASS = KPASS + 1
                  FRESID = FRESID - FRESID
                  QRESID = QRESID - QRESID
                  PASSED = 'BYPASS'
*
               END IF
*
*              Gather maximum of all CPU and WALL clock timings
*
               CALL SLCOMBINE( ICTXT, 'All', '>', 'W', 1, 1, WTIME )
               CALL SLCOMBINE( ICTXT, 'All', '>', 'C', 1, 1, CTIME )
*
*              Print results
*
               IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
*
*                 18 N^3 flops for PxLAHQR
*
                  NOPS = 18.0D+0*DBLE( N )**3
*
*                 Calculate total megaflops -- factorization only,
*                 -- for WALL and CPU time, and print output
*
*                 Print WALL time if machine supports it
*
                  IF( WTIME( 1 ).GT.0.0D+0 ) THEN
                     TMFLOPS = NOPS / ( WTIME( 1 )*1.0D+6 )
                  ELSE
                     TMFLOPS = 0.0D+0
                  END IF
                  IF( WTIME( 1 ).GE.0.0D+0 )
     $               WRITE( NOUT, FMT = 9993 )'WALL', N, NB, NPROW,
     $               NPCOL, WTIME( 1 ), TMFLOPS, PASSED
*
*                 Print CPU time if machine supports it
*
                  IF( CTIME( 1 ).GT.0.0D+0 ) THEN
                     TMFLOPS = NOPS / ( CTIME( 1 )*1.0D+6 )
                  ELSE
                     TMFLOPS = 0.0D+0
                  END IF
*
                  IF( CTIME( 1 ).GE.0.0D+0 )
     $               WRITE( NOUT, FMT = 9993 )'CPU ', N, NB, NPROW,
     $               NPCOL, CTIME( 1 ), TMFLOPS, PASSED
               END IF
*
   10       CONTINUE
*
   20    CONTINUE
*
         CALL BLACS_GRIDEXIT( ICTXT )
*
   30 CONTINUE
*
*     Print ending messages and close output file
*
      IF( IAM.EQ.0 ) THEN
         KTESTS = KPASS + KFAIL + KSKIP
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = 9992 )KTESTS
         IF( CHECK ) THEN
            WRITE( NOUT, FMT = 9991 )KPASS
            WRITE( NOUT, FMT = 9989 )KFAIL
         ELSE
            WRITE( NOUT, FMT = 9990 )KPASS
         END IF
         WRITE( NOUT, FMT = 9988 )KSKIP
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = 9987 )
         IF( NOUT.NE.6 .AND. NOUT.NE.0 )
     $      CLOSE ( NOUT )
      END IF
*
      CALL BLACS_EXIT( 0 )
*
 9999 FORMAT( 'ILLEGAL ', A6, ': ', A5, ' = ', I3,
     $      '; It should be at least 1' )
 9998 FORMAT( 'ILLEGAL GRID: nprow*npcol = ', I4, '. It can be at most',
     $      I4 )
 9997 FORMAT( 'Bad ', A6, ' parameters: going on to next test case.' )
 9996 FORMAT( 'Unable to perform ', A, ': need TOTMEM of at least',
     $      I11 )
 9995 FORMAT( 'TIME     N  NB    P    Q NEP Time   MFLOPS  CHECK' )
 9994 FORMAT( '---- ----- --- ---- ---- -------- -------- ------' )
 9993 FORMAT( A4, 1X, I5, 1X, I3, 1X, I4, 1X, I4, 1X, F8.2, 1X, F8.2,
     $      1X, A6 )
 9992 FORMAT( 'Finished ', I6, ' tests, with the following results:' )
 9991 FORMAT( I5, ' tests completed and passed residual checks.' )
 9990 FORMAT( I5, ' tests completed without checking.' )
 9989 FORMAT( I5, ' tests completed and failed residual checks.' )
 9988 FORMAT( I5, ' tests skipped because of illegal input values.' )
 9987 FORMAT( 'END OF TESTS.' )
 9986 FORMAT( '||H - Q*S*Q^T|| / (||H|| * N * eps) = ', G25.7 )
 9985 FORMAT( '||Q^T*Q - I|| / ( N * eps ) ', G25.7 )
*
      STOP
*
*     End of PCNEPDRIVER
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
	SUBROUTINE SNTBFE ( LUNT, IFXYN, LINE )

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    SNTBFE
C   PRGMMR: ATOR            ORG: NCEP       DATE: 2017-11-02
C
C ABSTRACT:  THIS SUBROUTINE READS A COMPLETE ENTRY (CORRESPONDING
C   TO THE INPUT FXY NUMBER) FROM AN ASCII MASTER CODE/FLAG TABLE AND
C   STORES THE OUTPUT INTO AN INTERNAL MEMORY STRUCTURE.
C
C PROGRAM HISTORY LOG:
C 2017-11-02  J. ATOR    -- ORIGINAL AUTHOR
C
C USAGE:    CALL SNTBFE ( LUNT, IFXYN, LINE )
C   INPUT ARGUMENT LIST:
C     LUNT     - INTEGER: FORTRAN LOGICAL UNIT NUMBER OF ASCII FILE
C                CONTAINING MASTER CODE/FLAG TABLE INFORMATION
C     IFXYN    - INTEGER: BIT-WISE REPRESENTATION OF FXY NUMBER
C     LINE     - CHARACTER*(*): FIRST LINE OF TABLE ENTRY
C
C REMARKS:
C    THIS ROUTINE CALLS:        ADN30    BORT2    IFXY     IGETFXY
C                               IGETNTBL JSTCHR   PARSTR   STRNUM
C                               STRTBFE
C    THIS ROUTINE IS CALLED BY: RDMTBF
C                               Normally not called by any application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

	CHARACTER*(*)	LINE
	CHARACTER*160	CLINE, TAGS(4), CDSTR(2), ADSC(10), CVAL(25)
	CHARACTER*128	BORT_STR1, BORT_STR2
	CHARACTER*6	ADN30, CLEMON, CDSC
	DIMENSION	IDFXY(10), IDVAL(25)

	LOGICAL	DONE, LSTNBLK

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

C	We already have the FXY number.  Now we need to read and parse
C	all of the remaining lines from the table entry for this FXY
C	number.  The information for each individual code figure or bit
C	number will then be stored as a separate entry within the
C	internal memory structure.

	DONE = .FALSE.
	NIDFXY = 0
	NIDVAL = 0

	DO WHILE ( .NOT. DONE )

	    IF ( IGETNTBL ( LUNT, CLINE ) .NE. 0 ) THEN
		BORT_STR2 = '                  IS INCOMPLETE'
		GOTO 900
	    ENDIF

	    CALL PARSTR ( CLINE, TAGS, 4, NTAG, '|', .FALSE. )
	    IF ( ( NTAG .LT. 2 ) .OR. ( NTAG .GT. 3 ) ) THEN
		BORT_STR2 = '                  HAS BAD CARD'
		GOTO 900
	    ENDIF

	    IF ( NTAG .EQ. 2 ) THEN

C		This line contains a list of dependencies.

		CALL PARSTR ( TAGS(2), CDSTR, 2, NTAG, '=', .FALSE. )
		IF ( NTAG .NE. 2 ) THEN
		  BORT_STR2 = '           HAS BAD DEPENDENCY CARD'
		  GOTO 900
		ENDIF

C		Parse the list of FXY numbers.

		CALL PARSTR ( CDSTR(1), ADSC, 10, NIDFXY, ',', .FALSE. )
		IF ( NIDFXY .EQ. 0 ) THEN
		  BORT_STR2 = '        HAS BAD DEPENDENCY LIST (FXY)'
		  GOTO 900
		ENDIF
		DO II = 1, NIDFXY
		  IF ( IGETFXY ( ADSC(II), CDSC ) .NE. 0 ) THEN
		    BORT_STR2 = '        HAS BAD DEPENDENCY (FXY)'
		    GOTO 900
		  ENDIF
		  IDFXY(II) = IFXY( CDSC )
		ENDDO

C		Parse the list of values.

		CALL PARSTR ( CDSTR(2), CVAL, 25, NIDVAL, ',', .FALSE. )
		IF ( NIDVAL .EQ. 0 ) THEN
		  BORT_STR2 = '        HAS BAD DEPENDENCY LIST (VAL)'
		  GOTO 900
		ENDIF
		DO II = 1, NIDVAL
		  CALL JSTCHR ( CVAL(II), IER )
		  CALL STRNUM ( CVAL(II), IVAL )
		  IDVAL(II) = IVAL
		ENDDO

	    ELSE

C		This line contains a value (code figure or bit number)
C		and corresponding meaning.

		IPT = INDEX ( TAGS(2), ' >' )
		IF ( IPT .EQ. 0 ) THEN

C		  This is the last line for this table entry.

		  DONE = .TRUE.
		ELSE
		  TAGS(2)(IPT+1:IPT+1) = ' '
		ENDIF

		CALL JSTCHR ( TAGS(2), IER )
		CALL STRNUM ( TAGS(2), IVAL )

C		Find the last non-blank character in the meaning string.

		CALL JSTCHR ( TAGS(3), IER )
		LT3 = LEN(TAGS(3))
		LSTNBLK = .FALSE.
		DO WHILE ( ( LT3 .GT. 0 ) .AND. ( .NOT. LSTNBLK ) )
		  IF ( TAGS(3)(LT3:LT3) .NE. ' ' ) THEN
		    LSTNBLK = .TRUE.
		  ELSE
		    LT3 = LT3 - 1
		  ENDIF
		ENDDO

C		Store the information for this value within the internal
C		memory structure.

		IF ( ( NIDFXY .EQ. 0 ) .AND. ( NIDVAL .EQ. 0 ) ) THEN
		  CALL STRTBFE ( IFXYN, IVAL, TAGS(3), LT3, -1, -1 )
		ELSE
		  DO II = 1, NIDFXY
		    DO JJ = 1, NIDVAL
		      CALL STRTBFE ( IFXYN, IVAL, TAGS(3), LT3,
     +				     IDFXY(II), IDVAL(JJ) )
		    ENDDO
		  ENDDO
		ENDIF

	    ENDIF

	ENDDO

	RETURN

 900	CLEMON = ADN30 ( IFXYN, 6 )
	WRITE(BORT_STR1,'("BUFRLIB: SNTBFE - TABLE F ENTRY FOR' //
     .     ' ELEMENT DESCRIPTOR: ",5A)')
     .     CLEMON(1:1), '-', CLEMON(2:3), '-', CLEMON(4:6)
	CALL BORT2(BORT_STR1,BORT_STR2)
	END

ext: f


text:
! This module is used both for binned and peeloff images

module type_image

  use core_lib
  use mpi_hdf5_io
  use type_photon

  implicit none
  save

  private
  public :: image
  public :: image_scale
  public :: image_bin
  public :: image_bin_raytraced
  public :: image_setup
  public :: image_write
  public :: in_image

  public :: image_raytracing_initialize
  public :: image_raytracing_set_spectrum
  public :: image_raytracing_set_blackbody
  public :: image_raytracing_set_emissivity
  public :: image_raytracing_set_opacity

  public :: source_helper
  public :: dust_helper

  type source_helper
     real(dp), allocatable :: spectrum(:)
  end type source_helper

  type dust_helper
     real(dp), allocatable :: chi(:)
     real(dp), allocatable :: emissivity_var(:)
     real(dp), allocatable :: log10_emissivity(:,:)
  end type dust_helper

  type image

     ! IMAGE PARAMETERS

     integer :: n_x,n_y
     real(dp) :: x_min,x_max,y_min,y_max
     ! number of pixels in each direction, and view limits

     integer :: n_nu
     logical :: use_exact_nu
     integer :: inu_min,inu_max
     real(dp) :: nu_min,nu_max
     real(dp) :: log10_nu_min,log10_nu_max
     real(dp),allocatable :: nu(:)
     ! number of frequencies, and frequency range

     integer :: n_ap
     real(dp) :: ap_min,ap_max
     real(dp) :: log10_ap_min,log10_ap_max
     ! number of apertures for the SED, and aperture range

     integer :: n_view
     type(angle3d_dp), allocatable :: a(:)
     ! number of viewing angles, and viewing angles

     ! IMAGE AND SED

     logical :: compute_image
     type(stokes_dp),allocatable :: img(:,:,:,:,:)
     type(stokes_dp),allocatable :: img2(:,:,:,:,:)
     real(dp),allocatable :: imgn(:,:,:,:,:)
     ! the stokes Image array and the summation array for the variance

     logical :: compute_sed
     type(stokes_dp),allocatable :: sed(:,:,:,:)
     type(stokes_dp),allocatable :: sed2(:,:,:,:)
     real(dp),allocatable :: sedn(:,:,:,:)
     ! the stokes SED array and the summation array for the variance

     ! RAYTRACING

     type(source_helper),allocatable :: sources(:)
     ! spectra of the sources, binned to the image spectral resolution

     type(dust_helper),allocatable :: dust(:)
     ! emissivity of the dust types, binned to the image spectral resolution

     real(dp),allocatable :: tmp_spectrum(:)

     character(len=10) :: track_origin
     logical :: uncertainties

     integer :: io_type

     integer :: n_orig

     integer :: n_sources
     integer :: n_dust

  end type image

  interface image_raytracing_set_spectrum
     module procedure image_raytracing_set_spectrum_pdf
     module procedure image_raytracing_set_spectrum_nu_fnu
  end interface image_raytracing_set_spectrum

contains

  integer function orig(p)
    ! Find the origin flag for a photon
    implicit none
    type(photon),intent(in) :: p
    if(p%scattered) then
       if(p%reprocessed) then
          orig = 4 ! scattered dust emission
       else
          orig = 3 ! scattered source emission
       end if
    else
       if(p%reprocessed) then
          orig = 2 ! dust emission
       else
          orig = 1 ! source emission
       end if
    end if
  end function orig

  subroutine image_scale(img,scale)

    implicit none

    type(image), intent(inout) :: img
    real(dp), intent(in) :: scale

    if(img%compute_image) then
       img%img%i = img%img%i * scale
       img%img%q = img%img%q * scale
       img%img%u = img%img%u * scale
       img%img%v = img%img%v * scale
    end if
    if(img%compute_sed) then
       img%sed%i = img%sed%i * scale
       img%sed%q = img%sed%q * scale
       img%sed%u = img%sed%u * scale
       img%sed%v = img%sed%v * scale
    end if

    if(img%uncertainties) then
       if(img%compute_image) then
          img%img2%i = img%img2%i * scale**2
          img%img2%q = img%img2%q * scale**2
          img%img2%u = img%img2%u * scale**2
          img%img2%v = img%img2%v * scale**2
       end if
       if(img%compute_sed) then
          img%sed2%i = img%sed2%i * scale**2
          img%sed2%q = img%sed2%q * scale**2
          img%sed2%u = img%sed2%u * scale**2
          img%sed2%v = img%sed2%v * scale**2
       end if
    end if

  end subroutine image_scale

  subroutine image_setup(handle,path,img,n_view,n_sources,n_dust,use_exact_nu,frequencies)

    implicit none
    integer(hid_t),intent(in) :: handle
    character(len=*),intent(in) :: path

    type(image),intent(out) :: img
    integer,intent(in) :: n_view, n_sources, n_dust
    logical,intent(in) :: use_exact_nu
    real(dp),intent(in),optional :: frequencies(:)
    real(dp) :: wav_min, wav_max
    integer :: io_bytes

    img%n_view = n_view
    img%n_sources = n_sources
    img%n_dust = n_dust

    call mp_read_keyword(handle, path, 'n_wav',img%n_nu)

    call mp_read_keyword(handle, path, 'compute_image',img%compute_image)
    if(img%compute_image) then
       call mp_read_keyword(handle, path, 'n_x',img%n_x)
       call mp_read_keyword(handle, path, 'n_y',img%n_y)
       call mp_read_keyword(handle, path, 'x_min',img%x_min)
       call mp_read_keyword(handle, path, 'x_max',img%x_max)
       call mp_read_keyword(handle, path, 'y_min',img%y_min)
       call mp_read_keyword(handle, path, 'y_max',img%y_max)
    end if

    call mp_read_keyword(handle, path, 'compute_sed',img%compute_sed)
    if(img%compute_sed) then
       call mp_read_keyword(handle, path, 'n_ap',img%n_ap)
       call mp_read_keyword(handle, path, 'ap_min',img%ap_min)
       call mp_read_keyword(handle, path, 'ap_max',img%ap_max)
       img%log10_ap_min = log10(img%ap_min)
       img%log10_ap_max = log10(img%ap_max)
    end if

    call mp_read_keyword(handle, path, 'track_origin',img%track_origin)

    select case(trim(img%track_origin))
    case('detailed')
       img%n_orig = 2 * (img%n_sources + img%n_dust)
    case('basic', 'yes')
       img%n_orig = 4
    case('no')
       img%n_orig = 1
    case default
       call error("image_setup", "unknown track_origin flag: "//trim(img%track_origin))
    end select

    call mp_read_keyword(handle, path, 'uncertainties',img%uncertainties)

    if(use_exact_nu) then

       img%use_exact_nu = .true.

       call mp_read_keyword(handle, path, 'inu_min',img%inu_min)
       call mp_read_keyword(handle, path, 'inu_max',img%inu_max)

       if(.not.present(frequencies)) call error('image_setup','frequencies should be given if use_exact_nu is .true.')

       allocate(img%nu(size(frequencies)))
       img%nu = frequencies(img%inu_min:img%inu_max)

    else

       img%use_exact_nu = .false.

       call mp_read_keyword(handle, path, 'wav_min',wav_min)
       call mp_read_keyword(handle, path, 'wav_max',wav_max)

       img%nu_min = c_cgs / (wav_max * 1.e-4)
       img%nu_max = c_cgs / (wav_min * 1.e-4)

       img%log10_nu_min = log10(img%nu_min)
       img%log10_nu_max = log10(img%nu_max)

    end if

    ! Allocate arrays

    if(img%compute_image) then
       allocate(img%img(img%n_nu,img%n_x,img%n_y,img%n_view,img%n_orig))
       if(img%uncertainties) then
          allocate(img%img2(img%n_nu,img%n_x,img%n_y,img%n_view,img%n_orig))
          allocate(img%imgn(img%n_nu,img%n_x,img%n_y,img%n_view,img%n_orig))
       end if
    end if

    if(img%compute_sed) then
       allocate(img%sed(img%n_nu,img%n_ap,img%n_view,img%n_orig))
       if(img%uncertainties) then
          allocate(img%sed2(img%n_nu,img%n_ap,img%n_view,img%n_orig))
          allocate(img%sedn(img%n_nu,img%n_ap,img%n_view,img%n_orig))
       end if
    end if

    ! Initialize all arrays

    if(img%compute_image) then
       img%img = stokes_dp(0.,0.,0.,0.)
       if(img%uncertainties) then
          img%img2 = stokes_dp(0.,0.,0.,0.)
          img%imgn = 0._dp
       end if
    end if

    if(img%compute_sed) then
       img%sed = stokes_dp(0.,0.,0.,0.)
       if(img%uncertainties) then
          img%sed2 = stokes_dp(0.,0.,0.,0.)
          img%sedn = 0._dp
       end if
    end if

    call mp_read_keyword(handle, path, 'io_bytes',io_bytes)

    select case(io_bytes)
    case(4)
       img%io_type = sp
    case(8)
       img%io_type = dp
    case default
       call error("image_setup", "unexpected value of io_bytes (should be 4 or 8)")
    end select

  end subroutine image_setup

  subroutine find_sed_bin(img,x_image,y_image,ir)

    implicit none

    type(image),intent(in) :: img  ! Image
    real(dp),intent(in) :: x_image, y_image
    integer,intent(out) :: ir
    real(dp) :: log10_r_image

    log10_r_image = log10(sqrt(x_image*x_image + y_image*y_image))

    if(log10_r_image < img%log10_ap_min .or. img%n_ap == 1) then
       ir = 1
    else
       ir = ipos(img%log10_ap_min,img%log10_ap_max,log10_r_image,img%n_ap-1) + 1
    end if

  end subroutine find_sed_bin

  subroutine find_image_bin(img,x_image,y_image,ix,iy)

    implicit none

    type(image),intent(in)    :: img  ! Image
    real(dp),intent(in) :: x_image, y_image
    integer,intent(out) :: ix,iy

    ix = ipos(img%x_min,img%x_max,x_image,img%n_x)
    iy = ipos(img%y_min,img%y_max,y_image,img%n_y)

  end subroutine find_image_bin

  logical function in_image(img, x_image, y_image, nu)

    implicit none

    type(image),intent(in) :: img
    real(dp),intent(in) :: x_image, y_image
    real(dp),intent(in),optional :: nu

    in_image = .false.

    ! Check if photon would end up inside the frequency range. Note that we
    ! can only know for sure that the photon is outside the range in which
    ! case in_image is already set to .false. so we just return. If it's
    ! inside the frequency range, then we still have to check if it would
    ! fall in the image or SEDs.
    if(present(nu)) then
       if(nu < img%nu_min .or. nu > img%nu_max) return
    end if

    ! Check if photon would fall in image
    if(img%compute_image) then
       if((x_image >= img%x_min .and. x_image <= img%x_max) .or. (x_image <= img%x_min .and. x_image >= img%x_max)) then
          if((y_image >= img%y_min .and. y_image <= img%y_max) .or. (y_image <= img%y_min .and. y_image >= img%y_max)) then
             in_image = .true.
             return
          end if
       end if
    end if

    ! Check if photon would fall in SED
    if(img%compute_sed) then
       if(x_image*x_image + y_image*y_image <= img%ap_max * img%ap_max) then
          in_image = .true.
          return
       end if
    end if

  end function in_image

  subroutine image_bin(img,p,x_image,y_image,im)

    implicit none

    type(image),intent(inout)    :: img  ! Image
    type(photon),intent(in) :: p
    real(dp),intent(in) :: x_image, y_image
    integer,intent(in) :: im ! sub-image to bin into
    real(dp) :: log10_nu_image
    integer :: ix,iy,ir,inu,io ! Bins

    if(img%compute_image.and..not.allocated(img%img)) call error('bin_photon','Image not allocated')
    if(img%compute_sed.and..not.allocated(img%sed)) call error('bin_photon','SED not allocated')

    if(is_nan(p%energy)) then
       call warn("image_bin","photon has NaN energy - ignoring")
       return
    end if

    if(is_nan(p%s%i)) then
       call warn("image_bin","photon has NaN flux - ignoring")
       return
    end if

    ! Find frequency bin
    if(img%use_exact_nu) then
       inu = p%inu - img%inu_min + 1
    else
       log10_nu_image = log10(p%nu)
       inu = ipos(img%log10_nu_min,img%log10_nu_max,log10_nu_image,img%n_nu)
    end if

    ! Find origin flag
    if(trim(img%track_origin) == 'detailed') then
       io = ((orig(p) - mod(orig(p),2)) * img%n_sources + (orig(p) - mod(orig(p)+1,2) - 1) * img%n_dust) / 2
       if(mod(orig(p),2)==0) then
          io = io + p%dust_id
       else
          io = io + p%source_id
       end if
    else if(trim(img%track_origin) == 'basic') then
       io = orig(p)
    else
       io = 1
    end if

    if(inu >= 1 .and. inu <= img%n_nu) then
       if(img%compute_image) then
          call find_image_bin(img,x_image,y_image,ix,iy)
          if(ix >= 1 .and. ix <= img%n_x) then
             if(iy >= 1 .and. iy <= img%n_y) then
                img%img(inu,ix,iy,im,io) = img%img(inu,ix,iy,im,io) + p%s * p%energy
                if(img%uncertainties) then
                   img%img2(inu,ix,iy,im,io) = img%img2(inu,ix,iy,im,io) + p%s**2._dp * p%energy**2._dp
                   img%imgn(inu,ix,iy,im,io) = img%imgn(inu,ix,iy,im,io) + 1._dp
                end if
             end if
          end if
       end if
       if(img%compute_sed) then
          call find_sed_bin(img,x_image,y_image,ir)
          if(ir >= 1 .and. ir <= img%n_ap) then
             img%sed(inu,ir,im,io) = img%sed(inu,ir,im,io) + p%s * p%energy
             if(img%uncertainties) then
                img%sed2(inu,ir,im,io) = img%sed2(inu,ir,im,io) + p%s**2._dp * p%energy**2._dp
                img%sedn(inu,ir,im,io) = img%sedn(inu,ir,im,io) + 1._dp
             end if
          end if
       end if
    end if

  end subroutine image_bin

  subroutine image_bin_raytraced(img,p,x_image,y_image,im,column_density)

    implicit none

    type(image),intent(inout)    :: img  ! Image
    type(photon),intent(in) :: p
    real(dp),intent(in) :: x_image, y_image
    integer,intent(in) :: im ! sub-image to bin into
    real(dp),intent(in) :: column_density(:)
    integer :: ix,iy,ir,iw,id,io ! Bins

    if(img%compute_image.and..not.allocated(img%img)) call error('bin_photon','Image not allocated')
    if(img%compute_sed.and..not.allocated(img%sed)) call error('bin_photon','SED not allocated')

    if(is_nan(p%energy)) then
       call warn("image_bin_raytraced","photon has NaN energy - ignoring")
       return
    end if

    if(is_nan(p%s%i)) then
       call warn("image_bin_raytraced","photon has NaN flux - ignoring")
       return
    end if

    ! Find origin flag
    if(trim(img%track_origin) == 'detailed') then
       io = ((orig(p) - mod(orig(p),2)) * img%n_sources + (orig(p) - mod(orig(p)+1,2) - 1) * img%n_dust) / 2
       if(mod(orig(p),2)==0) then
          io = io + p%dust_id
       else
          io = io + p%source_id
       end if
    else if(trim(img%track_origin) == 'basic') then
       io = orig(p)
    else
       io = 1
    end if

    ! Now the fun begins
    select case(p%emiss_type)
    case(1,2)
       img%tmp_spectrum = img%sources(p%source_id)%spectrum
    case(3)
       img%tmp_spectrum = (img%dust(p%dust_id)%log10_emissivity(:,p%emiss_var_id+1) - &
            &                   img%dust(p%dust_id)%log10_emissivity(:,p%emiss_var_id)) * p%emiss_var_frac + &
            &                   img%dust(p%dust_id)%log10_emissivity(:,p%emiss_var_id)
       img%tmp_spectrum = 10._dp**(img%tmp_spectrum)
       where(is_nan(img%tmp_spectrum))
          img%tmp_spectrum = 0.
       end where
    case default
       stop "unknown emiss_type"
    end select

    img%tmp_spectrum = img%tmp_spectrum * p%s%i * p%energy

    do id=1,size(column_density)
       img%tmp_spectrum = img%tmp_spectrum * exp(- column_density(id) * img%dust(id)%chi)
    end do

    if(img%compute_image) then
       call find_image_bin(img,x_image,y_image,ix,iy)
       if(ix >= 1 .and. ix <= img%n_x) then
          if(iy >= 1 .and. iy <= img%n_y) then
             do iw=1,img%n_nu
                img%img(iw,ix,iy,im,io)%i = img%img(iw,ix,iy,im,io)%i + img%tmp_spectrum(iw)
                if(img%uncertainties) then
                   img%img2(iw,ix,iy,im,io)%i = img%img2(iw,ix,iy,im,io)%i + img%tmp_spectrum(iw)**2._dp
                   img%imgn(iw,ix,iy,im,io) = img%imgn(iw,ix,iy,im,io) + 1._dp
                end if
             end do
          end if
       end if
    end if

    if(img%compute_sed) then
       call find_sed_bin(img,x_image,y_image,ir)
       if(ir >= 1 .and. ir <= img%n_ap) then
          do iw=1,img%n_nu
             img%sed(iw,ir,im,io)%i = img%sed(iw,ir,im,io)%i + img%tmp_spectrum(iw)
             if(img%uncertainties) then
                img%sed2(iw,ir,im,io)%i = img%sed2(iw,ir,im,io)%i + img%tmp_spectrum(iw)**2._dp
                img%sedn(iw,ir,im,io) = img%sedn(iw,ir,im,io) + 1._dp
             end if
          end do
       end if
    end if

  end subroutine image_bin_raytraced

  subroutine image_write(img,group)

    implicit none

    type(image),intent(in)      :: img  ! Image to write out
    integer(hid_t),intent(in) :: group

    real(dp),allocatable :: cube5d(:,:,:,:,:)
    real(dp),allocatable :: cube5de(:,:,:,:,:)
    real(dp),allocatable :: cube6d(:,:,:,:,:,:)
    real(dp),allocatable :: cube6de(:,:,:,:,:,:)

    integer  :: ia,inu

    real(dp) :: dnunorm

    ! The factor for SED normalization is derived in the following way:
    !
    ! Given n_nu equal bins in log10 space from nu_min to nu_max, the center
    ! of the bins is given by:
    !
    ! log10[nu(j)] = [log10(nu_max) - log10(nu_min)] * (j - 0.5) / n_nu
    !
    ! Thus:
    !
    ! nu(j) = (nu_max / nu_min) ** ((j - 0.5) / n_nu)
    !
    ! Therefore, the width of bins is given by
    !
    ! dnu(j) = nu(j + 0.5) - nu(j - 0.5)
    ! dnu(j) = (nu_max / nu_min) ** (j / n_nu)
    !        - (nu_max / nu_min) ** ((j - 1) / n_nu)
    !
    ! The SED normalization is done by converting F_nu * dnu to nu * Fnu, i.e.
    ! dividing by dnu(j) and multiplying by nu(j), i.e. dividing by dnunorm
    ! where:
    !
    ! dnunorm = dnu(j) / nu(j)
    ! dnunorm = (nu_max / nu_min) ** (+0.5 / n_nu)
    !         - (nu_max / nu_min) ** (-0.5 / n_nu)

    dnunorm = (img%nu_max / img%nu_min) ** (+0.5_dp / real(img%n_nu, dp)) &
         - (img%nu_max / img%nu_min) ** (-0.5_dp / real(img%n_nu, dp))

    if(img%compute_sed) then

       write(*,'(" [image_write] writing out SEDs")')

       allocate(cube5d(img%n_nu, img%n_ap, img%n_view,img%n_orig,4))
       if(img%uncertainties) allocate(cube5de(img%n_nu, img%n_ap, img%n_view,img%n_orig,4))

       cube5d(:,:,:,:,1) = img%sed%i
       cube5d(:,:,:,:,2) = img%sed%q
       cube5d(:,:,:,:,3) = img%sed%u
       cube5d(:,:,:,:,4) = img%sed%v

       if(img%uncertainties) then
          where(img%sedn > 1)
             cube5de(:,:,:,:,1) = sqrt((img%sed2%i + (img%sed%i)**2 / img%sedn) / (img%sedn - 1)) * sqrt(img%sedn)
             cube5de(:,:,:,:,2) = sqrt((img%sed2%q + (img%sed%q)**2 / img%sedn) / (img%sedn - 1)) * sqrt(img%sedn)
             cube5de(:,:,:,:,3) = sqrt((img%sed2%u + (img%sed%u)**2 / img%sedn) / (img%sedn - 1)) * sqrt(img%sedn)
             cube5de(:,:,:,:,4) = sqrt((img%sed2%v + (img%sed%v)**2 / img%sedn) / (img%sedn - 1)) * sqrt(img%sedn)
          elsewhere
             cube5de(:,:,:,:,1) = 0._dp
             cube5de(:,:,:,:,2) = 0._dp
             cube5de(:,:,:,:,3) = 0._dp
             cube5de(:,:,:,:,4) = 0._dp
          end where
       end if

       if(.not.img%use_exact_nu) then
          cube5d = cube5d / dnunorm
          if(img%uncertainties) cube5de = cube5de / dnunorm
       else
          do inu=1,img%n_nu
             cube5d(inu,:,:,:,:) = cube5d(inu,:,:,:,:) * img%nu(inu)
             if(img%uncertainties) cube5de(inu,:,:,:,:) = cube5de(inu,:,:,:,:) * img%nu(inu)
          end do
       end if

       do ia=2,img%n_ap
          cube5d(:,ia,:,:,:) = cube5d(:,ia-1,:,:,:) + cube5d(:,ia,:,:,:)
          if(img%uncertainties) cube5de(:,ia,:,:,:) = sqrt(cube5de(:,ia-1,:,:,:)**2 + cube5de(:,ia,:,:,:)**2)
       end do

       select case(img%io_type)
       case(sp)
          call mp_write_array(group, 'seds', real(cube5d, sp))
          if(img%uncertainties) call mp_write_array(group, 'seds_unc', real(cube5de, sp))
       case(dp)
          call mp_write_array(group, 'seds', real(cube5d, dp))
          if(img%uncertainties) call mp_write_array(group, 'seds_unc', real(cube5de, dp))
       case default
          call error("image_write","unexpected value of img%io_type (should be sp or dp)")
       end select

       if(.not.img%use_exact_nu) then
          call mp_write_keyword(group, 'seds','numin',img%nu_min)
          call mp_write_keyword(group, 'seds','numax',img%nu_max)
       end if
       call mp_write_keyword(group, 'seds','apmin',img%ap_min)
       call mp_write_keyword(group, 'seds','apmax',img%ap_max)

       call mp_write_keyword(group, 'seds', 'track_origin', img%track_origin)
       if(trim(img%track_origin) == 'detailed') then
          call mp_write_keyword(group, 'seds', 'n_sources', img%n_sources)
          call mp_write_keyword(group, 'seds', 'n_dust', img%n_dust)
       end if

    end if

    if(img%compute_image) then

       write(*,'(" [image_write] writing out images")')

       allocate(cube6d(img%n_nu, img%n_x, img%n_y, img%n_view, img%n_orig, 4))
       if(img%uncertainties) allocate(cube6de(img%n_nu, img%n_x, img%n_y, img%n_view, img%n_orig, 4))

       cube6d(:,:,:,:,:,1) = img%img%i
       cube6d(:,:,:,:,:,2) = img%img%q
       cube6d(:,:,:,:,:,3) = img%img%u
       cube6d(:,:,:,:,:,4) = img%img%v

       if(img%uncertainties) then
          where(img%imgn > 1)
             cube6de(:,:,:,:,:,1) = sqrt((img%img2%i + (img%img%i)**2 / img%imgn) / (img%imgn - 1)) * sqrt(img%imgn)
             cube6de(:,:,:,:,:,2) = sqrt((img%img2%q + (img%img%q)**2 / img%imgn) / (img%imgn - 1)) * sqrt(img%imgn)
             cube6de(:,:,:,:,:,3) = sqrt((img%img2%u + (img%img%u)**2 / img%imgn) / (img%imgn - 1)) * sqrt(img%imgn)
             cube6de(:,:,:,:,:,4) = sqrt((img%img2%v + (img%img%v)**2 / img%imgn) / (img%imgn - 1)) * sqrt(img%imgn)
          elsewhere
             cube6de(:,:,:,:,:,1) = 0._dp
             cube6de(:,:,:,:,:,2) = 0._dp
             cube6de(:,:,:,:,:,3) = 0._dp
             cube6de(:,:,:,:,:,4) = 0._dp
          end where
       end if

       if(.not.img%use_exact_nu) then
          cube6d = cube6d / dnunorm
          if(img%uncertainties) cube6de = cube6de / dnunorm
       else
          do inu=1,img%n_nu
             cube6d(inu,:,:,:,:,:) = cube6d(inu,:,:,:,:,:) * img%nu(inu)
             if(img%uncertainties) cube6de(inu,:,:,:,:,:) = cube6de(inu,:,:,:,:,:) * img%nu(inu)
          end do
       end if

       select case(img%io_type)
       case(sp)
          call mp_write_array(group, 'images', real(cube6d, sp))
          if(img%uncertainties) call mp_write_array(group, 'images_unc', real(cube6de, sp))
       case(dp)
          call mp_write_array(group, 'images', real(cube6d, dp))
          if(img%uncertainties) call mp_write_array(group, 'images_unc', real(cube6de, dp))
       case default
          call error("image_write","unexpected value of img%io_type (should be sp or dp)")
       end select

       if(.not.img%use_exact_nu) then
          call mp_write_keyword(group, 'images','numin',img%nu_min)
          call mp_write_keyword(group, 'images','numax',img%nu_max)
       end if
       call mp_write_keyword(group, 'images','xmin',img%x_min)
       call mp_write_keyword(group, 'images','xmax',img%x_max)
       call mp_write_keyword(group, 'images','ymin',img%y_min)
       call mp_write_keyword(group, 'images','ymax',img%y_max)

       call mp_write_keyword(group, 'images', 'track_origin', img%track_origin)
       if(trim(img%track_origin) == 'detailed') then
          call mp_write_keyword(group, 'images', 'n_sources', img%n_sources)
          call mp_write_keyword(group, 'images', 'n_dust', img%n_dust)
       end if

    end if

    if(img%use_exact_nu) then
       call mp_table_write_header(group, 'frequencies',img%inu_max - img%inu_min + 1,1,(/'nu'/),(/1/),(/h5t_ieee_f64le/))
       call mp_table_write_column(group, 'frequencies','nu',img%nu(img%inu_min:img%inu_max))
    end if

    write(*,'(" [image_write] done")')

  end subroutine image_write

  ! RAYTRACING

  subroutine image_raytracing_initialize(img,n_sources,n_dust)
    implicit none
    type(image), intent(inout) :: img
    integer,intent(in) :: n_sources, n_dust
    allocate(img%sources(n_sources))
    allocate(img%dust(n_dust))
    allocate(img%tmp_spectrum(img%n_nu))
  end subroutine image_raytracing_initialize

  subroutine image_raytracing_set_spectrum_pdf(img,source_id,spectrum)

    implicit none

    type(image), intent(inout) :: img
    integer,intent(in) :: source_id
    type(pdf_dp),intent(in) :: spectrum

    call image_raytracing_set_spectrum(img,source_id,spectrum%x,spectrum%pdf)

  end subroutine image_raytracing_set_spectrum_pdf

  subroutine image_raytracing_set_spectrum_nu_fnu(img,source_id,nu,fnu)

    implicit none

    type(image), intent(inout) :: img
    integer,intent(in) :: source_id
    real(dp),intent(in) :: nu(:), fnu(:)
    real(dp), dimension(img%n_nu) :: nu_new, e_new
    real(dp) :: numin, numax
    integer :: inu

    allocate(img%sources(source_id)%spectrum(img%n_nu))

    if(img%use_exact_nu) then

       img%sources(source_id)%spectrum = interp1d_loglog(nu, fnu, img%nu, bounds_error=.false., fill_value=0._dp)

    else

       do inu=1,img%n_nu

          numin = 10._dp**(img%log10_nu_min + (img%log10_nu_max - img%log10_nu_min) * real(inu-1, dp) / real(img%n_nu, dp))
          numax = 10._dp**(img%log10_nu_min + (img%log10_nu_max - img%log10_nu_min) * real(inu, dp) / real(img%n_nu, dp))

          nu_new(inu) = 10._dp**(0.5_dp * (log10(numin) + log10(numax)))
          e_new(inu) = integral_loglog(nu, fnu, numin, numax)

       end do

       e_new = e_new / integral_loglog(nu, fnu)

       img%sources(source_id)%spectrum = e_new

    end if

  end subroutine image_raytracing_set_spectrum_nu_fnu

  subroutine image_raytracing_set_blackbody(img,source_id,temperature)

    implicit none

    type(image), intent(inout) :: img
    integer,intent(in) :: source_id
    real(dp),intent(in) :: temperature
    integer :: inu

    integer :: n_nu
    real(dp),allocatable :: nu(:)

    real(dp) :: log10_nu_min, log10_nu_max

    if(img%use_exact_nu) then

       call image_raytracing_set_spectrum(img,source_id,img%nu,normalized_B_nu(img%nu,temperature))

    else

       log10_nu_min = log10(3.e9_dp)  ! 10 cm (0.05K)
       log10_nu_max = log10(3.e16_dp) ! 10 nm (1/2 million K)

       n_nu = ceiling((log10_nu_max - log10_nu_min) * 100000)
       allocate(nu(n_nu))
       do inu=1,n_nu
          nu(inu) = 10._dp**(real(inu-1,dp)/real(n_nu-1,dp)*(log10_nu_max - log10_nu_min) + log10_nu_min)
       end do

       call image_raytracing_set_spectrum(img,source_id,nu,normalized_B_nu(nu,temperature))

    end if

  contains

    elemental real(dp) function normalized_B_nu(nu,T)
      implicit none
      real(dp),intent(in) :: nu,T
      real(dp),parameter :: a = two * h_cgs / c_cgs / c_cgs / stef_boltz * pi
      real(dp),parameter :: b = h_cgs / k_cgs
      real(dp) :: T4
      T4 = T*T*T*T
      normalized_B_nu = a * nu * nu * nu / ( exp(b*nu/T) - one) / T4
    end function normalized_B_nu

  end subroutine image_raytracing_set_blackbody

  subroutine image_raytracing_set_emissivity(img,dust_id,emissivity_var,emissivity)

    implicit none

    type(image),intent(inout) :: img
    integer,intent(in) :: dust_id
    real(dp),intent(in) :: emissivity_var(:)
    type(pdf_dp),intent(in) :: emissivity(:)
    real(dp), dimension(img%n_nu) :: nu_new
    real(dp), dimension(img%n_nu, size(emissivity_var)) :: e_new
    real(dp) :: numin, numax
    integer :: inu
    integer :: j

    if(img%use_exact_nu) then

       do j=1,size(emissivity_var)
          e_new(:,j) = interp1d_loglog(emissivity(j)%x, emissivity(j)%pdf, img%nu, bounds_error=.false., fill_value=0._dp)
       end do

    else

       do j=1,size(emissivity_var)

          do inu=1,img%n_nu

             numin = 10._dp**(img%log10_nu_min + (img%log10_nu_max - img%log10_nu_min) * real(inu-1, dp) / real(img%n_nu, dp))
             numax = 10._dp**(img%log10_nu_min + (img%log10_nu_max - img%log10_nu_min) * real(inu, dp) / real(img%n_nu, dp))

             nu_new(inu) = 10._dp**(0.5_dp * (log10(numin) + log10(numax)))
             e_new(inu,j) = integral_loglog(emissivity(j)%x, emissivity(j)%pdf, numin, numax)

          end do

          e_new(:,j) = e_new(:,j) / integral_loglog(emissivity(j)%x, emissivity(j)%pdf)

       end do

    end if

    allocate(img%dust(dust_id)%emissivity_var(size(emissivity_var)))
    img%dust(dust_id)%emissivity_var = emissivity_var

    allocate(img%dust(dust_id)%log10_emissivity(img%n_nu,size(emissivity_var)))
    do j=1,size(emissivity_var)
       img%dust(dust_id)%log10_emissivity(:,j) = log10(e_new(:,j))
    end do

  end subroutine image_raytracing_set_emissivity

  subroutine image_raytracing_set_opacity(img,dust_id,nu,chi_nu)

    implicit none

    type(image),intent(inout) :: img
    integer,intent(in) :: dust_id
    real(dp),intent(in) :: nu(:), chi_nu(:)
    real(dp), dimension(img%n_nu) :: nu_new, chi_nu_new
    real(dp) :: numin, numax
    integer :: inu

    if(img%use_exact_nu) then

       chi_nu_new = interp1d_loglog(nu, chi_nu, img%nu, bounds_error=.false., fill_value=0._dp)

    else

       do inu=1,img%n_nu

          numin = 10._dp**(img%log10_nu_min + (img%log10_nu_max - img%log10_nu_min) * real(inu-1, dp) / real(img%n_nu, dp))
          numax = 10._dp**(img%log10_nu_min + (img%log10_nu_max - img%log10_nu_min) * real(inu, dp) / real(img%n_nu, dp))

          nu_new(inu) = 10._dp**(0.5_dp * (log10(numin) + log10(numax)))
          chi_nu_new(inu) = integral_loglog(nu, chi_nu, numin, numax) / (numax - numin)

       end do

    end if

    allocate(img%dust(dust_id)%chi(img%n_nu))
    img%dust(dust_id)%chi = chi_nu_new

  end subroutine image_raytracing_set_opacity

end module type_image

max_stars_repo_name: astrofrog/hyperion

max_issues_repo_name: astrofrog/hyperion

max_forks_repo_name: astrofrog/hyperion


text:
! RUN: %python %S/test_errors.py %s %flang_fc1
! Check that non-logical expressions are not allowed.
! Check that non-scalar expressions are not allowed.
! TODO: Insure all non-logicals are prohibited.

LOGICAL, DIMENSION (2) :: B

!ERROR: Must have LOGICAL type, but is REAL(4)
IF (A) A = LOG (A)
!ERROR: Must be a scalar value, but is a rank-1 array
IF (B) A = LOG (A)

END

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
  subroutine setPhotoChannelSigma (sDCM, egamma, photoData)

! ======================================================================
!    Main routine to extract ds/do for channel 1-4 (single pion
!    production):
!
!    Written by K. K. Gudima, Fall 2003?
!    Modified by AJS, May, 2004.
!    Modified by KKG, Sep., 2004
!    Edited by AJS, January, 2005.
!    Modified by AJS, February, 2005.
!    Edited by AJS, LANL T-2, December, 2011.
!    Edited by LMK, XCP-3, July 2013 (included error protection).
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64
    use standardDCMParams, only: zro, one, two, radianToDegree, emnucg
    use standardDCMData,   only: domo2, domo2i, gppipn, gppi0p, elg, &
         & thetai, photonEG

    implicit none
    class(StandardDCM),            intent(inout) :: sDCM
    real(real64),                  intent(in   ) :: egamma
    type(sDCMPhotonCrossSections), intent(  out) :: photoData

    integer(int32) :: ie, ieg1, ieg2, j, jch
    real(real64)   :: deleg, den1, eg1, eg2, sig1, sig2, sint, temp

    real(real64), dimension(2, 20) :: sig = zro
    real(real64), dimension(2, 19) :: r   = zro

! ======================================================================

    ! Obtain cross section for gamma energy in each channel at each angle
    do jch = 1,2

       ! Obtain ieg1/2 values for nearest cross section data
       if (egamma <= elg(jch,2)) then
          ieg1 = 2
          ieg2 = 3
       elseif (egamma >= elg(jch,50)) then
          ieg1 = 49
          ieg2 = 50
       else
          do ie = 3,50
             if (egamma >= elg(jch,ie-1) .and. egamma <= elg(jch,ie)) then
                ieg1 = ie - 1
                ieg2 = ie
                go to 10
             endif
          enddo
       endif

10     continue

       ! Check ieg1/2 values for validity (i.e. within range)
       if (ieg1 < 2 .or. ieg1 > 49 .or. ieg2 < 3 .or. ieg2 > 50) then
          write(sDCM%io%message,2000) ieg1, ieg2
          call sDCM%io%print(3, 3, sDCM%io%message)
          write(sDCM%io%message,2010)
          call sDCM%io%print(3, 3, sDCM%io%message)

          ! Obtain nearest index
          if ( ieg1 <  2 ) ieg1 =  2
          if ( ieg1 > 49 ) ieg1 = 49
          if ( ieg2 <  3 ) ieg2 =  3
          if ( ieg2 > 50 ) ieg2 = 50
       endif

       ! Obtain interpolation energies for the channel
       eg1 = elg(jch,ieg1)
       eg2 = elg(jch,ieg2)
       temp = eg2 - eg1
       if (temp < div0Lim .and. temp > -div0Lim) then
          temp = div0Lim
          write(sDCM%io%message,1000) "116"
          call sDCM%io%print(4, 3, sDCM%io%message)
       end if

       ! Obtain interpolation values (den1 and deleg)
       den1 = zro
       if (eg2.ne.eg1) den1 = one/(temp)
       deleg = egamma - eg1
       sint = zro

       ! Obtain cross section for all bins
       do j = 1,19

          ! Obtain points nearest data for the channel
          if (jch == 1)     then
             sig1 = gppipn(ieg1,j)
             sig2 = gppipn(ieg2,j)
          elseif (jch == 2) then
             sig1 = gppi0p(ieg1,j)
             sig2 = gppi0p(ieg2,j)
          endif

          ! Perform linear interpolation
          sig(jch,j) = sig1 + (sig2 - sig1)*(den1*deleg)

          ! Calculate the channel's angle-integrated cross section
          if (j >= 2) then
             sint = sint + domo2(j)*(sig(jch,j-1) + sig(jch,j))
          endif
          r(jch,j) = sint

       enddo
       sig(jch,20) = sint

       ! Normalize 'r" for all bins in the channel
       do j = 1,19
          if (sint > zro) r(jch,j) = r(jch,j)/sint
       enddo

    end do

    ! Obtain "fine" mesh cross section based on the "course" mesh cross sections
    do jch = 1,2
       sint = zro
       do j = 1,181
          ! Interpolate for angle "j" based on course mesh (theta, sig)
          photoData%si(jch,j) = photonEG%qintxs(thetai(j), sig, jch, 2, 19)

          ! Calculate angle-integrated cross section
          if (j >= 2) then
             sint = sint + domo2i(j)*(photoData%si(jch,j-1) + photoData%si(jch,j))
          endif
          photoData%ri(jch,j) = sint
       end do
       photoData%si(jch,182) = sint

       ! Normalize "ri" for each angle
       do j = 1,181
          if (sint > zro) photoData%ri(jch,j) = photoData%ri(jch,j)/sint
       enddo

    end do

    return

! ======================================================================
1000 format("Divide by zero error prevented in ", &
          & "'setPhotoChannelSigma.f90', line(s)", A)
2000 format("Photon cross section interpolation bounds (ieg1=", i3, &
          & ", ieg2=", i3, ") could not be obtained.")
2010 format("   Will round to nearest valid bound...")
! ======================================================================
  end subroutine setPhotoChannelSigma

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


text:
module riemann_mod
  use amr_parameters
  use hydro_parameters
  use const
  USE io_mod
  USE trace_mod
  USE mpi_mod
  implicit none
  private
  type riemann_t
  contains
    procedure, nopass:: approx
    procedure, nopass:: acoustic
    procedure, nopass:: llf
    procedure, nopass:: hll
    procedure, nopass:: hllc
    procedure, nopass:: hllc_a
    procedure, nopass:: hllc_v
  end type
  type(riemann_t), public:: rieman
contains

!===============================================================================
!===============================================================================
subroutine approx(qleft,qright,fgdnv,ngrid,nvar)
  integer::ngrid,nvar
  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,qgdnv,fgdnv

  ! local arrays
  real(dp),dimension(1:ngrid)::rl   ,ul   ,pl   ,cl
  real(dp),dimension(1:ngrid)::rr   ,ur   ,pr   ,cr
  real(dp),dimension(1:ngrid)::ro   ,uo   ,po   ,co
  real(dp),dimension(1:ngrid)::rstar,ustar,pstar,cstar
  real(dp),dimension(1:ngrid)::wl   ,wr   ,wo
  real(dp),dimension(1:ngrid)::sgnm ,spin ,spout,ushock
  real(dp),dimension(1:ngrid)::frac ,delp ,pold
  integer ,dimension(1:ngrid)::ind  ,ind2

  ! local variables
  real(dp)::smallp, gamma6, ql, qr, usr, usl, wwl, wwr, smallpp, entho, etot
  integer ::i, j, n, iter, n_new

  ! constants
  smallp = smallc**2/gamma
  smallpp = smallr*smallp
  gamma6 = (gamma+one)/(two*gamma)
  entho = one/(gamma-one)

  ! Pressure, density and velocity
  do i=3,ngrid-1
     rl(i)=MAX(qleft (i,1),smallr)
     ul(i)=    qleft (i,2)
     pl(i)=MAX(qleft (i,3),rl(i)*smallp)
     rr(i)=MAX(qright(i,1),smallr)
     ur(i)=    qright(i,2)
     pr(i)=MAX(qright(i,3),rr(i)*smallp)
  end do

  ! Lagrangian sound speed
  do i=3,ngrid-1
     cl(i) = gamma*pl(i)*rl(i)
     cr(i) = gamma*pr(i)*rr(i)
  end do

  ! First guess
  do i=3,ngrid-1
     wl(i)= sqrt(cl(i)); wr(i) = sqrt(cr(i))
     pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i)))/(wl(i)+wr(i))
     pstar(i) = MAX(pstar(i),0.0_dp)
     pold (i)= pstar(i)
  end do
  n = ngrid
  do i = 1,n
     ind(i)=i
  end do

  ! Newton-Raphson iterations to find pstar at the required accuracy
  ! for a two-shock Riemann problem
  do iter = 1,niter_riemann
     do i=1,n
        wwl=sqrt(cl(ind(i))*(one+gamma6*(pold(i)-pl(ind(i)))/pl(ind(i))))
        wwr=sqrt(cr(ind(i))*(one+gamma6*(pold(i)-pr(ind(i)))/pr(ind(i))))
        ql=two*wwl**3/(wwl**2+cl(ind(i)))
        qr=two*wwr**3/(wwr**2+cr(ind(i)))
        usl=ul(ind(i))-(pold(i)-pl(ind(i)))/wwl
        usr=ur(ind(i))+(pold(i)-pr(ind(i)))/wwr
        delp(i)=MAX(qr*ql/(qr+ql)*(usl-usr),-pold(i))
     end do
     do i=1,n
        pold(i)=pold(i)+delp(i)
     end do
     ! Convergence indicator
     do i=1,n
        uo(i)=ABS(delp(i)/(pold(i)+smallpp))
     end do
     n_new=0
     do i=1,n
        if(uo(i)>1.d-06)then
           n_new=n_new+1
           ind2(n_new)=ind (i)
           po  (n_new)=pold(i)
        end if
     end do
     j=n_new
     do i=1,n
        if(uo(i)<=1.d-06)then
           n_new=n_new+1
           ind2(n_new)=ind (i)
           po  (n_new)=pold(i)
        end if
     end do
     ind (1:n)=ind2(1:n)
     pold(1:n)=po  (1:n)
     n=j
  end do

  ! Star region pressure
  ! for a two-shock Riemann problem
  do i=3,ngrid-1
     pstar(ind(i))=pold(i)
  end do
  do i=3,ngrid-1
     wl(i) = sqrt(cl(i)*(one+gamma6*(pstar(i)-pl(i))/pl(i)))
     wr(i) = sqrt(cr(i)*(one+gamma6*(pstar(i)-pr(i))/pr(i)))
  end do

  ! Star region velocity
  ! for a two shock Riemann problem
  do i=3,ngrid-1
     ustar(i) = half*(ul(i) + (pl(i)-pstar(i))/wl(i) + &
          &           ur(i) - (pr(i)-pstar(i))/wr(i) )
  end do

  ! Left going or right going contact wave
  do i=3,ngrid-1
     sgnm(i) = sign(one,ustar(i))
  end do

  ! Left or right unperturbed state
  do i=3,ngrid-1
     if(sgnm(i)==one)then
        ro(i) = rl(i)
        uo(i) = ul(i)
        po(i) = pl(i)
        wo(i) = wl(i)
     else
        ro(i) = rr(i)
        uo(i) = ur(i)
        po(i) = pr(i)
        wo(i) = wr(i)
     end if
  end do
  do i=3,ngrid-1
     co(i) = max(smallc,sqrt(abs(gamma*po(i)/ro(i))))
  end do

  ! Star region density
  do i=3,ngrid-1
     if(pstar(i)>= po(i))then
        ! Shock
        rstar(i) = ro(i)/(one+ro(i)*(po(i)-pstar(i))/wo(i)**2)
     else
        ! Rarefaction
        rstar(i) = ro(i)*(pstar(i)/po(i))**(one/gamma)
     end if
  end do

  do i=3,ngrid-1
     ! Prevent vacuum formation in star region
     rstar(i) = max(rstar(i),smallr)
     ! Star region sound speed
     cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i)))
     cstar(i) = max(cstar(i),smallc)
     ! Compute rarefaction head and tail speed
     spout(i) = co   (i)-sgnm(i)*uo   (i)
     spin (i) = cstar(i)-sgnm(i)*ustar(i)
     ! Compute shock speed
     ushock(i) = wo(i)/ro(i)-sgnm(i)*uo(i)
  end do

  do i=3,ngrid-1
     if(pstar(i)>=po(i))then
        spout(i)=ushock(i)
        spin (i)=spout (i)
     end if
  end do

  ! Sample the solution at x/t=0
  do i=3,ngrid-1
     if(spout(i)<=zero)then
        qgdnv(i,1) = ro(i)
        qgdnv(i,2) = uo(i)
        qgdnv(i,3) = po(i)
     else if(spin(i)>=zero)then
        qgdnv(i,1) = rstar(i)
        qgdnv(i,2) = ustar(i)
        qgdnv(i,3) = pstar(i)
     else
        frac(i)=spout(i)/(spout(i)-spin(i))
        qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i)
        qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i)
        qgdnv(i,1) = ro(i)*(qgdnv(i,3)/po(i))**(one/gamma)
     end if
  end do

  ! Passive scalars
  do n = 4,nvar
     do i=3,ngrid-1
        if(sgnm(i)==one)then
           qgdnv(i,n) = qleft(i,n)
        else
           qgdnv(i,n) = qright(i,n)
        end if
     end do
  end do

  ! Compute fluxes
  do i = 3,ngrid-1
     fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2)  ! Mass density
     fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2  ! Normal momentum
     etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2
     etot = etot + half*qgdnv(i,1)*qgdnv(i,4)**2
     etot = etot + half*qgdnv(i,1)*qgdnv(i,5)**2
     fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3))     ! Total energy
  end do
  ! Other advected quantities
  do n = 4, nvar
     do i = 3,ngrid-1
        fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n)
     end do
  end do
end subroutine approx

!===============================================================================
!===============================================================================
subroutine acoustic(qleft,qright,fgdnv,ngrid,nvar)
  integer::ngrid,nvar
  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,qgdnv,fgdnv

  ! local variables
  integer::i,n
  real(dp)::smallp, entho, etot

  ! local arrays
  real(dp),dimension(1:ngrid)::rl   ,ul   ,pl   ,cl
  real(dp),dimension(1:ngrid)::rr   ,ur   ,pr   ,cr
  real(dp),dimension(1:ngrid)::ro   ,uo   ,po   ,co
  real(dp),dimension(1:ngrid)::rstar,ustar,pstar,cstar
  real(dp),dimension(1:ngrid)::wl   ,wr   ,wo
  real(dp),dimension(1:ngrid)::sgnm ,spin ,spout,ushock
  real(dp),dimension(1:ngrid)::frac

  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)

  ! Initial states pressure, density and velocity
  do i=3,ngrid-1
     rl(i)=max(qleft (i,1),smallr)
     ul(i)=    qleft (i,2)
     pl(i)=max(qleft (i,3),rl(i)*smallp)
     rr(i)=max(qright(i,1),smallr)
     ur(i)=    qright(i,2)
     pr(i)=max(qright(i,3),rr(i)*smallp)
  end do

  ! Acoustic star state
  do i=3,ngrid-1
     cl(i) = sqrt(gamma*pl(i)/rl(i))
     cr(i) = sqrt(gamma*pr(i)/rr(i))
     wl(i) = cl(i)*rl(i)
     wr(i) = cr(i)*rr(i)
     pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i))) &
          &   / (wl(i)+wr(i))
     ustar(i) = ((wr(i)*ur(i)+wl(i)*ul(i))+(pl(i)-pr(i))) &
          &   / (wl(i)+wr(i))
!!$     pstar(i) = MAX(pstar(i),zero)
  end do

  ! Left going or right going contact wave
  do i=3,ngrid-1
     sgnm(i) = sign(one,ustar(i))
  end do

  ! Left or right unperturbed state
  do i=3,ngrid-1
     if(sgnm(i)==one)then
        ro(i) = rl(i)
        uo(i) = ul(i)
        po(i) = pl(i)
        wo(i) = wl(i)
        co(i) = cl(i)
     else
        ro(i) = rr(i)
        uo(i) = ur(i)
        po(i) = pr(i)
        wo(i) = wr(i)
        co(i) = cr(i)
     end if
  end do

  ! Star region density and sound speed
  do i=3,ngrid-1
     rstar(i) = ro(i)+(pstar(i)-po(i))/co(i)**2
     rstar(i) = max(rstar(i),smallr)
     cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i)))
     cstar(i) = max(cstar(i),smallc)
  end do

  ! Head and tail speed of rarefaction
  do i=3,ngrid-1
     spout(i) = co   (i)-sgnm(i)*uo   (i)
     spin (i) = cstar(i)-sgnm(i)*ustar(i)
  end do

  ! Shock speed
  do i=3,ngrid-1
     ushock(i) = half*(spin(i)+spout(i))
     ushock(i) = max(ushock(i),-sgnm(i)*ustar(i))
  end do
  do i=3,ngrid-1
     if(pstar(i)>=po(i))then
        spout(i)=ushock(i)
        spin (i)=spout (i)
     end if
  end do

  ! Sample the solution at x/t=0
  do i=3,ngrid-1
     if(spout(i)<zero)then      ! Initial state
        qgdnv(i,1) = ro(i)
        qgdnv(i,2) = uo(i)
        qgdnv(i,3) = po(i)
     else if(spin(i)>=zero)then  ! Star region
        qgdnv(i,1) = rstar(i)
        qgdnv(i,2) = ustar(i)
        qgdnv(i,3) = pstar(i)
     else                        ! Rarefaction
        frac(i) = spout(i)/(spout(i)-spin(i))
        qgdnv(i,1) = frac(i)*rstar(i) + (one - frac(i))*ro(i)
        qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i)
        qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i)
     end if
  end do

  ! Passive scalars
  do n = 4,nvar
     do i=3,ngrid-1
        if(sgnm(i)==one)then
           qgdnv(i,n) = qleft (i,n)
        else
           qgdnv(i,n) = qright(i,n)
        end if
     end do
  end do

  ! Compute fluxes
  do i = 3,ngrid-1
     fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2)  ! Mass density
     fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2  ! Normal momentum
     etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2
     etot = etot             + half*qgdnv(i,1)*qgdnv(i,4)**2
     etot = etot             + half*qgdnv(i,1)*qgdnv(i,5)**2
     fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3))     ! Total energy
  end do
  ! Other advected quantities
  do n = 4, nvar
     do i = 3,ngrid-1
        fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n)
     end do
  end do
end subroutine acoustic

!===============================================================================
!===============================================================================
subroutine llf(qleft,qright,fgdnv,ngrid,nvar)
  integer::ngrid,nvar
  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv

  ! local arrays
  real(dp),dimension(1:ngrid,1:nvar)::fleft,fright,uleft,uright
  real(dp),dimension(1:ngrid)::cmax

  ! local variables
  integer::i,n
  real(dp)::smallp, entho
  real(dp)::rl   ,ul   ,pl   ,cl
  real(dp)::rr   ,ur   ,pr   ,cr

  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)

  ! Maximum wave speed
  do i=3,ngrid-1
     rl=max(qleft (i,1),smallr)
     ul=    qleft (i,2)
     pl=max(qleft (i,3),rl*smallp)
     rr=max(qright(i,1),smallr)
     ur=    qright(i,2)
     pr=max(qright(i,3),rr*smallp)
     cl= sqrt(gamma*pl/rl)
     cr= sqrt(gamma*pr/rr)
     cmax(i)=max(abs(ul)+cl,abs(ur)+cr)
  end do

  ! Compute conservative variables
  do i = 3,ngrid-1
     ! Mass density
     uleft (i,1) = qleft (i,1)
     uright(i,1) = qright(i,1)
     ! Normal momentum
     uleft (i,2) = qleft (i,1)*qleft (i,2)
     uright(i,2) = qright(i,1)*qright(i,2)
     ! Total energy
     uleft (i,3) = qleft (i,3)*entho + half*qleft (i,1)*qleft (i,2)**2
     uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2
     uleft (i,3) = uleft (i,3)       + half*qleft (i,1)*qleft (i,4)**2
     uright(i,3) = uright(i,3)       + half*qright(i,1)*qright(i,4)**2
     uleft (i,3) = uleft (i,3)       + half*qleft (i,1)*qleft (i,5)**2
     uright(i,3) = uright(i,3)       + half*qright(i,1)*qright(i,5)**2
  end do
  ! Other advected quantities
  do n = 4, nvar
     do i = 3,ngrid-1
        uleft (i,n) = qleft (i,1)*qleft (i,n)
        uright(i,n) = qright(i,1)*qright(i,n)
     end do
  end do

  ! Compute left and right fluxes
  do i = 3,ngrid-1
     ! Mass density
     fleft (i,1) = uleft (i,2)
     fright(i,1) = uright(i,2)
     ! Normal momentum
     fleft (i,2) = qleft (i,3)+uleft (i,2)*qleft (i,2)
     fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2)
     ! Total energy
     fleft (i,3) = qleft (i,2)*(uleft (i,3)+qleft (i,3))
     fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3))
  end do
  ! Other advected quantities
  do n = 4, nvar
     do i = 3,ngrid-1
        fleft (i,n) = fleft (i,1)*qleft (i,n)
        fright(i,n) = fright(i,1)*qright(i,n)
     end do
  end do

  ! Compute Lax-Friedrich fluxes
  do n = 1, nvar
     do i = 3,ngrid-1
        fgdnv(i,n) = half*(fleft(i,n)+fright(i,n)-cmax(i)*(uright(i,n)-uleft(i,n)))
     end do
  end do
end subroutine llf

!===============================================================================
!===============================================================================
subroutine hllc (qleft,qright,fgdnv,ngrid,nvar,isothermal,error,detailed_timer)
  integer::ngrid,nvar
  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
  REAL(dp)::SL,SR
  REAL(dp)::entho
  REAL(dp)::rl,pl,ul,ecinl,etotl,ptotl
  REAL(dp)::rr,pr,ur,ecinr,etotr,ptotr
  REAL(dp)::cfastl,rcl,rstarl
  REAL(dp)::cfastr,rcr,rstarr
  REAL(dp)::etotstarl,etotstarr
  REAL(dp)::ustar,ptotstar
  REAL(dp)::ro,uo,ptoto,etoto
  REAL(dp)::smallp
  INTEGER ::ivar, i
  LOGICAL :: isothermal,error,detailed_timer
  INTEGER, SAVE:: itimer=0
  !
  if (detailed_timer) call trace_begin ('riemann::hllc', itimer=itimer)
  ! constants
  smallp = smallc**2/gamma
  if (isothermal .or. gamma==1d0) then
    entho=one
  else
    entho=one/(gamma-one)
  end if
  error = .false.
  !
  if (io%verbose>1) then
    do i=3,ngrid-1
      !
      ! Left variables
      rl=max(qleft (i,1),smallr)
      Pl=max(qleft (i,3),rl*smallp)
      ul=    qleft (i,2)
      !
      ecinl = half*rl*(ul*ul+qleft(i,4)**2+qleft(i,5)**2)
      etotl = Pl*entho+ecinl
      Ptotl = Pl
      !
      ! Right variables
      rr=max(qright(i,1),smallr)
      Pr=max(qright(i,3),rr*smallp)
      ur=    qright(i,2)
      !
      ecinr = half*rr*(ur*ur+qright(i,4)**2+qright(i,5)**2)
      etotr = Pr*entho+ecinr
      Ptotr = Pr
      !
      ! Find the largest eigenvalues in the normal direction to the interface
      cfastl=sqrt(max(gamma*Pl/rl,smallc**2))
      cfastr=sqrt(max(gamma*Pr/rr,smallc**2))
      !
      ! Compute HLL wave speed
      SL=min(ul,ur)-max(cfastl,cfastr)
      SR=max(ul,ur)+max(cfastl,cfastr)
      !
      ! Compute lagrangian sound speed
      rcl=rl*(ul-SL)
      rcr=rr*(SR-ur)
      if ((rcl+rcl)==0.0) then
        print 1,'HLLC: rl,Pl,ul,ecinl,etotl,Ptotl,cfastl,SL =', &
          i,rl,Pl,ul,ecinl,etotl,Ptotl,cfastl,SL,SL-ul
        print 1,'HLLC: rr,Pr,ur,ecinr,etotr,Ptotr,cfastr,SR =', &
          i,rr,Pr,ur,ecinr,etotr,Ptotr,cfastr,SR,SR-ur
        1 format(a,i4,1p,10e12.3)
        error = .true.
        return
      end if
      !
      ! Compute acoustic star state
      ustar   =(rcr*ur   +rcl*ul   +  (Ptotl-Ptotr))/(rcr+rcl)
      Ptotstar=(rcr*Ptotl+rcl*Ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
      if ((SL-ustar)==0.0 .or. (SR-ustar)==0.0) then
        print 1,'HLLC: rl,Pl,ul,ecinl,etotl,Ptotl,cfastl,SL =', &
          i,rl,Pl,ul,ecinl,etotl,Ptotl,cfastl,SL,SL-ul
        print 1,'HLLC: rr,Pr,ur,ecinr,etotr,Ptotr,cfastr,SR =', &
          i,rr,Pr,ur,ecinr,etotr,Ptotr,cfastr,SR,SR-ur
        print 1, 'HLLC: ustar,Ptotstar', i,ustar,Ptotstar
        call mpi%abort ('HLLC')
        error = .true.
        return
      end if
      !
      ! Left star region variables
      rstarl=rl*(SL-ul)/(SL-ustar)
      etotstarl=((SL-ul)*etotl-Ptotl*ul+Ptotstar*ustar)/(SL-ustar)
      !
      ! Right star region variables
      rstarr=rr*(SR-ur)/(SR-ustar)
      etotstarr=((SR-ur)*etotr-Ptotr*ur+Ptotstar*ustar)/(SR-ustar)
      !
      ! Sample the solution at x/t=0
      if(SL>0d0)then
         ro=rl
         uo=ul
         Ptoto=Ptotl
         etoto=etotl
      else if(ustar>0d0)then
         ro=rstarl
         uo=ustar
         Ptoto=Ptotstar
         etoto=etotstarl
      else if (SR>0d0)then
         ro=rstarr
         uo=ustar
         Ptoto=Ptotstar
         etoto=etotstarr
      else
         ro=rr
         uo=ur
         Ptoto=Ptotr
         etoto=etotr
      end if
      !
      ! Compute the Godunov flux
      fgdnv(i,1) = ro*uo
      fgdnv(i,2) = ro*uo*uo+Ptoto
      fgdnv(i,3) = (etoto+Ptoto)*uo
    end do
  else
    do i=3,ngrid-1
      !
      ! Left variables
      rl=max(qleft (i,1),smallr)
      Pl=max(qleft (i,3),rl*smallp)
      ul=    qleft (i,2)
      !
      ecinl = half*rl*(ul*ul+qleft(i,4)**2+qleft(i,5)**2)
      etotl = Pl*entho+ecinl
      Ptotl = Pl
      !
      ! Right variables
      rr=max(qright(i,1),smallr)
      Pr=max(qright(i,3),rr*smallp)
      ur=    qright(i,2)
      !
      ecinr = half*rr*(ur*ur+qright(i,4)**2+qright(i,5)**2)
      etotr = Pr*entho+ecinr
      Ptotr = Pr
      !
      ! Find the largest eigenvalues in the normal direction to the interface
      cfastl=sqrt(max(gamma*Pl/rl,smallc**2))
      cfastr=sqrt(max(gamma*Pr/rr,smallc**2))
      !
      ! Compute HLL wave speed
      SL=min(ul,ur)-max(cfastl,cfastr)
      SR=max(ul,ur)+max(cfastl,cfastr)
      !
      ! Compute lagrangian sound speed
      rcl=rl*(ul-SL)
      rcr=rr*(SR-ur)
      !
      ! Compute acoustic star state
      ustar   =(rcr*ur   +rcl*ul   +  (Ptotl-Ptotr))/(rcr+rcl)
      Ptotstar=(rcr*Ptotl+rcl*Ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
      !
      ! Left star region variables
      rstarl=rl*(SL-ul)/(SL-ustar)
      etotstarl=((SL-ul)*etotl-Ptotl*ul+Ptotstar*ustar)/(SL-ustar)
      !
      ! Right star region variables
      rstarr=rr*(SR-ur)/(SR-ustar)
      etotstarr=((SR-ur)*etotr-Ptotr*ur+Ptotstar*ustar)/(SR-ustar)
      !
      ! Sample the solution at x/t=0
!define MERGE
#ifdef MERGE
      ro=rr
      uo=ur
      Ptoto=Ptotr
      etoto=etotr
      ro   =merge(rstarr   ,ro     ,SR>0d0)
      uo   =merge(ustar    ,uo     ,SR>0d0)
      Ptoto=merge(Ptotstar ,Ptoto  ,SR>0d0)
      etoto=merge(etotstarr,etoto  ,SR>0d0)
      ro   =merge(rstarr   ,ro     ,ustar>0d0)
      uo   =merge(ustar    ,uo     ,ustar>0d0)
      Ptoto=merge(Ptotstar ,Ptoto  ,ustar>0d0)
      etoto=merge(etotstarl,etoto  ,ustar>0d0)
      ro   =merge(rl       ,ro     ,SL>0d0)
      uo   =merge(ul       ,uo     ,SL>0d0)
      Ptoto=merge(Ptotl    ,Ptoto  ,SL>0d0)
      etoto=merge(etotl    ,etoto  ,SL>0d0)
#else
      if(SL>0d0)then
         ro=rl
         uo=ul
         Ptoto=Ptotl
         etoto=etotl
      else if(ustar>0d0)then
         ro=rstarl
         uo=ustar
         Ptoto=Ptotstar
         etoto=etotstarl
      else if (SR>0d0)then
         ro=rstarr
         uo=ustar
         Ptoto=Ptotstar
         etoto=etotstarr
      else
         ro=rr
         uo=ur
         Ptoto=Ptotr
         etoto=etotr
      end if
#endif
      !
      ! Compute the Godunov flux
      fgdnv(i,1) = ro*uo
      fgdnv(i,2) = ro*uo*uo+Ptoto
      fgdnv(i,3) = (etoto+Ptoto)*uo
    end do
  end if
  !
  do ivar=4,nvar
    do i=3,ngrid-1
      fgdnv(i,ivar) = fgdnv(i,1)*merge(qleft(i,ivar),qright(i,ivar),fgdnv(i,1)>0)
    end do
  end do
  if (detailed_timer) call trace_end (itimer)
end subroutine hllc

!===============================================================================
!===============================================================================
subroutine hllc_a (qleft,qright,fgdnv,ngrid,nvar)
  integer::ngrid,nvar
  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
  REAL(dp)::entho
  REAL(dp)::pl,ecinl
  REAL(dp)::pr,ecinr
  REAL(dp)::cfastl,rcl
  REAL(dp)::cfastr,rcr
  REAL(dp)::smallp
  REAL(dp),dimension(ngrid)::SL,SR
  REAL(dp),dimension(ngrid)::ro,uo,Ptoto,Etoto
  REAL(dp),dimension(ngrid)::rl,ul,Ptotl,Etotl
  REAL(dp),dimension(ngrid)::rr,ur,Ptotr,Etotr
  REAL(dp),dimension(ngrid)::rstarl,Etotstarl
  REAL(dp),dimension(ngrid)::rstarr,Etotstarr
  REAL(dp),dimension(ngrid)::ustar,Ptotstar
  INTEGER ::ivar, i
  !
  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)
  !
  do i=3,ngrid-1
     !
     ! Left variables
     rl(i)=max(qleft (i,1),smallr)
     Pl=max(qleft (i,3),rl(i)*smallp)
     ul(i)=    qleft (i,2)
     !
     ecinl = half*rl(i)*(ul(i)*ul(i)+qleft(i,4)**2+qleft(i,5)**2)
     etotl(i) = Pl*entho+ecinl
     Ptotl(i) = Pl
     !
     ! Right variables
     rr(i)=max(qright(i,1),smallr)
     Pr=max(qright(i,3),rr(i)*smallp)
     ur(i)=    qright(i,2)
     !
     ecinr = half*rr(i)*(ur(i)*ur(i)+qright(i,4)**2+qright(i,5)**2)
     etotr(i) = Pr*entho+ecinr
     Ptotr(i) = Pr
     !
     ! Find the largest eigenvalues in the normal direction to the interface
     cfastl=sqrt(max(gamma*Pl/rl(i),smallc**2))
     cfastr=sqrt(max(gamma*Pr/rr(i),smallc**2))
     !
     ! Compute HLL wave speed
     SL(i)=min(ul(i),ur(i))-max(cfastl,cfastr)
     SR(i)=max(ul(i),ur(i))+max(cfastl,cfastr)
     !
     ! Compute lagrangian sound speed
     rcl=rl(i)*(ul(i)-SL(i))
     rcr=rr(i)*(SR(i)-ur(i))
     !
     ! Compute acoustic star state
     ustar(i)   =(rcr*ur(i)   +rcl*ul(i)   +  (Ptotl(i)-Ptotr(i)))/(rcr+rcl)
     Ptotstar(i)=(rcr*Ptotl(i)+rcl*Ptotr(i)+rcl*rcr*(ul(i)-ur(i)))/(rcr+rcl)
  end do
  do i=3,ngrid-1
     !
     ! Left star region variables
     rstarl(i)=rl(i)*(SL(i)-ul(i))/(SL(i)-ustar(i))
     etotstarl(i)=((SL(i)-ul(i))*etotl(i)-Ptotl(i)*ul(i)+Ptotstar(i)*ustar(i))/(SL(i)-ustar(i))
     !
     ! Right star region variables
     rstarr(i)=rr(i)*(SR(i)-ur(i))/(SR(i)-ustar(i))
     etotstarr(i)=((SR(i)-ur(i))*etotr(i)-Ptotr(i)*ur(i)+Ptotstar(i)*ustar(i))/(SR(i)-ustar(i))
  end do
     !
     ! Sample the solution at x/t=0
  do i=3,ngrid-1
     if(SL(i)>0d0)then
        ro(i)=rl(i)
        uo(i)=ul(i)
        Ptoto(i)=Ptotl(i)
        etoto(i)=etotl(i)
     else if(ustar(i)>0d0)then
        ro(i)=rstarl(i)
        uo(i)=ustar(i)
        Ptoto(i)=Ptotstar(i)
        etoto(i)=etotstarl(i)
     else if (SR(i)>0d0)then
        ro(i)=rstarr(i)
        uo(i)=ustar(i)
        Ptoto(i)=Ptotstar(i)
        etoto(i)=etotstarr(i)
     else
        ro(i)=rr(i)
        uo(i)=ur(i)
        Ptoto(i)=Ptotr(i)
        etoto(i)=etotr(i)
     end if
  end do
     !
     ! Compute the Godunov flux
  do i=3,ngrid-1
     fgdnv(i,1) = ro(i)*uo(i)
     fgdnv(i,2) = ro(i)*uo(i)*uo(i)+Ptoto(i)
     fgdnv(i,3) = (etoto(i)+Ptoto(i))*uo(i)
  end do
     !
  do ivar=4,nvar
  do i=3,ngrid-1
     fgdnv(i,ivar) = merge(fgdnv(i,1)*qleft (i,ivar), &
                           fgdnv(i,1)*qright(i,ivar), fgdnv(i,1)>0)
  end do
  end do
end subroutine hllc_a

!===============================================================================
!===============================================================================
subroutine hllc_v (qleft,qright,fgdnv,ngrid,nvar)
  integer::ngrid,nvar
  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
  REAL(dp),dimension(ngrid)::SL,SR
  REAL(dp),dimension(ngrid)::entho
  REAL(dp),dimension(ngrid)::rl,pl,ul,ecinl,etotl,ptotl
  REAL(dp),dimension(ngrid)::rr,pr,ur,ecinr,etotr,ptotr
  REAL(dp),dimension(ngrid)::cfastl,rcl,rstarl
  REAL(dp),dimension(ngrid)::cfastr,rcr,rstarr
  REAL(dp),dimension(ngrid)::etotstarl,etotstarr
  REAL(dp),dimension(ngrid)::ustar,ptotstar
  REAL(dp),dimension(ngrid)::ro,uo,ptoto,etoto
  REAL(dp),dimension(ngrid)::smallp
  INTEGER ::ivar, i
  !
  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)
  !
     !
     ! Left variables
     rl=max(qleft (:,1),smallr)
     Pl=max(qleft (:,3),rl*smallp)
     ul=    qleft (:,2)
     !
     ecinl = half*rl*(ul*ul+qleft(:,4)**2+qleft(:,5)**2)
     etotl = Pl*entho+ecinl
     Ptotl = Pl
     !
     ! Right variables
     rr=max(qright(:,1),smallr)
     Pr=max(qright(:,3),rr*smallp)
     ur=    qright(:,2)
     !
     ecinr = half*rr*(ur*ur+qright(:,4)**2+qright(:,5)**2)
     etotr = Pr*entho+ecinr
     Ptotr = Pr
     !
     ! Find the largest eigenvalues in the normal direction to the interface
     cfastl=sqrt(max(gamma*Pl/rl,smallc**2))
     cfastr=sqrt(max(gamma*Pr/rr,smallc**2))
     !
     ! Compute HLL wave speed
     SL=min(ul,ur)-max(cfastl,cfastr)
     SR=max(ul,ur)+max(cfastl,cfastr)
     !
     ! Compute lagrangian sound speed
     rcl=rl*(ul-SL)
     rcr=rr*(SR-ur)
     !
     ! Compute acoustic star state
     ustar   =(rcr*ur   +rcl*ul   +  (Ptotl-Ptotr))/(rcr+rcl)
     Ptotstar=(rcr*Ptotl+rcl*Ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
     !
     ! Left star region variables
     rstarl=rl*(SL-ul)/(SL-ustar)
     etotstarl=((SL-ul)*etotl-Ptotl*ul+Ptotstar*ustar)/(SL-ustar)
     !
     ! Right star region variables
     rstarr=rr*(SR-ur)/(SR-ustar)
     etotstarr=((SR-ur)*etotr-Ptotr*ur+Ptotstar*ustar)/(SR-ustar)
     !
     ! Sample the solution at x/t=0
     where (SL>0d0)
        ro=rl
        uo=ul
        Ptoto=Ptotl
        etoto=etotl
     else where (ustar>0d0)
        ro=rstarl
        uo=ustar
        Ptoto=Ptotstar
        etoto=etotstarl
     else where  (SR>0d0)
        ro=rstarr
        uo=ustar
        Ptoto=Ptotstar
        etoto=etotstarr
     else where
        ro=rr
        uo=ur
        Ptoto=Ptotr
        etoto=etotr
     end where
     !
     ! Compute the Godunov flux
     fgdnv(:,1) = ro*uo
     fgdnv(:,2) = ro*uo*uo+Ptoto
     fgdnv(:,3) = (etoto+Ptoto)*uo
     do ivar = 4,nvar
        where (fgdnv(:,1)>0)
           fgdnv(:,ivar) = fgdnv(:,1)*qleft (:,ivar)
        else where
           fgdnv(:,ivar) = fgdnv(:,1)*qright(:,ivar)
        end where
     end do
end subroutine hllc_v

!===============================================================================
!===============================================================================
subroutine hll(qleft,qright,fgdnv,ngrid,nvar)
  integer::ngrid,nvar
  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
  real(dp),dimension(1:ngrid,1:nvar)::fleft,fright,uleft,uright
  real(dp),dimension(1:ngrid)::SL,SR
  integer::i,n
  real(dp)::smallp, entho
  real(dp)::rl   ,ul   ,pl   ,cl
  real(dp)::rr   ,ur   ,pr   ,cr

  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)

  ! Maximum wave speed
  do i=3,ngrid-1
     rl=max(qleft (i,1),smallr)
     ul=    qleft (i,2)
     pl=max(qleft (i,3),rl*smallp)
     rr=max(qright(i,1),smallr)
     ur=    qright(i,2)
     pr=max(qright(i,3),rr*smallp)
     cl= sqrt(gamma*pl/rl)
     cr= sqrt(gamma*pr/rr)
     SL(i)=min(min(ul,ur)-max(cl,cr),zero)
     SR(i)=max(max(ul,ur)+max(cl,cr),zero)
  end do

  ! Compute conservative variables
  do i = 3,ngrid-1
     uleft (i,1) = qleft (i,1)
     uright(i,1) = qright(i,1)
     uleft (i,2) = qleft (i,1)*qleft (i,2)
     uright(i,2) = qright(i,1)*qright(i,2)
     uleft (i,3) = qleft (i,3)*entho + half*qleft (i,1)*qleft (i,2)**2
     uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2
     uleft (i,3) = uleft (i,3)       + half*qleft (i,1)*qleft (i,4)**2
     uright(i,3) = uright(i,3)       + half*qright(i,1)*qright(i,4)**2
     uleft (i,3) = uleft (i,3)       + half*qleft (i,1)*qleft (i,5)**2
     uright(i,3) = uright(i,3)       + half*qright(i,1)*qright(i,5)**2
  end do
  ! Other advected quantities
  do n = 4, nvar
     do i = 3,ngrid-1
        uleft (i,n) = qleft (i,1)*qleft (i,n)
        uright(i,n) = qright(i,1)*qright(i,n)
     end do
  end do

  ! Compute left and right fluxes
  do i = 3,ngrid-1
     fleft (i,1) = uleft (i,2)
     fright(i,1) = uright(i,2)
     fleft (i,2) = qleft (i,3)+uleft (i,2)*qleft (i,2)
     fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2)
     fleft (i,3) = qleft (i,2)*(uleft (i,3)+qleft (i,3))
     fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3))
  end do
  ! Other advected quantities
  do n = 4, nvar
     do i = 3,ngrid-1
        fleft (i,n) = fleft (i,1)*qleft (i,n)
        fright(i,n) = fright(i,1)*qright(i,n)
     end do
  end do

  ! Compute HLL fluxes
  do n = 1, nvar
     do i = 3,ngrid-1
        fgdnv(i,n) = (SR(i)*fleft(i,n)-SL(i)*fright(i,n) &
             & + SR(i)*SL(i)*(uright(i,n)-uleft(i,n)))/(SR(i)-SL(i))
     end do
  end do
end subroutine hll

end module riemann_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
      SUBROUTINE free_mem_ns_par(lreset)
      USE vmec_main
      USE realspace
      USE vforces
      USE xstuff
      USE csplinx
      USE fbal

      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      LOGICAL, INTENT(in) :: lreset

C-----------------------------------------------
      IF (ALLOCATED(pshalf))  DEALLOCATE(pshalf)
      IF (ALLOCATED(psqrts))  DEALLOCATE(psqrts)
      IF (ALLOCATED(pwint))   DEALLOCATE(pwint)
C     <MJL 2021-02-12>
      IF (ALLOCATED(pwint_ns))   DEALLOCATE(pwint_ns)
      IF (ALLOCATED(ireflect_par)) DEALLOCATE(ireflect_par)
C     </MJL>
      IF (ALLOCATED(pfaclam)) DEALLOCATE(pfaclam)
      IF (ALLOCATED(pchip))   DEALLOCATE(pchip)
      IF (ALLOCATED(pphip))   DEALLOCATE(pphip)
      IF (ALLOCATED(pfaclam)) DEALLOCATE(pfaclam)

      IF (ALLOCATED(pxc)     .and. lreset) DEALLOCATE(pxc)
      IF (ALLOCATED(pscalxc) .and. lreset) DEALLOCATE (pscalxc)

      IF (ALLOCATED(pxstore))    DEALLOCATE (pxstore)
      IF (ALLOCATED(pxcdot))     DEALLOCATE (pxcdot)
      IF (ALLOCATED(pxsave))     DEALLOCATE (pxsave)
      IF (ALLOCATED(pgc))        DEALLOCATE (pgc)
      IF (ALLOCATED(pcol_scale)) DEALLOCATE(pcol_scale)

      END SUBROUTINE free_mem_ns_par

      SUBROUTINE free_mem_ns(lreset)
      USE vmec_main
      USE realspace
      USE vforces
      USE xstuff
      USE csplinx
      USE fbal

      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      LOGICAL, INTENT(in) :: lreset
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: istat1 = 0, istat2 = 0, istat3 = 0, istat4 = 0,
     1           istat5 = 0, istat6 = 0, istat7 = 0, istat8 = 0,
     2           istat9 = 0, istat10 = 0
C-----------------------------------------------

      IF (ALLOCATED(phip)) THEN
         DEALLOCATE(phip, chip, shalf, sqrts, wint, stat=istat1)
      END IF

      IF (ALLOCATED(ireflect)) THEN
         DEALLOCATE(ireflect, stat=istat2)
      END IF

      IF (ALLOCATED(ard)) THEN
         DEALLOCATE(ard, arm, brd, brm, crd, azd, azm, bzd, bzm, sm,
     &              sp, bmin, bmax, stat=istat5)
      END IF

!      IF (ALLOCATED(rru_fac)) THEN
!        DEALLOCATE(rru_fac,frcc_fac,fzsc_fac,icurv,vpphi,presgrad,
!     &             r01,z01,bdamp, stat=istat6)
      IF (ALLOCATED(iotaf)) THEN
         DEALLOCATE(iotaf, mass, phi, presf, jcuru, jcurv, jdotb, buco)
         DEALLOCATE(bvco, bucof, bvcof, chi)
!#ifdef _ANIMEC
!     &              phot, pmap, pppr, papr, tpotb, pd,
!#endif
         DEALLOCATE(bdotgradv, equif, specw, tcon, psi, yellip, yinden)
         DEALLOCATE(ytrian, yshift, ygeo, overr, faclam, iotas, phips)
         DEALLOCATE(chips, pres, vp, beta_vol, jperp2, jpar2, bdotb)
         DEALLOCATE(clam, blam, dlam, phipf, chipf, rru_fac, rzu_fac)
         DEALLOCATE(frcc_fac, fzsc_fac, icurv, vpphi, presgrad, r01)
         DEALLOCATE(z01, bdamp, stat=istat6)
      END IF

      IF (ALLOCATED(rmidx)) THEN
         DEALLOCATE(rmidx, hmidx, wmidx, qmidx, tenmidx, ymidx, y2midx,
     &              stat=istat7)
      END IF

      IF (ALLOCATED(gc)) THEN
         DEALLOCATE(gc, xsave, xstore, xcdot, col_scale, stat=istat8)
      END IF
      IF (ALLOCATED(xc) .AND. lreset) THEN
         DEALLOCATE(xc, scalxc)
      END IF

      IF ((istat1 + istat2 + istat3 + istat4 + istat5 + istat6 +
     &     istat7 + istat8) .ne. 0) THEN
          PRINT *,' deallocation problem in free_mem_ns'
          PRINT *,' istat1 = ',istat1,' istat2 = ',istat2
          PRINT *,' istat3 = ',istat3,' istat4 = ',istat4
          PRINT *,' istat5 = ',istat5,' istat6 = ',istat6
          PRINT *,' istat7 = ',istat7,' istat8 = ',istat8
      END IF

      END SUBROUTINE free_mem_ns


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE DBMIRR (NELBS, NELBE, IDELB, NUMELB, NUMLNK, LINK,
     *  BLKTYP, NDIM, NONQUD)
C=======================================================================

C   --*** DBMIRR *** (GEN3D) Fixup connectivity for reflections
C   --   Written by Greg Sjaardema - revised 02/10/89
C   --   Modified from DBOELB Written by Amy Gilkey
C   --
C   --Parameters:
C   --   NELBS, NELBE - IN - the number of first and last element blocks
C   --   IDELB - IN - the element block IDs for each block
C   --   NUMELB - IN - the number of elements in each block
C   --   NUMLNK - IN - the number of nodes per element in each block
C   --   LINK - IN/OUT - the connectivity for each block
C   --

      include 'exodusII.inc'
      INTEGER IDELB(*)
      INTEGER NUMELB(*)
      INTEGER NUMLNK(*)
      INTEGER LINK(*)
      CHARACTER*(MXSTLN) BLKTYP(*)
      LOGICAL NONQUD

      NONQUD = .FALSE.
      IELNK = 0

      DO 100 NELB = NELBS, NELBE
         IELB = NELB-NELBS+1

         ISLNK = IELNK + 1
         IELNK = IELNK + NUMLNK(IELB) * NUMELB(IELB)

         CALL DBMIR1 (IDELB(IELB), NUMELB(IELB),
     *     NUMLNK(IELB), LINK(ISLNK), BLKTYP(IELB), NDIM, NONQUD)
  100 CONTINUE

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas