!
!     ---------------------------------------------------------------
!      I T R E X G 2
!     ---------------------------------------------------------------
!
!     Written by G. Gaigalas,                                      *
!     Vanderbilt University,  Nashville             October 1996   *
!
      INTEGER FUNCTION ITREXG2 (I1, I2, I3, I4, K)
!...Translated by Pacific-Sierra Research 77to90  4.3E  10:51:45  11/16/01
!...Switches:
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: I1
      INTEGER , INTENT(IN) :: I2
      INTEGER , INTENT(IN) :: I3
      INTEGER , INTENT(IN) :: I4
      INTEGER , INTENT(OUT) :: K
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: J
!-----------------------------------------------
      J = MIN0(IABS(I1 - I2),IABS(I3 - I4))
      K = MIN0(IABS(I1 + I2),IABS(I3 + I4)) - J + 1
      ITREXG2 = J
      RETURN
      END FUNCTION ITREXG2


text:
Creative Media is a web and graphic design firm based at the University of California, Davis.  Its mission is to help students and businesses with their publicity and communication needs through lowcost and memorable websites, eyecatching print materials, and creative marketing strategies. The Graphics office is located on the third floor of the Memorial Union next to the north elevator. The rest of Creative Media is down the hall.

 Projects

Creative Media manages websites for:
  ASUCD (http://asucd.ucdavis.edu site)  and every subsite of ASUCD.  See also ASUCD Website.
  AggieTV (http://aggietv.org)
  The California Aggie (http://www.californiaaggie.com site)
  Community Housing Listing (http://chl.ucdavis.edu/ site)
  Experimental College (http://asucd.ucdavis.edu/experimentalcollege/ site)
  Unitrans (http://unitrans.ucdavis.edu/ site)
  And many more.

Creative Media also handles:
  Everything for AGTV
  Calculating results for ASUCD Elections
  Designing logos, fliers, posters, Tshirts, buttons, signs, banners, and print publications for ASUCD units and a select few outside organizations.

 Staff 20092010
  Director: Alex Park  the man on top; a career employee
  Graphics Director: Stephanie Pan
  Web Programmers: Byron Dover, Bryan Stitt, Andrew Theiss, and Chris Thielan
  Graphic Artists:  Leslie Cheng, Kevin Chu, and Caitlin Flint

 History

Creative Media was the combination of two different former units of ASUCD:  Network Services (name?) and AS Graphics.  Thus it combined the normal IT department with the group that designed ASUCDs publications and publicity, and put under a staff director (Alex Park).  The unit used to be called Media Link, but was renamed to Creative Media a few years ago.  (Despite the name, its still not considered a media unit in the ASUCD Bylaws.)

Due to the lack of checks or balances on the unit, Jonathon Leathers called for an audit of Creative Media during his winter 2005 senate campaign. Ballsy. He did not win a seat on senate, and the audit did not occur. AGTV spent some time under the California Aggie journalistic microscope in winter 2005 after the ASUCD ASUCD Senate Senate passed a Shadiness Factor bill that would drop $33,000 on Apple G5 computers and plasmascreen monitors for AGTV.  CM never actually spent the money, and ended up buying cheaper screens (~$11,000) from another Senate Bill.

Creative Media and the Student Government Administrative Office were the subject of an ASUCD Court http://www.ucdgreens.org/complaint_aggie_article.shtml complaint filed in March of 2004. The complaint alleged Creative Media and SGAO failed to update the ASUCD Website with documents.  However, its not actually Creative Medias responsibility, its SGAOs responsibility to keep student government uptodate.

Controller Paul Harms found an unused $8,000 line item in the Creative Media budget that past controllers had forgotten to remove.  This sort of small problem is typical in ASUCDs finances and  should not be considered a reflection on the Unit Director.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module xsec

  use sdata, only: dp
  implicit none
  save

contains

  !******************************************************************************!

  SUBROUTINE XS_updt (xbcon, xftem, xmtem, xcden, xbpos)
  !
  ! Purpose:
  !    To update XS for given TH paramaters and rod position
  !

  USE sdata, ONLY : nnod, mat, sigtr, siga, nuf, sigf, sigs, dc, &
  get_time, xs_time
  use io, only: bbcon, bftem, bmtem, bcden, bcrod, bcbcs, bcrod, bxtab

  IMPLICIT NONE

  REAL(DP), INTENT(IN) :: xbcon  ! Provided Boron Concentration
  REAL(DP), DIMENSION(:), INTENT(IN) :: xftem  ! Provided fuel temperature
  REAL(DP), DIMENSION(:), INTENT(IN) :: xmtem  ! Provided moderator temperature
  REAL(DP), DIMENSION(:), INTENT(IN) :: xcden  ! Provided coolant density
  REAL(DP), DIMENSION(:), INTENT(IN) :: xbpos  ! Provided control rod bank position

  integer :: n

  REAL(DP) :: st, fn

  st = get_time()

  if (bxtab == 0) then
    CALL base_updt()
    IF (bbcon == 1 .OR. bcbcs == 1) CALL bcon_updt(xbcon)
    IF (bftem == 1) CALL ftem_updt(xftem)
    IF (bmtem == 1) CALL mtem_updt(xmtem)
    IF (bcden == 1) CALL cden_updt(xcden)
    IF (bcrod == 1) CALL crod_updt(xbpos)
    CALL Dsigr_updt()
  else
    DO n = 1, nnod
      CALL brInterp(0, mat(n), xcden(n),  xbcon, xftem(n), xmtem(n), sigtr(n,:), &
      siga(n,:), nuf(n,:), sigf(n,:), sigs(n,:,:), dc(n,:,:))
    END DO
    IF (bcrod == 1) CALL crod_tab_updt(xbpos)
    CALL Dsigr_updt()
  end if

  call check_xs()

  fn = get_time()
  xs_time = xs_time + (fn-st)

  END SUBROUTINE XS_updt

  !******************************************************************************!

  SUBROUTINE check_xs ()
  !
  ! Purpose:
  !    To check errors in xsec
  !

  use sdata, only : D, sigr, nuf, chi, sigs, nnod, ng, mat, nmat
  use io,    only : ounit

  IMPLICIT NONE

  INTEGER :: n, g, h

  DO g = 1, ng
    DO n = 1, nnod
      IF (D(n,g) < 1.e-20) then
        WRITE(*,*)
        WRITE(*,3541) mat(n), g
        WRITE(*,*) ' DIFFUSION COEF. IS CLOSE TO ZERO OR NEGATIVE'
        WRITE(ounit,3541) mat(n), g
        WRITE(ounit,*) ' DIFFUSION COEF. IS CLOSE TO ZERO OR NEGATIVE'
        STOP
      END IF
      IF (sigr(n,g) < 0.) then
        WRITE(*,*)
        WRITE(*,3542) mat(n), g
        WRITE(*,*) ' REMOVAL XS IS  NEGATIVE'
        WRITE(ounit,3542) mat(n), g
        WRITE(ounit,*) ' REMOVAL XS IS NEGATIVE'
        STOP
      END IF
      IF (nuf(n,g) < 0.) then
        WRITE(*,*)
        WRITE(*,3544) mat(n), g
        WRITE(*,*) ' NU*FISSION XS IS  NEGATIVE'
        WRITE(ounit,3544) mat(n), g
        WRITE(ounit,*) ' NU*FISSION XS IS NEGATIVE'
        STOP
      END IF
    END DO
    DO n = 1, nmat
      IF (chi(n,g) < 0.) then
        WRITE(*,*)
        WRITE(*,3543) n, g
        WRITE(*,*) ' FISSION SPECTRUM IS  NEGATIVE'
        WRITE(ounit,3543) n, g
        WRITE(ounit,*) ' FISSION SPECTRUM IS NEGATIVE'
        STOP
      END IF
    END DO
    DO h = 1, ng
      DO n = 1, nnod
        IF (sigs(n,g,h) < 0.) then
          WRITE(*,*)
          WRITE(*,3545) mat(n), g, h
          WRITE(*,*) ' SCATTERING XS IS  NEGATIVE'
          WRITE(ounit,3545) mat(n), g, h
          WRITE(ounit,*) ' SCATTERING XS IS NEGATIVE'
          STOP
        END IF
      END DO
    END DO
  END DO

  3541 FORMAT (2X, 'ERROR IN THE DIFFUSION COEFFICIENT ', &
  'FOR MATERIAL NUMBER : ', I3, ' ENERGY GROUP : ', I2)
  3542 FORMAT (2X, 'ERROR IN THE REMOVAL CROSS SECTION ', &
  'FOR MATERIAL NUMBER : ', I3, ' ENERGY GROUP : ', I2)
  3543 FORMAT (2X, 'ERROR IN THE FISSION SPECTRUM ', &
  'FOR MATERIAL NUMBER : ', I3, ' ENERGY GROUP : ', I2)
  3544 FORMAT (2X, 'ERROR IN THE NU*FISSION XSEC ', &
  'FOR MATERIAL NUMBER : ', I3, ' ENERGY GROUP : ', I2)
  3545 FORMAT (2X, 'ERROR IN THE SCATTERING CROSS SECTION ', &
  'FOR MATERIAL NUMBER : ', I3, ' ENERGY GROUP : ', I2, &
   ' TO ENERGY GROUP : ', I2)

  END SUBROUTINE check_xs

  !******************************************************************************!

  SUBROUTINE base_updt ()
  !
  ! Purpose:
  !    To update current XS to base XS
  !

  USE sdata, ONLY: nnod, sigtr, siga, nuf, sigf, sigs, mat, &
                   xsigtr, xsiga, xnuf, xsigf, xsigs

  IMPLICIT NONE

  INTEGER :: n

  DO n = 1, nnod
      sigtr(n,1:)   = xsigtr(mat(n),1:)
      siga (n,1:)   = xsiga (mat(n),1:)
      nuf  (n,1:)   = xnuf  (mat(n),1:)
      sigf (n,1:)   = xsigf (mat(n),1:)
      sigs (n,1:,1:) = xsigs (mat(n),1:,1:)
  END DO

  END SUBROUTINE base_updt

  !******************************************************************************!

  SUBROUTINE Dsigr_updt ()
  !
  ! Purpose:
  !    To update diffusion constant and removal XS
  !

  USE sdata, ONLY: nnod, ng, sigtr, siga,  &
                   sigs, D, sigr

  IMPLICIT NONE

  INTEGER :: i, g, h
  REAL(DP) :: dum

  DO i = 1, nnod
    DO g = 1, ng
      IF (sigtr(i,g) < 1.e-5 ) STOP "Negative diffusion coefficient encountered"
      D(i,g) = 1._dp/(3._dp*sigtr(i,g))
      dum = 0.
      DO h= 1, ng
          IF (g /= h) dum = dum + sigs(i,g,h)
      END DO
      sigr(i,g) =  siga(i,g) + dum
    END DO
  END DO

  END SUBROUTINE Dsigr_updt

  !******************************************************************************!

  SUBROUTINE crod_updt (bpos)
  !
  ! Purpose: TO UPDATE AND CALCUALTE VOLUME WEIGHTED HOMOGENIZED CX FOR RODDED NODE
  !

  USE sdata, ONLY: ng, nxx, nyy, nzz, xyz, zdel, mat, &
                   sigtr, siga, nuf, sigf, sigs, &
                   dsigtr, dsiga, dnuf, dsigf, dsigs, &
                   coreh, fbmap, pos0, ssize

  IMPLICIT NONE

  REAL(DP), DIMENSION(:), INTENT(IN) :: bpos

  INTEGER ::i, j, k, g, h
  REAL(DP) :: rodh, vfrac
  REAL(DP) :: dum

  ! For each node
  DO j = 1, nyy
    DO i = 1, nxx
       IF (fbmap(i,j) > 0) THEN
          !!!(rodh -> posistion the tip of the control rod the top of core)
           rodh = coreh - pos0  - bpos(fbmap(i,j))*ssize
           dum = 0._dp
           DO k = nzz, 1, -1
             ! For partially rodded node, get volume weighted homogenized CX (0 < vfrac < 1._DP)
             IF (rodh >= dum .AND. rodh <= dum+zdel(k)) THEN   ! If this node partially rodded
               vfrac = (rodh - dum) / zdel(k)
               sigtr(xyz(i,j,k),:) = sigtr(xyz(i,j,k),:) + &
                                  vfrac * dsigtr(mat(xyz(i,j,k)),:)
               siga(xyz(i,j,k),:)  = siga(xyz(i,j,k),:) + &
                                  vfrac * dsiga(mat(xyz(i,j,k)),:)
               nuf(xyz(i,j,k),:)   = nuf(xyz(i,j,k),:) + &
                                  vfrac * dnuf(mat(xyz(i,j,k)),:)
               sigf(xyz(i,j,k),:)  = sigf(xyz(i,j,k),:) + &
                                  vfrac * dsigf(mat(xyz(i,j,k)),:)
               sigs(xyz(i,j,k),:,:)  = sigs(xyz(i,j,k),:,:) + &
                                    vfrac * dsigs(mat(xyz(i,j,k)),:,:)
               EXIT
             END IF
             ! For fully rodded node, vfrac = 1.
             sigtr(xyz(i,j,k),:) = sigtr(xyz(i,j,k),:) + dsigtr(mat(xyz(i,j,k)),:)
             siga(xyz(i,j,k),:)  = siga(xyz(i,j,k),:) + dsiga(mat(xyz(i,j,k)),:)
             nuf(xyz(i,j,k),:)   = nuf(xyz(i,j,k),:) + dnuf(mat(xyz(i,j,k)),:)
             sigf(xyz(i,j,k),:)  = sigf(xyz(i,j,k),:) + dsigf(mat(xyz(i,j,k)),:)
             sigs(xyz(i,j,k),:,:)  = sigs(xyz(i,j,k),:,:) + dsigs(mat(xyz(i,j,k)),:,:)

             dum = dum + zdel(k)
           END DO
           ! if negative CX found, Surpress CX to zero and calculate D and sigr
           DO k = nzz, 1, -1
             DO g = 1, ng
                IF (siga(xyz(i,j,k),g) < 0.) siga(xyz(i,j,k),g) = 0.
                IF (nuf(xyz(i,j,k),g) < 0.) nuf(xyz(i,j,k),g) = 0.
                IF (sigf(xyz(i,j,k),g) < 0.) sigf(xyz(i,j,k),g) = 0.
                DO h = 1, ng
                  IF (sigs(xyz(i,j,k),g,h) < 0.) sigs(xyz(i,j,k),g,h) = 0.
                END DO
             END DO
           END DO
       END IF
    END DO
  END DO

  END SUBROUTINE crod_updt

  !******************************************************************************!

  SUBROUTINE crod_tab_updt (bpos)
  !
  ! Purpose: TO UPDATE AND CALCUALTE VOLUME WEIGHTED HOMOGENIZED CX FOR RODDED NODE
  !

  USE sdata, ONLY: ng, nxx, nyy, nzz, xyz, zdel, mat, &
                   sigtr, siga, nuf, sigf, sigs, dc, &
                   coreh, fbmap, pos0, ssize, m, &
                   nnod, cden, ftem, mtem, bcon

  IMPLICIT NONE

  REAL(DP), DIMENSION(:), INTENT(IN) :: bpos

  INTEGER ::i, j, k, g, h
  REAL(DP) :: rodh, vfrac
  REAL(DP) :: dum

  INTEGER :: nd

  REAL(DP), DIMENSION(nnod,ng) :: rsigtr, rsiga, rnuf, rsigf
  REAL(DP), DIMENSION(nnod,ng,ng) :: rsigs
  REAL(DP), DIMENSION(nnod,ng,6) :: rdc

  DO j = 1, nyy
    DO i = 1, nxx
      IF (fbmap(i,j) > 0) THEN
         !(rodh -> posistion the tip of the control rod the top of core)
         rodh = coreh - pos0  - bpos(fbmap(i,j))*ssize
         dum = 0._DP
           DO k = nzz, 1, -1
             nd = xyz(i,j,k)                                  ! Node number
             IF (m(mat(nd))%trod == 1) THEN
               CALL brInterp(1, mat(nd), cden(nd), bcon, ftem(nd), &
               mtem(nd), rsigtr(nd,:), rsiga(nd,:), rnuf(nd,:), &
               rsigf(nd,:), rsigs(nd,:,:), rdc(nd,:,:))
             ELSE
               WRITE(101,1671) fbmap(i,j), mat(nd)
               WRITE(*,1671) fbmap(i,j), mat(nd)
               STOP
               1671 FORMAT (2X, 'CONTROL ROD BANK NUMBER ', I4, &
               ' COINCIDES WITH MATERIAL NUMBER ', I4, &
               ' THAT DOES NOT HAVE CONTROL ROD DATA IN XTAB FILE')
             END IF
             ! For partially rodded node, get volume weighted homogenized CX (0 < vfrac < 1._DP)
             IF (rodh >= dum .AND. rodh <= dum+zdel(k)) THEN   ! If this node partially rodded
               vfrac = (rodh - dum) / zdel(k)                    ! Rodded fraction
               sigtr(nd,:)  = (1._DP - vfrac) * sigtr(nd,:)  + &
               vfrac * rsigtr(nd,:)
               siga(nd,:)   = (1._DP - vfrac) * siga(nd,:)   + &
               vfrac * rsiga(nd,:)
               nuf(nd,:)    = (1._DP - vfrac) * nuf(nd,:)    + &
               vfrac * rnuf(nd,:)
               sigf(nd,:)   = (1._DP - vfrac) * sigf(nd,:)   + &
               vfrac * rsigf(nd,:)
               sigs(nd,:,:) = (1._DP - vfrac) * sigs(nd,:,:) + &
               vfrac * rsigs(nd,:,:)
               dc(nd,:,:)   = (1._DP - vfrac) * dc(nd,:,:)   + &
               vfrac * rdc(nd,:,:)
               EXIT
             END IF
             ! For fully rodded node, vfrac = 1.
             sigtr(nd,:)  = rsigtr(nd,:)
             siga(nd,:)   = rsiga(nd,:)
             nuf(nd,:)    = rnuf(nd,:)
             sigf(nd,:)   = rsigf(nd,:)
             sigs(nd,:,:) = rsigs(nd,:,:)
             dc(nd,:,:)   = rdc(nd,:,:)

             dum = dum + zdel(k)
           END DO
           ! if negative CX found, Surpress CX to zero  and calculate D and sigr
           DO k = nzz, 1, -1
             nd = xyz(i,j,k)
             DO g = 1, ng
               IF (siga(nd,g) < 0.) siga(nd,g)  = 0.
               IF (nuf(nd,g)  < 0.) nuf(nd,g)   = 0.
               IF (sigf(nd,g) < 0.) sigf(nd,g)  = 0.
               DO h = 1, ng
                 IF (sigs(nd,g,h) < 0.) sigs( nd,g,h) = 0.
               END DO
               DO h = 1, 6
                 IF (dc(nd,g,h) < 0.) dc(nd,g,h) = 0.
               END DO
             END DO
           END DO
       END IF
    END DO
  END DO

  END SUBROUTINE crod_tab_updt

  !******************************************************************************!

  SUBROUTINE bcon_updt (bcon)

  !
  ! Purpose:
  !    To update CX for given boron concentration

  USE sdata, ONLY: nnod, ng, sigtr, siga, nuf, sigf, sigs, mat, &
                   csigtr, csiga, cnuf, csigf, csigs, rbcon

  IMPLICIT NONE

  REAL(DP), INTENT(IN) :: bcon
  INTEGER :: n, g, h

  DO n = 1, nnod
      DO g = 1, ng
          sigtr(n,g) = sigtr(n,g) + csigtr(mat(n),g) * (bcon - rbcon)
          siga(n,g)  = siga(n,g)  + csiga(mat(n),g)  * (bcon - rbcon)
          nuf(n,g)   = nuf(n,g)   + cnuf(mat(n),g)   * (bcon - rbcon)
          sigf(n,g)  = sigf(n,g)  + csigf(mat(n),g)  * (bcon - rbcon)
          DO h = 1, ng
              sigs(n,g,h) = sigs(n,g,h) + csigs(mat(n),g,h) * (bcon - rbcon)
          END DO
      END DO
  END DO

  END SUBROUTINE bcon_updt

  !******************************************************************************!

  SUBROUTINE ftem_updt (ftem)

  !
  ! Purpose:
  !    To update CX for given fuel temp

  USE sdata, ONLY: nnod, ng, sigtr, siga, nuf, sigf, sigs, mat, &
                   fsigtr, fsiga, fnuf, fsigf, fsigs, rftem

  IMPLICIT NONE

  REAL(DP), DIMENSION(:), INTENT(IN) :: ftem
  INTEGER :: n, g, h

  DO n = 1, nnod
      DO g = 1, ng
          sigtr(n,g) = sigtr(n,g) + fsigtr(mat(n),g) * (SQRT(ftem(n))- SQRT(rftem))
          siga(n,g)  = siga(n,g)  + fsiga(mat(n),g)  * (SQRT(ftem(n)) - SQRT(rftem))
          nuf(n,g)   = nuf(n,g)   + fnuf(mat(n),g)   * (SQRT(ftem(n)) - SQRT(rftem))
          sigf(n,g)  = sigf(n,g)  + fsigf(mat(n),g)  * (SQRT(ftem(n)) - SQRT(rftem))
          DO h = 1, ng
             sigs(n,g,h) = sigs(n,g,h) + fsigs(mat(n),g,h) * (SQRT(ftem(n)) - SQRT(rftem))
          END DO
        END DO
  END DO

  END SUBROUTINE ftem_updt

  !******************************************************************************!

  SUBROUTINE mtem_updt (mtem)

  !
  ! Purpose:
  !    To update CX for given moderator temperature

  USE sdata, ONLY: nnod, ng, sigtr, siga, nuf, sigf, sigs, mat, &
                   msigtr, msiga, mnuf, msigf, msigs, rmtem

  IMPLICIT NONE

  REAL(DP), DIMENSION(:), INTENT(IN) :: mtem
  INTEGER :: n, g, h

  DO n = 1, nnod
      DO g = 1, ng
          sigtr(n,g) = sigtr(n,g) + msigtr(mat(n),g) * (mtem(n) - rmtem)
          siga(n,g)  = siga(n,g)  + msiga(mat(n),g)  * (mtem(n) - rmtem)
          nuf(n,g)   = nuf(n,g)   + mnuf(mat(n),g)   * (mtem(n) - rmtem)
          sigf(n,g)  = sigf(n,g)  + msigf(mat(n),g)  * (mtem(n) - rmtem)
          DO h = 1, ng
              sigs(n,g,h) = sigs(n,g,h) + msigs(mat(n),g,h) * (mtem(n) - rmtem)
          END DO
      END DO
  END DO

  END SUBROUTINE mtem_updt

  !******************************************************************************!

  SUBROUTINE cden_updt (cden)

  !
  ! Purpose:
  !    To update CX for given coolant density

  USE sdata, ONLY: nnod, ng, sigtr, siga, nuf, sigf, sigs, mat, &
                   lsigtr, lsiga, lnuf, lsigf, lsigs, rcden

  IMPLICIT NONE

  REAL(DP), DIMENSION(:), INTENT(IN) :: cden
  INTEGER :: n, g, h

  DO n = 1, nnod
      DO g = 1, ng
          sigtr(n,g) = sigtr(n,g) + lsigtr(mat(n),g) * (cden(n) - rcden)
          siga(n,g)  = siga(n,g)  + lsiga(mat(n),g)  * (cden(n) - rcden)
          nuf(n,g)   = nuf(n,g)   + lnuf(mat(n),g)   * (cden(n) - rcden)
          sigf(n,g)  = sigf(n,g)  + lsigf(mat(n),g)  * (cden(n) - rcden)
          DO h = 1, ng
              sigs(n,g,h) = sigs(n,g,h) + lsigs(mat(n),g,h) * (cden(n) - rcden)
          END DO
      END DO
  END DO

  END SUBROUTINE cden_updt

  !******************************************************************************!

  SUBROUTINE brInterp (rod, mn, xcden, xbcon, xftem, xmtem, sigtr, siga, nuf, &
    sigf, sigs, dc)
  !Purpose: To interpolate the xsec data from xtab file for given bcon,
  ! ftem, mtem and cden

  USE sdata, ONLY: m, XBRANCH, ng
  use io,    only: ounit

  IMPLICIT NONE

  INTEGER, INTENT(IN) :: rod, mn  ! CR indicator and material number
  REAL(DP), INTENT(IN) :: xbcon, xftem, xmtem, xcden  ! TH Parameters
  REAL(DP), DIMENSION(:), INTENT(OUT) :: sigtr, siga, nuf, sigf
  REAL(DP), DIMENSION(:,:), INTENT(OUT) :: dc, sigs

  INTEGER, PARAMETER :: nx = 8
  TYPE(XBRANCH), DIMENSION(nx) :: xs   !Temporary xsec for interpolation
  INTEGER :: s, t, u, v, mx
  INTEGER :: s1, s2, t1, t2, u1, u2, v1, v2  ! Dimnesion Position of the given parameters
  INTEGER :: i
  REAL(DP) :: radx

  !Define + and - operators for XBRANCH type addition and substitution respectively
  INTERFACE OPERATOR (+)
    MODULE PROCEDURE brAdd
  END INTERFACE
  INTERFACE OPERATOR (-)
    MODULE PROCEDURE brSubst
  END INTERFACE
  INTERFACE OPERATOR (*)
    MODULE PROCEDURE brRealMult
  END INTERFACE

  ! Set to default
  s1=1; s2=1; t1=1; t2=1; u1=1; u2=1; v1=1; v2=1

  ! Get 2 closest points for interpolation
  ! FOR COOLANT DENSITY
  IF (m(mn)%nd > 1) THEN
    mx = m(mn)%nd
    IF (xcden >= m(mn)%pd(1) .AND. xcden <= m(mn)%pd(mx)) THEN
      DO s = 2, mx
        IF (xcden >= m(mn)%pd(s-1) .AND. xcden <= m(mn)%pd(s)) THEN
          s1 = s-1
          s2 = s
          EXIT
        END IF
      END DO
    ELSE IF (xcden < m(mn)%pd(1) .AND. (m(mn)%pd(1) - xcden) / m(mn)%pd(1) < 0.2) THEN
      s1 = 1
      s2 = 2
    ELSE IF (xcden > m(mn)%pd(mx) .AND. (xcden - m(mn)%pd(mx)) / m(mn)%pd(mx) < 0.2) THEN
      s1 = mx - 1
      s2 = mx
    ELSE
      WRITE(ounit,1567) xcden
      WRITE(*,1567) xcden
      STOP
      1567 FORMAT(2X, '  ERROR: COOLANT DENSITY ', F7.3 ,' g/cm3 &
      & IS OUT OF THE RANGE OF THE BRANCH PARAMETER')
    END IF
  END IF

  ! FOR BORON CONCENTRATION
  IF (m(mn)%nb > 1) THEN
    mx = m(mn)%nb
    IF (xbcon >= m(mn)%pb(1) .AND. xbcon <= m(mn)%pb(mx)) THEN
      DO t = 2, mx
        IF (xbcon >= m(mn)%pb(t-1) .AND. xbcon <= m(mn)%pb(t)) THEN
          t1 = t-1
          t2 = t
          EXIT
        END IF
      END DO
    ELSE IF (xbcon < m(mn)%pb(1) .AND. (m(mn)%pb(1) - xbcon) < 100.) THEN
      t1 = 1
      t2 = 2
    ELSE IF (xbcon > m(mn)%pb(mx) .AND. (xbcon - m(mn)%pb(mx)) < 100.) THEN
      t1 = mx - 1
      t2 = mx
    ELSE
      WRITE(ounit,1568) xbcon
      WRITE(*,1568) xbcon
      STOP
      1568 FORMAT(2X, '  ERROR: BORON CONCENTRATION ', F8.1 ,' PPM &
      & IS OUT OF THE RANGE OF THE BRANCH PARAMETER')
    END IF
  END IF

  ! FOR FUEL TEMPERATURE
  IF (m(mn)%nf > 1) THEN
    mx = m(mn)%nf
    IF (xftem >= m(mn)%pf(1) .AND. xftem <= m(mn)%pf(mx)) THEN
      DO u = 2, mx
        IF (xftem >= m(mn)%pf(u-1) .AND. xftem <= m(mn)%pf(u)) THEN
          u1 = u-1
          u2 = u
          EXIT
        END IF
      END DO
    ELSE IF (xftem < m(mn)%pf(1) .AND. (m(mn)%pf(1) - xftem) / m(mn)%pf(1) < 0.2) THEN
      u1 = 1
      u2 = 2
    ELSE IF (xftem > m(mn)%pf(mx) .AND. (xftem - m(mn)%pf(mx)) / m(mn)%pf(mx) < 0.2) THEN
      u1 = mx - 1
      u2 = mx
    ELSE
      WRITE(ounit,1570) xftem
      WRITE(*,1570) xftem
      STOP
      1570 FORMAT(2X, '  ERROR: FUEL TEMPERATURE ', F7.1 ,' K &
      & IS OUT OF THE RANGE OF THE BRANCH PARAMETER')
    END IF
  END IF

  ! FOR MODERATOR TEMPERATURE
  IF (m(mn)%nm > 1) THEN
    mx = m(mn)%nm
    IF (xmtem >= m(mn)%pm(1) .AND. xmtem <= m(mn)%pm(mx)) THEN
      DO v = 2, mx
        IF (xmtem >= m(mn)%pm(v-1) .AND. xmtem <= m(mn)%pm(v)) THEN
          v1 = v-1
          v2 = v
          EXIT
        END IF
      END DO
    ELSE IF (xmtem < m(mn)%pm(1) .AND. (m(mn)%pm(1) - xmtem) / m(mn)%pm(1) < 0.2) THEN
      v1 = 1
      v2 = 2
    ELSE IF (xmtem > m(mn)%pm(mx) .AND. (xmtem - m(mn)%pm(mx)) / m(mn)%pm(mx) < 0.2) THEN
      v1 = mx - 1
      v2 = mx
    ELSE
      WRITE(ounit,1569) xmtem
      WRITE(*,1569) xmtem
      STOP
      1569 FORMAT(2X, '  ERROR: MODERATOR TEMPERATURE ', F7.1 ,' K &
      & IS OUT OF THE RANGE OF THE BRANCH PARAMETER')
    END IF
  END IF

  !Start doing interpolation
  DO i = 1, nx   !Allocate memory to xs
    ALLOCATE(xs(i)%sigtr(ng))
    ALLOCATE(xs(i)%siga(ng))
    ALLOCATE(xs(i)%nuf(ng))
    ALLOCATE(xs(i)%sigf(ng))
    ALLOCATE(xs(i)%dc(ng,6))
    ALLOCATE(xs(i)%sigs(ng,ng))
  END DO

  IF (rod == 0) THEN   ! For Unrodded XSEC
    !interpolation on Moderator Temperature
    IF (m(mn)%nm > 1) THEN
      radx = (xmtem - m(mn)%pm(v1)) / (m(mn)%pm(v2) - m(mn)%pm(v1))
      xs(1) = m(mn)%xsec(s1,t1,u1,v1) + &
      radx * (m(mn)%xsec(s1,t1,u1,v2) - m(mn)%xsec(s1,t1,u1,v1))
      xs(2) = m(mn)%xsec(s1,t1,u2,v1) + &
      radx * (m(mn)%xsec(s1,t1,u2,v2) - m(mn)%xsec(s1,t1,u2,v1))
      xs(3) = m(mn)%xsec(s1,t2,u1,v1) + &
      radx * (m(mn)%xsec(s1,t2,u1,v2) - m(mn)%xsec(s1,t2,u1,v1))
      xs(4) = m(mn)%xsec(s1,t2,u2,v1) + &
      radx * (m(mn)%xsec(s1,t2,u2,v2) - m(mn)%xsec(s1,t2,u2,v1))
      xs(5) = m(mn)%xsec(s2,t1,u1,v1) + &
      radx * (m(mn)%xsec(s2,t1,u1,v2) - m(mn)%xsec(s2,t1,u1,v1))
      xs(6) = m(mn)%xsec(s2,t1,u2,v1) + &
      radx * (m(mn)%xsec(s2,t1,u2,v2) - m(mn)%xsec(s2,t1,u2,v1))
      xs(7) = m(mn)%xsec(s2,t2,u1,v1) + &
      radx * (m(mn)%xsec(s2,t2,u1,v2) - m(mn)%xsec(s2,t2,u1,v1))
      xs(8) = m(mn)%xsec(s2,t2,u2,v1) + &
      radx * (m(mn)%xsec(s2,t2,u2,v2) - m(mn)%xsec(s2,t2,u2,v1))
    ELSE
      xs(1) = m(mn)%xsec(s1,t1,u1,v1)
      xs(2) = m(mn)%xsec(s1,t1,u2,v1)
      xs(3) = m(mn)%xsec(s1,t2,u1,v1)
      xs(4) = m(mn)%xsec(s1,t2,u2,v1)
      xs(5) = m(mn)%xsec(s2,t1,u1,v1)
      xs(6) = m(mn)%xsec(s2,t1,u2,v1)
      xs(7) = m(mn)%xsec(s2,t2,u1,v1)
      xs(8) = m(mn)%xsec(s2,t2,u2,v1)
    END IF
    !interpolation on Fuel Temperature
    IF (m(mn)%nf > 1) THEN
      radx = (xftem - m(mn)%pf(u1)) / (m(mn)%pf(u2) - m(mn)%pf(u1))
      xs(1) = xs(1) + radx * (xs(2) - xs(1))
      xs(3) = xs(3) + radx * (xs(4) - xs(3))
      xs(5) = xs(5) + radx * (xs(6) - xs(5))
      xs(7) = xs(7) + radx * (xs(8) - xs(7))
    END IF
    !interpolation on Boron concentration
    IF (m(mn)%nb > 1) THEN
      radx = (xbcon - m(mn)%pb(t1)) / (m(mn)%pb(t2) - m(mn)%pb(t1))
      xs(1) = xs(1) + radx * (xs(3) - xs(1))
      xs(5) = xs(5) + radx * (xs(7) - xs(5))
    END IF
    !interpolation on coolant density
    IF (m(mn)%nd > 1) THEN
      xs(1) = xs(1) + (xcden - m(mn)%pd(s1)) / (m(mn)%pd(s2) - m(mn)%pd(s1)) * &
      (xs(5) - xs(1))
    END IF
  ELSE   ! For Rodded XSEC
    !interpolation on Moderator Temperature
    IF (m(mn)%nm > 1) THEN
      radx = (xmtem - m(mn)%pm(v1)) / (m(mn)%pm(v2) - m(mn)%pm(v1))
      xs(1) = m(mn)%rxsec(s1,t1,u1,v1) + &
      radx * (m(mn)%rxsec(s1,t1,u1,v2) - m(mn)%rxsec(s1,t1,u1,v1))
      xs(2) = m(mn)%rxsec(s1,t1,u2,v1) + &
      radx * (m(mn)%rxsec(s1,t1,u2,v2) - m(mn)%rxsec(s1,t1,u2,v1))
      xs(3) = m(mn)%rxsec(s1,t2,u1,v1) + &
      radx * (m(mn)%rxsec(s1,t2,u1,v2) - m(mn)%rxsec(s1,t2,u1,v1))
      xs(4) = m(mn)%rxsec(s1,t2,u2,v1) + &
      radx * (m(mn)%rxsec(s1,t2,u2,v2) - m(mn)%rxsec(s1,t2,u2,v1))
      xs(5) = m(mn)%rxsec(s2,t1,u1,v1) + &
      radx * (m(mn)%rxsec(s2,t1,u1,v2) - m(mn)%rxsec(s2,t1,u1,v1))
      xs(6) = m(mn)%rxsec(s2,t1,u2,v1) + &
      radx * (m(mn)%rxsec(s2,t1,u2,v2) - m(mn)%rxsec(s2,t1,u2,v1))
      xs(7) = m(mn)%rxsec(s2,t2,u1,v1) + &
      radx * (m(mn)%rxsec(s2,t2,u1,v2) - m(mn)%rxsec(s2,t2,u1,v1))
      xs(8) = m(mn)%rxsec(s2,t2,u2,v1) + &
      radx * (m(mn)%rxsec(s2,t2,u2,v2) - m(mn)%rxsec(s2,t2,u2,v1))
    ELSE
      xs(1) = m(mn)%rxsec(s1,t1,u1,v1)
      xs(2) = m(mn)%rxsec(s1,t1,u2,v1)
      xs(3) = m(mn)%rxsec(s1,t2,u1,v1)
      xs(4) = m(mn)%rxsec(s1,t2,u2,v1)
      xs(5) = m(mn)%rxsec(s2,t1,u1,v1)
      xs(6) = m(mn)%rxsec(s2,t1,u2,v1)
      xs(7) = m(mn)%rxsec(s2,t2,u1,v1)
      xs(8) = m(mn)%rxsec(s2,t2,u2,v1)
    END IF
    !interpolation on Fuel Temperature
    IF (m(mn)%nf > 1) THEN
      radx = (xftem - m(mn)%pf(u1)) / (m(mn)%pf(u2) - m(mn)%pf(u1))
      xs(1) = xs(1) + radx * (xs(2) - xs(1))
      xs(3) = xs(3) + radx * (xs(4) - xs(3))
      xs(5) = xs(5) + radx * (xs(6) - xs(5))
      xs(7) = xs(7) + radx * (xs(8) - xs(7))
    END IF
    !interpolation on Boron concentration
    IF (m(mn)%nb > 1) THEN
      radx = (xbcon - m(mn)%pb(t1)) / (m(mn)%pb(t2) - m(mn)%pb(t1))
      xs(1) = xs(1) + radx * (xs(3) - xs(1))
      xs(5) = xs(5) + radx * (xs(7) - xs(5))
    END IF

    !interpolation on coolant density
    IF (m(mn)%nd > 1) THEN
      xs(1) = xs(1) + (xcden - m(mn)%pd(s1)) / (m(mn)%pd(s2) - m(mn)%pd(s1)) * &
      (xs(5) - xs(1))
    END IF
  END IF
  sigtr = xs(1)%sigtr
  siga = xs(1)%siga
  nuf = xs(1)%nuf
  sigf = xs(1)%sigf
  sigs = xs(1)%sigs
  dc   = xs(1)%dc

  DO i = 1, nx   !DeAllocate memory to xs
    DEALLOCATE(xs(i)%sigtr)
    DEALLOCATE(xs(i)%siga)
    DEALLOCATE(xs(i)%nuf)
    DEALLOCATE(xs(i)%sigf)
    DEALLOCATE(xs(i)%dc)
    DEALLOCATE(xs(i)%sigs)
  END DO

  END SUBROUTINE brInterp

  !******************************************************************************!

  FUNCTION brAdd(A, B) RESULT (C)

    ! To perform XBRANCH data type addition

  USE sdata, ONLY: XBRANCH

  IMPLICIT NONE

  TYPE(XBRANCH), INTENT(IN) :: A, B
  TYPE(XBRANCH) :: C

  C%sigtr = A%sigtr + B%sigtr
  C%siga = A%siga + B%siga
  C%nuf = A%nuf + B%nuf
  C%sigf = A%sigf + B%sigf
  C%sigs = A%sigs + B%sigs
  C%dc = A%dc + B%dc

  END FUNCTION brAdd

  !******************************************************************************!

  FUNCTION brSubst(A, B) RESULT (C)

      ! To perform XBRANCH data type substraction

  USE sdata, ONLY: XBRANCH

  IMPLICIT NONE

  TYPE(XBRANCH), INTENT(IN) :: A, B
  TYPE(XBRANCH) :: C

  C%sigtr = A%sigtr - B%sigtr
  C%siga = A%siga - B%siga
  C%nuf = A%nuf - B%nuf
  C%sigf = A%sigf - B%sigf
  C%sigs = A%sigs - B%sigs
  C%dc = A%dc - B%dc

  END FUNCTION brSubst

  !******************************************************************************!

  FUNCTION brRealMult(Re, A) RESULT (B)

      ! To perform XBRANCH data type substraction

  USE sdata, ONLY: XBRANCH, DP

  IMPLICIT NONE

  REAL(DP), INTENT(IN) :: Re
  TYPE(XBRANCH), INTENT(IN) :: A
  TYPE(XBRANCH) :: B

  B%sigtr = Re * A%sigtr
  B%siga = Re * A%siga
  B%nuf = Re * A%nuf
  B%sigf = Re * A%sigf
  B%sigs = Re * A%sigs
  B%dc = Re * A%dc

  END FUNCTION brRealMult

  !******************************************************************************!

end module

max_stars_repo_name: imronuke/KOMODO

max_issues_repo_name: imronuke/KOMODO

max_forks_repo_name: imronuke/KOMODO


text:
! These routines do the characteristic tracing under the parabolic
! profiles in each zone to the edge / half-time.

module trace_ppm_module

  implicit none

  private

  public trace_ppm

contains

  subroutine trace_ppm(q,c,flatn,qd_l1,qd_l2,qd_h1,qd_h2, &
                       dloga,dloga_l1,dloga_l2,dloga_h1,dloga_h2, &
                       qxm,qxp,qym,qyp,qpd_l1,qpd_l2,qpd_h1,qpd_h2, &
                       srcQ,src_l1,src_l2,src_h1,src_h2, &
                       gamc,gc_l1,gc_l2,gc_h1,gc_h2, &
                       ilo1,ilo2,ihi1,ihi2,dx,dy,dt)

    use network, only : nspec
    use eos_type_module
    use eos_module
    use bl_constants_module
    use meth_params_module, only : QVAR, QRHO, QU, QV, QREINT, QPRES, &
         QTEMP, QFS, QFX, QGAME, &
         small_dens, small_pres, &
         ppm_type, ppm_trace_sources, ppm_temp_fix, &
         ppm_reference_eigenvectors, &
         ppm_predict_gammae, &
         npassive, qpass_map
    use ppm_module, only : ppm

    implicit none

    integer ilo1,ilo2,ihi1,ihi2
    integer qd_l1,qd_l2,qd_h1,qd_h2
    integer dloga_l1,dloga_l2,dloga_h1,dloga_h2
    integer qpd_l1,qpd_l2,qpd_h1,qpd_h2
    integer src_l1,src_l2,src_h1,src_h2
    integer gc_l1,gc_l2,gc_h1,gc_h2

    double precision     q(qd_l1:qd_h1,qd_l2:qd_h2,QVAR)
    double precision     c(qd_l1:qd_h1,qd_l2:qd_h2)
    double precision flatn(qd_l1:qd_h1,qd_l2:qd_h2)
    double precision dloga(dloga_l1:dloga_h1,dloga_l2:dloga_h2)

    double precision qxm(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR)
    double precision qxp(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR)
    double precision qym(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR)
    double precision qyp(qpd_l1:qpd_h1,qpd_l2:qpd_h2,QVAR)

    double precision  srcQ(src_l1:src_h1,src_l2:src_h2,QVAR)
    double precision gamc(gc_l1:gc_h1,gc_l2:gc_h2)

    double precision dx, dy, dt

    ! Local variables
    integer :: i, j, iwave, idim
    integer :: n, ipassive

    double precision :: hdt, dtdx, dtdy

    ! To allow for easy integration of radiation, we adopt the
    ! following conventions:
    !
    ! rho : mass density
    ! u, v, w : velocities
    ! p : gas (hydro) pressure
    ! ptot : total pressure (note for pure hydro, this is
    !        just the gas pressure)
    ! rhoe_g : gas specific internal energy
    ! cgas : sound speed for just the gas contribution
    ! cc : total sound speed (including radiation)
    ! h_g : gas specific enthalpy / cc**2
    ! gam_g : the gas Gamma_1
    ! game : gas gamma_e
    !
    ! for pure hydro, we will only consider:
    !   rho, u, v, w, ptot, rhoe_g, cc, h_g

    double precision :: cc, csq, Clag
    double precision :: rho, u, v, p, rhoe_g, h_g
    double precision :: gam_g, game

    double precision :: drho, dptot, drhoe_g
    double precision :: dge, dtau
    double precision :: dup, dvp, dptotp
    double precision :: dum, dvm, dptotm

    double precision :: rho_ref, u_ref, v_ref, p_ref, rhoe_g_ref, h_g_ref
    double precision :: tau_ref

    double precision :: cc_ref, csq_ref, Clag_ref, gam_g_ref, game_ref, gfactor
    double precision :: cc_ev, csq_ev, Clag_ev, rho_ev, p_ev, h_g_ev, tau_ev

    double precision :: alpham, alphap, alpha0r, alpha0e_g
    double precision :: sourcr,sourcp,source,courn,eta,dlogatmp

    double precision :: tau_s

    double precision, allocatable :: Ip(:,:,:,:,:)
    double precision, allocatable :: Im(:,:,:,:,:)

    double precision, allocatable :: Ip_src(:,:,:,:,:)
    double precision, allocatable :: Im_src(:,:,:,:,:)

    ! gamma_c/1 on the interfaces
    double precision, allocatable :: Ip_gc(:,:,:,:,:)
    double precision, allocatable :: Im_gc(:,:,:,:,:)

    type (eos_t) :: eos_state

    if (ppm_type == 0) then
       print *,'Oops -- shouldnt be in trace_ppm with ppm_type = 0'
       call bl_error("Error:: ppm_2d.f90 :: trace_ppm")
    end if

    dtdx = dt/dx
    dtdy = dt/dy

    ! indices: (x, y, dimension, wave, variable)
    allocate(Ip(ilo1-1:ihi1+1,ilo2-1:ihi2+1,2,3,QVAR))
    allocate(Im(ilo1-1:ihi1+1,ilo2-1:ihi2+1,2,3,QVAR))

    if (ppm_trace_sources == 1) then
       allocate(Ip_src(ilo1-1:ihi1+1,ilo2-1:ihi2+1,2,3,QVAR))
       allocate(Im_src(ilo1-1:ihi1+1,ilo2-1:ihi2+1,2,3,QVAR))
    endif

    allocate(Ip_gc(ilo1-1:ihi1+1,ilo2-1:ihi2+1,2,3,1))
    allocate(Im_gc(ilo1-1:ihi1+1,ilo2-1:ihi2+1,2,3,1))

    hdt = HALF * dt

    !=========================================================================
    ! PPM CODE
    !=========================================================================

    ! This does the characteristic tracing to build the interface
    ! states using the normal predictor only (no transverse terms).
    !
    ! We first fill the Im and Ip arrays -- these are the averages of
    ! the various primitive state variables under the parabolic
    ! interpolant over the region swept out by one of the 3 different
    ! characteristic waves.
    !
    ! Im is integrating to the left interface of the current zone
    ! (which will be used to build the right state at that interface)
    ! and Ip is integrating to the right interface of the current zone
    ! (which will be used to build the left state at that interface).
    !
    ! The indices are: Ip(i, j, dim, wave, var)
    !
    ! The choice of reference state is designed to minimize the
    ! effects of the characteristic projection.  We subtract the I's
    ! off of the reference state, project the quantity such that it is
    ! in terms of the characteristic varaibles, and then add all the
    ! jumps that are moving toward the interface to the reference
    ! state to get the full state on that interface.

    ! Compute Ip and Im -- this does the parabolic reconstruction,
    ! limiting, and returns the integral of each profile under each
    ! wave to each interface
    do n=1,QVAR
       call ppm(q(:,:,n),qd_l1,qd_l2,qd_h1,qd_h2, &
                q(:,:,QU:QV),c,qd_l1,qd_l2,qd_h1,qd_h2, &
                flatn, &
                Ip(:,:,:,:,n),Im(:,:,:,:,n), &
                ilo1,ilo2,ihi1,ihi2,dx,dy,dt)
    end do

    ! temperature-based PPM -- if desired, take the Ip(T)/Im(T)
    ! constructed above and use the EOS to overwrite Ip(p)/Im(p)
    if (ppm_temp_fix == 1) then
       call build(eos_state)
       do j = ilo2-1, ihi2+1
          do i = ilo1-1, ihi1+1
             do idim = 1, 2
                do iwave = 1, 3
                   eos_state%rho         = Ip(i,j,idim,iwave,QRHO)
                   eos_state%T           = Ip(i,j,idim,iwave,QTEMP)
                   eos_state%massfrac(:) = Ip(i,j,idim,iwave,QFS:QFS-1+nspec)
                   eos_state%aux         = Ip(i,j,idim,iwave,QFX:QFX-1+naux)

                   call eos_rt(eos_state)

                   Ip(i,j,idim,iwave,QPRES) = eos_state%p
                   Ip(i,j,idim,iwave,QREINT) = Ip(i,j,idim,iwave,QRHO)*eos_state%e
                   Ip_gc(i,j,idim,iwave,1) = eos_state%gam1

                   eos_state%rho         = Im(i,j,idim,iwave,QRHO)
                   eos_state%T           = Im(i,j,idim,iwave,QTEMP)
                   eos_state%massfrac(:) = Im(i,j,idim,iwave,QFS:QFS-1+nspec)
                   eos_state%aux         = Im(i,j,idim,iwave,QFX:QFX-1+naux)

                   call eos_rt(eos_state)

                   Im(i,j,idim,iwave,QPRES) = eos_state%p
                   Im(i,j,idim,iwave,QREINT) = Im(i,j,idim,iwave,QRHO)*eos_state%e
                   Im_gc(i,j,idim,iwave,1) = eos_state%gam1
                enddo
             enddo
          enddo
       enddo
       call destroy(eos_state)
    endif

    ! get an edge-based gam1 here if we didn't get it from the EOS
    ! call above (for ppm_temp_fix = 1)
    if (ppm_temp_fix /= 1) then
          call ppm(gamc(:,:),gc_l1,gc_l2,gc_h1,gc_h2, &
                   q(:,:,QU:QV),c,qd_l1,qd_l2,qd_h1,qd_h2, &
                   flatn, &
                   Ip_gc(:,:,:,:,1),Im_gc(:,:,:,:,1), &
                   ilo1,ilo2,ihi1,ihi2,dx,dy,dt)
       endif

    if (ppm_trace_sources == 1) then
       do n = 1, QVAR
          call ppm(srcQ(:,:,n),src_l1,src_l2,src_h1,src_h2, &
                   q(:,:,QU:QV),c,qd_l1,qd_l2,qd_h1,qd_h2, &
                   flatn, &
                   Ip_src(:,:,:,:,n),Im_src(:,:,:,:,n), &
                   ilo1,ilo2,ihi1,ihi2,dx,dy,dt)
       enddo
    endif

    !-------------------------------------------------------------------------
    ! x-direction
    !-------------------------------------------------------------------------

    ! Trace to left and right edges using upwind PPM
    do j = ilo2-1, ihi2+1
       do i = ilo1-1, ihi1+1

          cc = c(i,j)
          csq = cc**2

          rho = q(i,j,QRHO)
          u = q(i,j,QU)
          v = q(i,j,QV)

          p = q(i,j,QPRES)
          rhoe_g = q(i,j,QREINT)
          h_g = ( (p+rhoe_g)/rho )/csq

          Clag = rho*cc

          gam_g = gamc(i,j)
          game = q(i,j,QGAME)

          !-------------------------------------------------------------------
          ! plus state on face i
          !-------------------------------------------------------------------

          ! set the reference state
          ! this will be the fastest moving state to the left --
          ! this is the method that Miller & Colella and Colella &
          ! Woodward use
          rho_ref  = Im(i,j,1,1,QRHO)
          u_ref    = Im(i,j,1,1,QU)

          p_ref    = Im(i,j,1,1,QPRES)
          rhoe_g_ref = Im(i,j,1,1,QREINT)

          tau_ref  = ONE/Im(i,j,1,1,QRHO)

          gam_g_ref  = Im_gc(i,j,1,1,1)
          game_ref  = Im(i,j,1,1,QGAME)

          rho_ref = max(rho_ref,small_dens)
          p_ref = max(p_ref,small_pres)

          ! for tracing (optionally)
          cc_ref = sqrt(gam_g_ref*p_ref/rho_ref)
          csq_ref = cc_ref**2
          Clag_ref = rho_ref*cc_ref
          h_g_ref = ( (p_ref + rhoe_g_ref)/rho_ref )/csq_ref

          ! *m are the jumps carried by u-c
          ! *p are the jumps carried by u+c

          dum    = u_ref    - Im(i,j,1,1,QU)
          dptotm    = p_ref    - Im(i,j,1,1,QPRES)

          drho  = rho_ref  - Im(i,j,1,2,QRHO)
          dptot    = p_ref    - Im(i,j,1,2,QPRES)
          drhoe_g = rhoe_g_ref - Im(i,j,1,2,QREINT)
          dtau  = tau_ref  - ONE/Im(i,j,1,2,QRHO)

          dup    = u_ref    - Im(i,j,1,3,QU)
          dptotp    = p_ref    - Im(i,j,1,3,QPRES)

          ! if we are doing source term tracing, then we add the force to
          ! the velocity here, otherwise we will deal with this in the
          ! trans_X routines
          if (ppm_trace_sources == 1) then
             dum = dum - hdt*Im_src(i,j,1,1,QU)
             dup = dup - hdt*Im_src(i,j,1,3,QU)
          endif

          ! optionally use the reference state in evaluating the
          ! eigenvectors
          if (ppm_reference_eigenvectors == 0) then
             rho_ev  = rho
             cc_ev   = cc
             csq_ev  = csq
             Clag_ev = Clag
             h_g_ev = h_g
             p_ev    = p
             tau_ev  = ONE/rho
          else
             rho_ev  = rho_ref
             cc_ev   = cc_ref
             csq_ev  = csq_ref
             Clag_ev = Clag_ref
             h_g_ev = h_g_ref
             p_ev    = p_ref
             tau_ev  = tau_ref
          endif

          if (ppm_predict_gammae == 0) then

             ! (rho, u, p, (rho e)) eigensystem

             ! these are analogous to the beta's from the original
             ! PPM paper (except we work with rho instead of tau).
             ! This is simply (l . dq), where dq = qref - I(q)

             alpham = HALF*(dptotm/(rho_ev*cc_ev) - dum)*rho_ev/cc_ev
             alphap = HALF*(dptotp/(rho_ev*cc_ev) + dup)*rho_ev/cc_ev
             alpha0r = drho - dptot/csq_ev
             alpha0e_g = drhoe_g - dptot*h_g_ev  ! note h_g has a 1/c**2 in it

          else

             ! (tau, u, p, game) eigensystem

             ! this is the way things were done in the original PPM
             ! paper -- here we work with tau in the characteristic
             ! system.

             alpham = HALF*( dum - dptotm/Clag_ev)/Clag_ev
             alphap = HALF*(-dup - dptotp/Clag_ev)/Clag_ev
             alpha0r = dtau + dptot/Clag_ev**2

             dge = game_ref - Im(i,j,1,2,QGAME)
             gfactor = (game - ONE)*(game - gam_g)  ! why not reference ev here?
             alpha0e_g = gfactor*dptot/(tau_ev*Clag_ev**2) + dge

          endif ! which tracing method

          if (u-cc > ZERO) then
             alpham = ZERO
          else if (u-cc < ZERO) then
             alpham = -alpham
          else
             alpham = -HALF*alpham
          endif

          if (u+cc > ZERO) then
             alphap = ZERO
          else if (u+cc < ZERO) then
             alphap = -alphap
          else
             alphap = -HALF*alphap
          endif

          if (u > ZERO) then
             alpha0r = ZERO
             alpha0e_g = ZERO
          else if (u < ZERO) then
             alpha0r = -alpha0r
             alpha0e_g = -alpha0e_g
          else
             alpha0r = -HALF*alpha0r
             alpha0e_g = -HALF*alpha0e_g
          endif

          ! the final interface states are just
          ! q_s = q_ref - sum (l . dq) r
          if (i >= ilo1) then

             if (ppm_predict_gammae == 0) then
                qxp(i,j,QRHO)   = rho_ref + alphap + alpham + alpha0r
                qxp(i,j,QU)     = u_ref + (alphap - alpham)*cc_ev/rho_ev
                qxp(i,j,QREINT) = rhoe_g_ref + (alphap + alpham)*h_g_ev*csq_ev + alpha0e_g
                qxp(i,j,QPRES)  = p_ref + (alphap + alpham)*csq_ev
             else
                tau_s = tau_ref + alphap + alpham + alpha0r
                qxp(i,j,QRHO)   = ONE/tau_s

                qxp(i,j,QU)     = u_ref + (alpham - alphap)*Clag_ev
                qxp(i,j,QPRES)  = p_ref - (alphap + alpham)*Clag_ev**2

                qxp(i,j,QGAME) = game_ref + gfactor*(alpham + alphap)/tau_ev + alpha0e_g
                qxp(i,j,QREINT) = qxp(i,j,QPRES )/(qxp(i,j,QGAME) - ONE)
             end if

             ! enforce small_*
             qxp(i,j,QRHO) = max(small_dens,qxp(i,j,QRHO))
             qxp(i,j,QPRES) = max(qxp(i,j,QPRES), small_pres)

             ! transverse velocity -- there is no projection here, so
             ! we don't need a reference state.  We only care about
             ! the state traced under the middle wave

             ! Recall that I already takes the limit of the parabola
             ! in the event that the wave is not moving toward the
             ! interface
             qxp(i,j,QV) = Im(i,j,1,2,QV)

             if (ppm_trace_sources == 1) then
                qxp(i,j,QV)  = qxp(i,j,QV)  + hdt*Im_src(i,j,1,2,QV)
             endif

          endif

          !-------------------------------------------------------------------
          ! minus state on face i+1
          !-------------------------------------------------------------------

          ! set the reference state
          ! this will be the fastest moving state to the right
          rho_ref  = Ip(i,j,1,3,QRHO)
          u_ref    = Ip(i,j,1,3,QU)

          p_ref    = Ip(i,j,1,3,QPRES)
          rhoe_g_ref = Ip(i,j,1,3,QREINT)

          tau_ref  = ONE/Ip(i,j,1,3,QRHO)

          gam_g_ref    = Ip_gc(i,j,1,3,1)
          game_ref    = Ip(i,j,1,3,QGAME)

          rho_ref = max(rho_ref,small_dens)
          p_ref = max(p_ref,small_pres)

          ! for tracing (optionally)
          cc_ref = sqrt(gam_g_ref*p_ref/rho_ref)
          csq_ref = cc_ref**2
          Clag_ref = rho_ref*cc_ref
          h_g_ref = ( (p_ref + rhoe_g_ref)/rho_ref )/csq_ref

          ! *m are the jumps carried by u-c
          ! *p are the jumps carried by u+c

          dum    = u_ref    - Ip(i,j,1,1,QU)
          dptotm    = p_ref    - Ip(i,j,1,1,QPRES)

          drho  = rho_ref  - Ip(i,j,1,2,QRHO)
          dptot    = p_ref    - Ip(i,j,1,2,QPRES)
          drhoe_g = rhoe_g_ref - Ip(i,j,1,2,QREINT)
          dtau  = tau_ref  - ONE/Ip(i,j,1,2,QRHO)

          dup    = u_ref    - Ip(i,j,1,3,QU)
          dptotp    = p_ref    - Ip(i,j,1,3,QPRES)

          ! if we are doing source term tracing, then we add the force to
          ! the velocity here, otherwise we will deal with this in the
          ! trans_X routines
          if (ppm_trace_sources == 1) then
             dum = dum - hdt*Ip_src(i,j,1,1,QU)
             dup = dup - hdt*Ip_src(i,j,1,3,QU)
          endif

          ! optionally use the reference state in evaluating the
          ! eigenvectors
          if (ppm_reference_eigenvectors == 0) then
             rho_ev  = rho
             cc_ev   = cc
             csq_ev  = csq
             Clag_ev = Clag
             h_g_ev = h_g
             p_ev    = p
             tau_ev  = ONE/rho
          else
             rho_ev  = rho_ref
             cc_ev   = cc_ref
             csq_ev  = csq_ref
             Clag_ev = Clag_ref
             h_g_ev = h_g_ref
             p_ev    = p_ref
             tau_ev  = tau_ref
          endif

          if (ppm_predict_gammae == 0) then

             ! (rho, u, p, (rho e)) eigensystem

             ! these are analogous to the beta's from the original
             ! PPM paper (except we work with rho instead of tau).
             ! This is simply (l . dq), where dq = qref - I(q)
             alpham = HALF*(dptotm/(rho_ev*cc_ev) - dum)*rho_ev/cc_ev
             alphap = HALF*(dptotp/(rho_ev*cc_ev) + dup)*rho_ev/cc_ev
             alpha0r = drho - dptot/csq_ev
             alpha0e_g = drhoe_g - dptot*h_g_ev  ! h_g has a 1/c**2 in it

          else

             ! (tau, u, p, game) eigensystem

             ! this is the way things were done in the original PPM
             ! paper -- here we work with tau in the characteristic
             ! system.

             alpham = HALF*( dum - dptotm/Clag_ev)/Clag_ev
             alphap = HALF*(-dup - dptotp/Clag_ev)/Clag_ev
             alpha0r = dtau + dptot/Clag_ev**2

             dge = game_ref - Ip(i,j,1,2,QGAME)
             gfactor = (game - ONE)*(game - gam_g)
             alpha0e_g = gfactor*dptot/(tau_ev*Clag_ev**2) + dge

          endif

          if (u-cc > ZERO) then
             alpham = -alpham
          else if (u-cc < ZERO) then
             alpham = ZERO
          else
             alpham = -HALF*alpham
          endif

          if (u+cc > ZERO) then
             alphap = -alphap
          else if (u+cc < ZERO) then
             alphap = ZERO
          else
             alphap = -HALF*alphap
          endif

          if (u > ZERO) then
             alpha0r = -alpha0r
             alpha0e_g = -alpha0e_g
          else if (u < ZERO) then
             alpha0r = ZERO
             alpha0e_g = ZERO
          else
             alpha0r = -HALF*alpha0r
             alpha0e_g = -HALF*alpha0e_g
          endif

          ! the final interface states are just
          ! q_s = q_ref - sum (l . dq) r
          if (i <= ihi1) then

             if (ppm_predict_gammae == 0) then

                qxm(i+1,j,QRHO)   = rho_ref + alphap + alpham + alpha0r
                qxm(i+1,j,QU)     = u_ref + (alphap - alpham)*cc_ev/rho_ev
                qxm(i+1,j,QREINT) = rhoe_g_ref + (alphap + alpham)*h_g_ev*csq_ev + alpha0e_g
                qxm(i+1,j,QPRES)  = p_ref + (alphap + alpham)*csq_ev

             else
                tau_s = tau_ref + alphap + alpham + alpha0r
                qxm(i+1,j,QRHO)   = ONE/tau_s

                qxm(i+1,j,QU)     = u_ref + (alpham - alphap)*Clag_ev
                qxm(i+1,j,QPRES)  = p_ref - (alphap + alpham)*Clag_ev**2

                qxm(i+1,j,QGAME) = game_ref + gfactor*(alpham + alphap)/tau_ev + alpha0e_g
                qxm(i+1,j,QREINT) = qxm(i+1,j,QPRES )/(qxm(i+1,j,QGAME) - ONE)

             end if

             ! enforce small_*
             qxm(i+1,j,QRHO) = max(qxm(i+1,j,QRHO),small_dens)
             qxm(i+1,j,QPRES) = max(qxm(i+1,j,QPRES), small_pres)

             ! transverse velocity -- there is no projection here, so
             ! we don't need a reference state.  We only care about
             ! the state traced under the middle wave
             qxm(i+1,j,QV) = Ip(i,j,1,2,QV)

             if (ppm_trace_sources == 1) then
                qxm(i+1,j,QV)  = qxm(i+1,j,QV)  + hdt*Ip_src(i,j,1,2,QV)
             endif

          endif

          !-------------------------------------------------------------------
          ! geometry source terms
          !-------------------------------------------------------------------

          if (dloga(i,j) /= 0) then
             courn = dtdx*(cc+abs(u))
             eta = (ONE-courn)/(cc*dt*abs(dloga(i,j)))
             dlogatmp = min(eta,ONE)*dloga(i,j)
             sourcr = -HALF*dt*rho*dlogatmp*u
             sourcp = sourcr*csq
             source = sourcp*h_g

             if (i <= ihi1) then
                qxm(i+1,j,QRHO) = qxm(i+1,j,QRHO) + sourcr
                qxm(i+1,j,QRHO) = max(qxm(i+1,j,QRHO),small_dens)
                qxm(i+1,j,QPRES) = qxm(i+1,j,QPRES) + sourcp
                qxm(i+1,j,QREINT) = qxm(i+1,j,QREINT) + source
             end if

             if (i >= ilo1) then
                qxp(i,j,QRHO) = qxp(i,j,QRHO) + sourcr
                qxp(i,j,QRHO) = max(qxp(i,j,QRHO),small_dens)
                qxp(i,j,QPRES) = qxp(i,j,QPRES) + sourcp
                qxp(i,j,QREINT) = qxp(i,j,QREINT) + source
             end if

          endif

       end do
    end do

    !-------------------------------------------------------------------------
    ! Now do the passively advected quantities
    !-------------------------------------------------------------------------

    ! We do all passively advected quantities in one loop
    do ipassive = 1, npassive
       n = qpass_map(ipassive)
       do j = ilo2-1, ihi2+1

          ! plus state on face i
          do i = ilo1, ihi1+1
             u = q(i,j,QU)

             ! We have
             !
             ! q_l = q_ref - Proj{(q_ref - I)}
             !
             ! and Proj{} represents the characteristic projection.
             ! But for these, there is only 1-wave that matters, the u
             ! wave, so no projection is needed.  Since we are not
             ! projecting, the reference state doesn't matter.

             if (u > ZERO) then
                qxp(i,j,n) = q(i,j,n)    ! we might want to change this to
                                         ! the limit of the parabola
             else if (u < ZERO) then
                qxp(i,j,n) = Im(i,j,1,2,n)
             else
                qxp(i,j,n) = q(i,j,n) + HALF*(Im(i,j,1,2,n) - q(i,j,n))
             endif
          enddo

          ! minus state on face i+1
          do i = ilo1-1, ihi1
             u = q(i,j,QU)

             if (u > ZERO) then
                qxm(i+1,j,n) = Ip(i,j,1,2,n)
             else if (u < ZERO) then
                qxm(i+1,j,n) = q(i,j,n)
             else
                qxm(i+1,j,n) = q(i,j,n) + HALF*(Ip(i,j,1,2,n) - q(i,j,n))
             endif
          enddo

       enddo
    enddo

    !-------------------------------------------------------------------------
    ! y-direction
    !-------------------------------------------------------------------------

    ! Trace to bottom and top edges using upwind PPM
    do j = ilo2-1, ihi2+1
       do i = ilo1-1, ihi1+1

          cc = c(i,j)
          csq = cc**2

          rho = q(i,j,QRHO)
          u = q(i,j,QU)
          v = q(i,j,QV)

          p = q(i,j,QPRES)
          rhoe_g = q(i,j,QREINT)
          h_g = ( (p+rhoe_g)/rho )/csq

          Clag = rho*cc

          gam_g = gamc(i,j)
          game = q(i,j,QGAME)

          !-------------------------------------------------------------------
          ! plus state on face j
          !-------------------------------------------------------------------

          ! set the reference state
          ! this will be the fastest moving state to the left
          rho_ref  = Im(i,j,2,1,QRHO)
          v_ref    = Im(i,j,2,1,QV)

          p_ref    = Im(i,j,2,1,QPRES)
          rhoe_g_ref = Im(i,j,2,1,QREINT)

          tau_ref  = ONE/Im(i,j,2,1,QRHO)

          gam_g_ref  = Im_gc(i,j,2,1,1)
          game_ref = Im(i,j,2,1,QGAME)

          rho_ref = max(rho_ref,small_dens)
          p_ref = max(p_ref,small_pres)

          ! for tracing (optionally)
          cc_ref = sqrt(gam_g_ref*p_ref/rho_ref)
          csq_ref = cc_ref**2
          Clag_ref = rho_ref*cc_ref
          h_g_ref = ( (rhoe_g_ref+p_ref)/rho_ref )/csq_ref

          ! *m are the jumps carried by v-c
          ! *p are the jumps carried by v+c

          dvm    = v_ref    - Im(i,j,2,1,QV)
          dptotm    = p_ref    - Im(i,j,2,1,QPRES)

          drho  = rho_ref  - Im(i,j,2,2,QRHO)
          dptot    = p_ref    - Im(i,j,2,2,QPRES)
          drhoe_g = rhoe_g_ref - Im(i,j,2,2,QREINT)
          dtau  = tau_ref  - ONE/Im(i,j,2,2,QRHO)

          dvp    = v_ref    - Im(i,j,2,3,QV)
          dptotp    = p_ref    - Im(i,j,2,3,QPRES)

          ! if we are doing source term tracing, then we add the force to
          ! the velocity here, otherwise we will deal with this in the
          ! trans_X routines
          if (ppm_trace_sources == 1) then
             dvm = dvm - hdt*Im_src(i,j,2,1,QV)
             dvp = dvp - hdt*Im_src(i,j,2,3,QV)
          endif

          ! optionally use the reference state in evaluating the
          ! eigenvectors
          if (ppm_reference_eigenvectors == 0) then
             rho_ev  = rho
             cc_ev   = cc
             csq_ev  = csq
             Clag_ev = Clag
             h_g_ev = h_g
             p_ev    = p
          else
             rho_ev  = rho_ref
             cc_ev   = cc_ref
             csq_ev  = csq_ref
             Clag_ev = Clag_ref
             h_g_ev = h_g_ref
             p_ev    = p_ref
          endif

          if (ppm_predict_gammae == 0) then

             ! (rho, u, p, (rho e)) eigensystem

             ! these are analogous to the beta's from the original PPM
             ! paper (except we work with rho instead of tau).  This
             ! is simply (l . dq), where dq = qref - I(q)
             alpham = HALF*(dptotm/(rho_ev*cc_ev) - dvm)*rho_ev/cc_ev
             alphap = HALF*(dptotp/(rho_ev*cc_ev) + dvp)*rho_ev/cc_ev
             alpha0r = drho - dptot/csq_ev
             alpha0e_g = drhoe_g - dptot*h_g_ev  ! h_g has 1/c**2 in it

          else

             ! (tau, u, p, game) eigensystem

             ! this is the way things were done in the original PPM
             ! paper -- here we work with tau in the characteristic
             ! system.

             alpham = HALF*( dvm - dptotm/Clag_ev)/Clag_ev
             alphap = HALF*(-dvp - dptotp/Clag_ev)/Clag_ev
             alpha0r = dtau + dptot/Clag_ev**2

             dge = game_ref - Im(i,j,2,2,QGAME)
             gfactor = (game - ONE)*(game - gam_g)
             alpha0e_g = gfactor*dptot/(tau_ev*Clag_ev**2) + dge

          endif

          if (v-cc > ZERO) then
             alpham = ZERO
          else if (v-cc < ZERO) then
             alpham = -alpham
          else
             alpham = -HALF*alpham
          endif

          if (v+cc > ZERO) then
             alphap = ZERO
          else if (v+cc < ZERO) then
             alphap = -alphap
          else
             alphap = -HALF*alphap
          endif

          if (v > ZERO) then
             alpha0r = ZERO
             alpha0e_g = ZERO
          else if (v < ZERO) then
             alpha0r = -alpha0r
             alpha0e_g = -alpha0e_g
          else
             alpha0r = -HALF*alpha0r
             alpha0e_g = -HALF*alpha0e_g
          endif

          ! the final interface states are just
          ! q_s = q_ref - sum (l . dq) r
          if (j >= ilo2) then
             if (ppm_predict_gammae == 0) then
                qyp(i,j,QRHO)   = rho_ref + alphap + alpham + alpha0r
                qyp(i,j,QV)     = v_ref + (alphap - alpham)*cc_ev/rho_ev
                qyp(i,j,QREINT) = rhoe_g_ref + (alphap + alpham)*h_g_ev*csq_ev + alpha0e_g
                qyp(i,j,QPRES)  = p_ref + (alphap + alpham)*csq_ev

             else
                tau_s = tau_ref + alphap + alpham + alpha0r
                qyp(i,j,QRHO)   = ONE/tau_s

                qyp(i,j,QV)     = v_ref + (alpham - alphap)*Clag_ev
                qyp(i,j,QPRES)  = p_ref - (alphap + alpham)*Clag_ev**2

                qyp(i,j,QGAME) = game_ref + gfactor*(alpham + alphap)/tau_ev + alpha0e_g
                qyp(i,j,QREINT) = qyp(i,j,QPRES )/(qyp(i,j,QGAME) - ONE)

             end if

             ! enforce small_*
             qyp(i,j,QRHO) = max(small_dens, qyp(i,j,QRHO))
             qyp(i,j,QPRES) = max(qyp(i,j,QPRES), small_pres)

             ! transverse velocity -- there is no projection here, so
             ! we don't need a reference state.  We only care about
             ! the state traced under the middle wave
             qyp(i,j,QU)     = Im(i,j,2,2,QU)

             if (ppm_trace_sources == 1) then
                qyp(i,j,QU)  = qyp(i,j,QU)  + hdt*Im_src(i,j,2,2,QU)
             endif

          endif

          !-------------------------------------------------------------------
          ! minus state on face j+1
          !-------------------------------------------------------------------

          ! set the reference state
          ! this will be the fastest moving state to the right
          rho_ref  = Ip(i,j,2,3,QRHO)
          v_ref    = Ip(i,j,2,3,QV)

          p_ref    = Ip(i,j,2,3,QPRES)
          rhoe_g_ref = Ip(i,j,2,3,QREINT)

          tau_ref  = ONE/Ip(i,j,2,3,QRHO)

          gam_g_ref    = Ip_gc(i,j,2,3,1)
          game_ref    = Ip(i,j,2,3,QGAME)

          rho_ref = max(rho_ref,small_dens)
          p_ref = max(p_ref,small_pres)

          ! for tracing (optionally)
          cc_ref = sqrt(gam_g_ref*p_ref/rho_ref)
          csq_ref = cc_ref**2
          Clag_ref = rho_ref*cc_ref
          h_g_ref = ( (rhoe_g_ref+p_ref)/rho_ref )/csq_ref

          ! *m are the jumps carried by v-c
          ! *p are the jumps carried by v+c

          dvm    = v_ref    - Ip(i,j,2,1,QV)
          dptotm    = p_ref    - Ip(i,j,2,1,QPRES)

          drho  = rho_ref  - Ip(i,j,2,2,QRHO)
          dptot    = p_ref    - Ip(i,j,2,2,QPRES)
          drhoe_g = rhoe_g_ref - Ip(i,j,2,2,QREINT)
          dtau  = tau_ref  - ONE/Ip(i,j,2,2,QRHO)

          dvp    = v_ref    - Ip(i,j,2,3,QV)
          dptotp    = p_ref    - Ip(i,j,2,3,QPRES)

          ! if we are doing source term tracing, then we add the force to
          ! the velocity here, otherwise we will deal with this in the
          ! trans_X routines
          if (ppm_trace_sources == 1) then
             dvm = dvm - hdt*Ip_src(i,j,2,1,QV)
             dvp = dvp - hdt*Ip_src(i,j,2,3,QV)
          endif

          ! optionally use the reference state in evaluating the
          ! eigenvectors
          if (ppm_reference_eigenvectors == 0) then
             rho_ev  = rho
             cc_ev   = cc
             csq_ev  = csq
             Clag_ev = Clag
             h_g_ev = h_g
             p_ev    = p
          else
             rho_ev  = rho_ref
             cc_ev   = cc_ref
             csq_ev  = csq_ref
             Clag_ev = Clag_ref
             h_g_ev = h_g_ref
             p_ev    = p_ref
          endif

          if (ppm_predict_gammae == 0) then

             ! (rho, u, p, (rho e) eigensystem

             ! these are analogous to the beta's from the original PPM
             ! paper.  This is simply (l . dq), where dq = qref - I(q)
             alpham = HALF*(dptotm/(rho_ev*cc_ev) - dvm)*rho_ev/cc_ev
             alphap = HALF*(dptotp/(rho_ev*cc_ev) + dvp)*rho_ev/cc_ev
             alpha0r = drho - dptot/csq_ev
             alpha0e_g = drhoe_g - dptot*h_g_ev

          else

             ! (tau, u, p, game) eigensystem

             ! this is the way things were done in the original PPM
             ! paper -- here we work with tau in the characteristic
             ! system.

             alpham = HALF*( dvm - dptotm/Clag_ev)/Clag_ev
             alphap = HALF*(-dvp - dptotp/Clag_ev)/Clag_ev
             alpha0r = dtau + dptot/Clag_ev**2

             dge = game_ref - Ip(i,j,2,2,QGAME)
             gfactor = (game - ONE)*(game - gam_g)
             alpha0e_g = gfactor*dptot/(tau_ev*Clag_ev**2) + dge

          endif

          if (v-cc > ZERO) then
             alpham = -alpham
          else if (v-cc < ZERO) then
             alpham = ZERO
          else
             alpham = -HALF*alpham
          endif

          if (v+cc > ZERO) then
             alphap = -alphap
          else if (v+cc < ZERO) then
             alphap = ZERO
          else
             alphap = -HALF*alphap
          endif

          if (v > ZERO) then
             alpha0r = -alpha0r
             alpha0e_g = -alpha0e_g
          else if (v < ZERO) then
             alpha0r = ZERO
             alpha0e_g = ZERO
          else
             alpha0r = -HALF*alpha0r
             alpha0e_g = -HALF*alpha0e_g
          endif

          ! the final interface states are just
          ! q_s = q_ref - sum (l . dq) r
          if (j <= ihi2) then
             if (ppm_predict_gammae == 0) then
                qym(i,j+1,QRHO)   = rho_ref + alphap + alpham + alpha0r
                qym(i,j+1,QV)     = v_ref + (alphap - alpham)*cc_ev/rho_ev
                qym(i,j+1,QREINT) = rhoe_g_ref + (alphap + alpham)*h_g_ev*csq_ev + alpha0e_g
                qym(i,j+1,QPRES)  = p_ref + (alphap + alpham)*csq_ev
             else
                tau_s = tau_ref + alphap + alpham + alpha0r
                qym(i,j+1,QRHO)   = ONE/tau_s

                qym(i,j+1,QV)     = v_ref + (alpham - alphap)*Clag_ev
                qym(i,j+1,QPRES)  = p_ref - (alphap + alpham)*Clag_ev**2

                qym(i,j+1,QGAME) = game_ref + gfactor*(alpham + alphap)/tau_ev + alpha0e_g
                qym(i,j+1,QREINT) = qym(i,j+1,QPRES )/(qym(i,j+1,QGAME) - ONE)

             end if

             ! enforce small_*
             qym(i,j+1,QRHO) = max(small_dens, qym(i,j+1,QRHO))
             qym(i,j+1,QPRES) = max(qym(i,j+1,QPRES), small_pres)

             ! transverse velocity -- there is no projection here, so
             ! we don't need a reference state.  We only care about
             ! the state traced under the middle wave
             qym(i,j+1,QU)    = Ip(i,j,2,2,QU)

             if (ppm_trace_sources == 1) then
                qym(i,j+1,QU)  = qym(i,j+1,QU)  + hdt*Ip_src(i,j,2,2,QU)
             endif

          endif

       end do
    end do

    !-------------------------------------------------------------------------
    ! Now do the passively advected quantities
    !-------------------------------------------------------------------------

    ! do all of the passively advected quantities in one loop
    do ipassive = 1, npassive
       n = qpass_map(ipassive)
       do i = ilo1-1, ihi1+1

          ! plus state on face j
          do j = ilo2, ihi2+1
             v = q(i,j,QV)

             if (v > ZERO) then
                qyp(i,j,n) = q(i,j,n)
             else if (v < ZERO) then
                qyp(i,j,n) = Im(i,j,2,2,n)
             else
                qyp(i,j,n) = q(i,j,n) + HALF*(Im(i,j,2,2,n) - q(i,j,n))
             endif
          enddo

          ! minus state on face j+1
          do j = ilo2-1, ihi2
             v = q(i,j,QV)

             if (v > ZERO) then
                qym(i,j+1,n) = Ip(i,j,2,2,n)
             else if (v < ZERO) then
                qym(i,j+1,n) = q(i,j,n)
             else
                qym(i,j+1,n) = q(i,j,n) + HALF*(Ip(i,j,2,2,n) - q(i,j,n))
             endif
          enddo

       enddo
    enddo

    deallocate(Ip,Im)
    deallocate(Ip_gc,Im_gc)
    if (ppm_trace_sources == 1) then
       deallocate(Ip_src,Im_src)
    endif

  end subroutine trace_ppm

end module trace_ppm_module

max_stars_repo_name: wjge/PeleC

max_issues_repo_name: marient/PeleC

max_forks_repo_name: marient/PeleC


text:
      subroutine bs2(n0,n1,xi,yi,ui,vi,n2,n3,xj,yj,gj,sj,mp,lm)
      implicit real*8(a-h,o-z)

      include 'memory.f'

      real tim(2),tfmm(9)

      complex*16  a(mpmax,nbne),ao(mpmax,nbne), b(mpmax,nbne,lmax)
      complex*16 ss(mpmax,mpmax),sr(mpmax,mpmax),rr(mpmax,mpmax)
      complex*16 aa(mpmax),bb(mpmax),c(mpmax),d(mpmax)
      complex*16 invt,invt2,t,ziic,zjjc

      dimension xi(npmax),yi(npmax),ui(npmax),vi(npmax)
      dimension xj(npmax),yj(npmax),gj(npmax),sj(npmax)

      dimension xdi(npmax),ydi(npmax),udi(npmax),vdi(npmax)
      dimension xdj(npmax),ydj(npmax),gdj(npmax),sdj(npmax)

      dimension ndi(npmax,2),ndbi(npmax),ndj(npmax,2),ndbj(npmax)
      dimension nbi(npmax),nfi(nbmax),nei(nbmax),nfio(nbmax),neio(nbmax)
      dimension nbj(npmax),nfj(nbmax),nej(nbmax),nfjo(nbmax),nejo(nbmax)
      dimension nii(nbmax),njj(nbmax),nc(3),nfn(9)

      call cpu_time(tic)
      do i = 1,9
       tfmm(i) = 0
      end do
      pi = 2*acos(0.0)
      xmax = -1e10
      xmin = 1e10
      ymax = -1e10
      ymin = 1e10

      do i = n0,n1
       xmax = max(xi(i),xmax)
       xmin = min(xi(i),xmin)
       ymax = max(yi(i),ymax)
       ymin = min(yi(i),ymin)
      end do

      do i = n2,n3
       xmax = max(xj(i),xmax)
       xmin = min(xj(i),xmin)
       ymax = max(yj(i),ymax)
       ymin = min(yj(i),ymin)
      end do

      rdx = (xmax-xmin)*(1+1e-5)
      rdy = (ymax-ymin)*(1+1e-5)
      rd = max(rdx,rdy)
      nmax = max(n1-n0+1,n3-n2+1)

      lev = lm
      lml = 2

      call boxdata(n0,n1,xi,yi,rd,lev,nbi,ndbi,ndi,nei,nfi,
     *    lbi,rb,xmin,ymin)
      call boxdata(n2,n3,xj,yj,rd,lev,nbj,ndbj,ndj,nej,nfj,
     *    lbj,rb,xmin,ymin)

      n = 0
      do i = n0,n1
       n = n+1
       xdi(n) = xi(ndbi(n))
       ydi(n) = yi(ndbi(n))
      end do

      n = 0
      do i = n2,n3
       n = n+1
       xdj(n) = xj(ndbj(n))
       ydj(n) = yj(ndbj(n))
       gdj(n) = gj(ndbj(n))
       sdj(n) = sj(ndbj(n))
      end do

      toc = tic
      call cpu_time(tic)
      tfmm(1) = tfmm(1)+tic-toc

c Step 1. S-expansion

      do i = 1,lbj
       call boxc(nfj(i),2,nc)
       xjc = xmin+(nc(1)+.5)*rb
       yjc = ymin+(nc(2)+.5)*rb
       do j = 1,mp
        bb(j) = 0
       end do

       do j = ndj(i,1),ndj(i,2)
        xjjc = xdj(j)-xjc
        yjjc = ydj(j)-yjc
        zjjc = dcmplx(xjjc,yjjc)
        do k = 1,mp
         bb(k) = bb(k)+gdj(j)*zjjc**(k-1)
        end do
       end do

       do j = 1,mp
        b(j,i,lev) = bb(j)
       end do
      end do

      toc = tic
      call cpu_time(tic)
      tfmm(2) = tfmm(2)+tic-toc

c Step 2. S|S-translation

      if(lm.gt.lml)then
       do lev = lm-1,lml,-1
        do i = 1,4**(lev+1)
         nejo(i) = nej(i)
        end do
        do i = 1,lbj
         nfjo(i) = nfj(i)
        end do
        lbjo = lbj
        call boxdata(n2,n3,xj,yj,rd,lev,nbj,ndbj,ndj,nej,nfj,
     *    lbj,rb,xmin,ymin)
        do i = 1,lbj
         call boxc(nfj(i),2,nc)
         xjc = (nc(1)+.5)*rb
         yjc = (nc(2)+.5)*rb
         do j = 1,mp
          bb(j) = dcmplx(0.0,0.0)
         end do
         do j = 1,4
          nfjc = nfj(i)*4+j-1
          do k = 1,lbjo
           if(nfjo(k).eq.nfjc)then
            call boxc(nfjc,2,nc)
            xjco = (nc(1)+.5)*rb/2
            yjco = (nc(2)+.5)*rb/2
            t = dcmplx(xjc-xjco,yjc-yjco)

            ss(1,1) = 1
            nminus = 1
            do l = 1,mp
             c(l) = t**(l-1)
            end do
            do n = 2,mp
             nminus = -nminus
             ss(n,1) = 0
             ss(1,n) = nminus*c(n)
             mminus = nminus
             d(1) = 1
             do m = 2,n-1
              d(m) = real(n-m+1)/(m-1)*d(m-1)
              mminus = -mminus
              ss(n,m) = 0
              ss(m,n) = mminus*c(n-m+1)*d(m)
             end do
             ss(n,n) = ss(n-1,n-1)
            end do

            do l = 1,mp
             do m = 1,mp
              bb(l) = bb(l)+ss(m,l)*b(m,nejo(nfjc+1),lev+1)
             end do
            end do

           end if
          end do
         end do
         do j = 1,mp
          b(j,i,lev) = bb(j)
         end do
        end do
       end do
       lev = lml
      end if
      call boxdata(n0,n1,xi,yi,rd,lev,nbi,ndbi,ndi,nei,nfi,
     *    lbi,rb,xmin,ymin)

      toc = tic
      call cpu_time(tic)
      tfmm(3) = tfmm(3)+tic-toc

c Step 3. S|R-translation

      do i = 1,lbi
       call boxc(nfi(i),2,nc)
       xic = (nc(1)+.5)*rb
       yic = (nc(2)+.5)*rb
       do j = 1,mp
        aa(j) = 0
       end do
       call e3b(nfi(i),nfj,lbj,nii,li,lev)
       do j = 1,li
        jj = nii(j)
        call boxc(jj,2,nc)
        xjc = (nc(1)+.5)*rb
        yjc = (nc(2)+.5)*rb
        t = dcmplx(xic-xjc,yic-yjc)

        invt = 1/t
        sr(1,1) = invt
        invt2 = invt*invt
        nminus = 1
        do k = 1,mp
         c(k) = invt**k
        end do
        do n = 2,mp
         nminus = -nminus
         sr(n,1) = c(n)
         sr(1,n) = nminus*c(n)
         d(1) = c(n)
         mminus = nminus
         do m = 2,n-1
          mminus = -mminus
          d(m) = -(1+real(n-1)/(m-1))*d(m-1)*invt
          sr(n,m) = d(m)
          sr(m,n) = mminus*d(m)
         end do
         sr(n,n) = -2*invt2*(2-1.0/(n-1))*sr(n-1,n-1)
        end do

        do k = 1,mp
         do l = 1,mp
          aa(k) = aa(k)+sr(l,k)*b(l,nej(jj+1),lev)
         end do
        end do

       end do
       do j = 1,mp
        a(j,i) = aa(j)
       end do
      end do

      toc = tic
      call cpu_time(tic)
      tfmm(4) = tfmm(4)+tic-toc

c Step 4. R|R-translation

      if(lm.gt.lml)then
       do lev = lml+1,lm

c R|R-translation

        do i = 1,4**(lev-1)
         neio(i) = nei(i)
        end do
        do i = 1,lbi
         nfio(i) = nfi(i)
         do j = 1,mp
          ao(j,i) = a(j,i)
         end do
        end do
        lbio = lbi
        call boxdata(n0,n1,xi,yi,rd,lev,nbi,ndbi,ndi,nei,nfi,
     *    lbi,rb,xmin,ymin)
        do i = 1,lbi
         call boxc(nfi(i),2,nc)
         xic = (nc(1)+.5)*rb
         yic = (nc(2)+.5)*rb
         nfip = nfi(i)/4
         call boxc(nfip,2,nc)
         xico = (nc(1)+.5)*rb*2
         yico = (nc(2)+.5)*rb*2
         t = dcmplx(xic-xico,yic-yico)

         rr(1,1) = 1
         do j = 1,mp
          c(j) = t**(j-1)
         end do
         do m = 2,mp
          rr(m,1) = c(m)
          rr(1,m) = 0
          d(1) = 1
          do n = 2,m-1
           d(n) = real(m-n+1)/(n-1)*d(n-1)
           rr(n,m) = 0
           rr(m,n) = c(m-n+1)*d(n)
          end do
          rr(m,m) = rr(m-1,m-1)
         end do

         do j = 1,mp
          t = 0
          do k = 1,mp
           t = t+rr(k,j)*ao(k,neio(nfip+1))
          end do
          a(j,i) = t
         end do

        end do
        call boxdata(n2,n3,xj,yj,rd,lev,nbj,ndbj,ndj,nej,nfj,
     *    lbj,rb,xmin,ymin)

        toc = tic
        call cpu_time(tic)
        tfmm(5) = tfmm(5)+tic-toc

c S|R-translation

        do i = 1,lbi
         call boxc(nfi(i),2,nc)
         xic = (nc(1)+.5)*rb
         yic = (nc(2)+.5)*rb
         do j = 1,mp
          aa(j) = a(j,i)
         end do
         call e4b(nfi(i),nfj,lbj,nii,li,lev)
         do j = 1,li
          jj = nii(j)
          call boxc(jj,2,nc)
          xjc = (nc(1)+.5)*rb
          yjc = (nc(2)+.5)*rb
          t = dcmplx(xic-xjc,yic-yjc)

          invt = 1/t
          sr(1,1) = invt
          invt2 = invt*invt
          nminus = 1
          do k = 1,mp
           c(k) = invt**k
          end do
          do n = 2,mp
           nminus = -nminus
           sr(n,1) = c(n)
           sr(1,n) = nminus*c(n)
           d(1) = c(n)
           mminus = nminus
           do m = 2,n-1
            mminus = -mminus
            d(m) = -(1+real(n-1)/(m-1))*d(m-1)*invt
            sr(n,m) = d(m)
            sr(m,n) = mminus*d(m)
           end do
           sr(n,n) = -2*invt2*(2-1.0/(n-1))*sr(n-1,n-1)
          end do

          do k = 1,mp
           do l = 1,mp
            aa(k) = aa(k)+sr(l,k)*b(l,nej(jj+1),lev)
           end do
          end do

         end do
         do j = 1,mp
          a(j,i) = aa(j)
         end do
        end do

        toc = tic
        call cpu_time(tic)
        tfmm(6) = tfmm(6)+tic-toc

       end do
       lev = lm
      end if

c Step 5. R-expansion and Final summation

      do i = 1,lbi
       call boxc(nfi(i),2,nc)
       xic = xmin+(nc(1)+.5)*rb
       yic = ymin+(nc(2)+.5)*rb
       call e2b(n2,n3,nbj,nfi(i),nfj,lbj,nfn,ln,lev)
       lj = 0
       do j = 1,ln
        jn = nej(nfn(j)+1)
        do k = ndj(jn,1),ndj(jn,2)
         lj = lj+1
         njj(lj) = k
        end do
       end do

       toc = tic
       call cpu_time(tic)
       tfmm(7) = tfmm(7)+tic-toc

       do j = ndi(i,1),ndi(i,2)
        xiic = xdi(j)-xic
        yiic = ydi(j)-yic
        ziic = dcmplx(xiic,yiic)

        t = 0
        do k = 1,mp
         t = t+a(k,i)*ziic**(k-1)
         uidu = -0.5/pi*imag(t)
         vidu = -0.5/pi*dble(t)
        end do

        toc = tic
        call cpu_time(tic)
        tfmm(8) = tfmm(8)+tic-toc

        do k = 1,lj
         kk = njj(k)
         xij = xdi(j)-xdj(kk)
         yij = ydi(j)-ydj(kk)
         rij = xij**2+yij**2
         if(rij.lt.eps)then
          rij = eps
         end if
         sij = sdj(kk)**2

         cutoff = 1/rij*(1-exp(-rij/2/sij))
c         cutoff = (rij**2+3*rij*sij+4*sij**2)/(rij+sij)**3
c         cutoff = (rij+2*sij)/(rij+sij)**2
c         cutoff = 1/(rij+sij)
         uidu = uidu+0.5/pi*gdj(kk)*yij*cutoff
         vidu = vidu-0.5/pi*gdj(kk)*xij*cutoff
        end do

        toc = tic
        call cpu_time(tic)
        tfmm(9) = tfmm(9)+tic-toc

        udi(j) = uidu
        vdi(j) = vidu
       end do
      end do
      n = 0
      do i = n0,n1
       n = n+1
       ui(ndbi(n)) = udi(n)
       vi(ndbi(n)) = vdi(n)
      end do
      print*,'init   : ',tfmm(1)
      print*,'index  : ',tfmm(7)
      print*,'p2p    : ',tfmm(9)
      print*,'p2m    : ',tfmm(2)
      print*,'m2m    : ',tfmm(3)
      print*,'m2l    : ',tfmm(4)+tfmm(6)
      print*,'l2l    : ',tfmm(5)
      print*,'l2p    : ',tfmm(8)
      print*,'------------------'
      print*,'total  : ',tfmm(1)+tfmm(2)+tfmm(3)+tfmm(4)+tfmm(5)+tfmm(6)
     *    +tfmm(7)+tfmm(8)+tfmm(9)
      return
      end

c Bookmarking box data structure

      subroutine boxdata(n0,n1,x,y,rd,lev,nb,ndb,nd,ne,nf,
     *    lb,rb,xmin,ymin)
      implicit real*8(a-h,o-z)

      include 'memory.f'

      dimension x(npmax),y(npmax)
      dimension nb(npmax),nbd(npmax),ndb(npmax),nd(npmax,2)
      dimension ne(nbmax),nf(nbmax),nxs(3,npmax)

      kl = 2**lev
      rb = rd/kl
      n = 0
      do i = n0,n1
       n = n+1
       nxs(1,n) = (x(i)-xmin)/rb
       nxs(2,n) = (y(i)-ymin)/rb
      end do

      call boxn(n,nxs,2,lev,nb)
      do i = 1,n
       nbd(i) = nb(i)
       ndb(i) = i+n0-1
      end do
      call sort(n,nbd,ndb)

      lb = 0
      nbo = -1
      do i = 1,n
       if(nbd(i).ne.nbo)then
        lb = lb+1
        ne(nbd(i)+1) = lb
        nf(lb) = nbd(i)
        nd(lb,1) = i
        if(lb.ge.2)then
         nd(lb-1,2) = i-1
        end if
        nbo = nbd(i)
       end if
      end do
      nd(lb,2) = n
      return
      end

c Retrive index of E3 boxes

      subroutine e3b(nf,nff,lb,nii,li,lev)
      implicit real*8(a-h,o-z)

      include 'memory.f'

      dimension nff(nbmax),nii(nbmax),nc(3),nfn(9)
      dimension nxs(3,npmax)

      kl = 2**lev
      call boxc(nf,2,nc)
      call neighbor(nc,kl,nxs,n)
      call boxn(n,nxs,2,lev,nfn)
      li = 0
      do i = 1,lb
       jj = 0
       do j = 1,n
        if(nff(i).eq.nfn(j))then
         jj = jj+1
        end if
       end do
       if(jj.eq.0)then
        li = li+1
        nii(li) = nff(i)
       end if
      end do
      return
      end

c Retrive index of E4 boxes

      subroutine e4b(nf,nff,lb,nii,li,lev)
      implicit real*8(a-h,o-z)

      include 'memory.f'

      dimension nff(nbmax),nii(nbmax),nc(3),nfn(9),nfpn(9)
      dimension nxs(3,npmax)

      kl = 2**lev
      nfp = nf/4
      call boxc(nfp,2,nc)
      call neighbor(nc,2*kl,nxs,np)
      call boxn(np,nxs,2,lev+1,nfpn)
      call boxc(nf,2,nc)
      call neighbor(nc,kl,nxs,n)
      call boxn(n,nxs,2,lev,nfn)
      li = 0
      do i = 1,np
       do j = 1,4
        nfnpn = nfpn(i)*4+j-1
        ii = 0
        do k = 1,n
         if(nfn(k).eq.nfnpn)then
          ii = ii+1
         end if
        end do
        if(ii.eq.0)then
         do k = 1,lb
          if(nff(k).eq.nfnpn)then
           li = li+1
           nii(li) = nfnpn
          end if
         end do
        end if
       end do
      end do
      return
      end

c Retrive index of E2 boxes

      subroutine e2b(n0,n1,nb,nf,nff,lb,nfn,n,lev)
      implicit real*8(a-h,o-z)

      include 'memory.f'

      dimension nb(npmax),nff(nbmax),nc(3),nfn(9),nfnn(9)
      dimension nxs(3,npmax)

      kl = 2**lev
      call boxc(nf,2,nc)
      call neighbor(nc,kl,nxs,nn)
      call boxn(nn,nxs,2,lev,nfnn)
      n = 0
      do i = 1,nn
       do j = 1,lb
        if(nff(j).eq.nfnn(i))then
         n = n+1
         nfn(n) = nfnn(i)
        end if
       end do
      end do
      return
      end

c Retrive box index from coordinate numbers

      subroutine boxn(n,nxs,ndi,lev,nb)
      implicit real*8(a-h,o-z)

      include 'memory.f'

      dimension nb(npmax),nxs(3,npmax)

      do i = 1,n
       nb(i) = 0
      end do
      do i = 1,lev
       do j = 1,n
        nbx = mod(nxs(1,j),2)
        nxs(1,j) = nxs(1,j)/2
        nb(j) = nb(j)+nbx*2**(ndi*(i-1)+1)
       end do

       if(ndi.ge.2)then
        do j = 1,n
         nby = mod(nxs(2,j),2)
         nxs(2,j) = nxs(2,j)/2
         nb(j) = nb(j)+nby*2**(ndi*(i-1))
        end do
       end if

       if(ndi.eq.3)then
        do j = 1,n
         nbz = mod(nxs(3,j),2)
         nxs(3,j) = nxs(3,j)/2
         nb(j) = nb(j)+nbz*2**(ndi*(i-1)+2)
        end do
       end if
      end do
      return
      end

c Retrive coordinate numbers from box index

      subroutine boxc(nf,ndi,nc)
      implicit real*8(a-h,o-z)

      dimension nc(3)

      do i = 1,3
       nc(i) = 0
      end do
      nb = nf
      k = 0
      i = 0
      do while(nb.ne.0)
       j = ndi-k
       nc(j) = nc(j)+mod(nb,2)*2**i
       nb = nb/2
       k = mod(k+1,ndi)
       if(k.eq.0)then
        i = i+1
       end if
      end do

      if(ndi.eq.3)then
       du = nc(1)
       nc(1) = nc(2)
       nc(2) = nc(3)
       nc(3) = du
      end if
      return
      end

c Retrive coordinate numbers of neighbor boxes

      subroutine neighbor(nc,kl,nxs,n)
      implicit real*8(a-h,o-z)

      include 'memory.f'

      dimension nc(3),nxs(3,npmax)

      if(nc(1).eq.0)then
       nx = 2
       nxi = 0
      elseif(nc(1).eq.kl-1)then
       nx = 2
       nxi = kl-2
      else
       nx = 3
       nxi = nc(1)-1
      end if

      if(nc(2).eq.0)then
       ny = 2
       nyi = 0
      elseif(nc(2).eq.kl-1)then
       ny = 2
       nyi = kl-2
      else
       ny = 3
       nyi = nc(2)-1
      end if

      n = 0
      do i = 1,nx
       do j = 1,ny
        n = n+1
        nxs(1,n) = nxi+i-1
        nxs(2,n) = nyi+j-1
       end do
      end do
      return
      end

c Quicksort

      subroutine sort(n,na,nb)
      implicit real*8(a-h,o-z)

      include 'memory.f'

      dimension na(npmax),nb(npmax),istack(50)

      jstack = 0
      m = 7
      l = 1
      ir = n
    1 if(ir-l.lt.m)then
       do j = l+1,ir
        naa = na(j)
        nbb = nb(j)
        do i = j-1,1,-1
         if(na(i).le.naa) goto 2
         na(i+1) = na(i)
         nb(i+1) = nb(i)
        end do
        i = 0
    2   na(i+1) = naa
        nb(i+1) = nbb
       end do
       if(jstack.eq.0) return
       ir = istack(jstack)
       l = istack(jstack-1)
       jstack = jstack-2
      else
       k = (l+ir)/2
       nd = na(k)
       na(k) = na(l+1)
       na(l+1) = nd
       nd = nb(k)
       nb(k) = nb(l+1)
       nb(l+1) = nd

       if(na(l+1).gt.na(ir))then
        nd = na(l+1)
        na(l+1) = na(ir)
        na(ir) = nd
        nd = nb(l+1)
        nb(l+1) = nb(ir)
        nb(ir) = nd
       end if

       if(na(l).gt.na(ir))then
        nd = na(l)
        na(l) = na(ir)
        na(ir) = nd
        nd = nb(l)
        nb(l) = nb(ir)
        nb(ir) = nd
       end if

       if(na(l+1).gt.na(l))then
        nd = na(l+1)
        na(l+1) = na(l)
        na(l) = nd
        nd = nb(l+1)
        nb(l+1) = nb(l)
        nb(l) = nd
       end if

       i = l+1
       j = ir
       naa = na(l)
       nbb = nb(l)
    3  continue
       i = i+1
       if(na(i).lt.naa) goto 3
    4  continue
       j = j-1
       if(na(j).gt.naa) goto 4
       if(j.lt.i) goto 5
       nd = na(i)
       na(i) = na(j)
       na(j) = nd
       nd = nb(i)
       nb(i) = nb(j)
       nb(j) = nd
       goto 3
    5  na(l) = na(j)
       na(j) = naa
       nb(l) = nb(j)
       nb(j) = nbb
       jstack = jstack+2
       if(jstack.gt.50) pause 'NSTACK too small in sort'
       if(ir-i+1.ge.j-1)then
        istack(jstack) = ir
        istack(jstack-1) = i
        ir = j-1
       else
        istack(jstack) = j-1
        istack(jstack-1) = l
        l = i
       end if
      end if
      goto 1
      return
      end

max_stars_repo_name: gforsyth/pyrbf

max_issues_repo_name: gforsyth/pyrbf

max_forks_repo_name: gforsyth/pyrbf


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine shape8hr(xl,xsj,shp,gs,a)
!
!     shape functions and derivatives for a 8-node linear isoparametric
!     mean strain solid element with hourglass control
!
!     Reference: Flanagan, D.P., Belytschko, T.; "Uniform  strain hexahedron
!     and quadrilateral with orthogonal Hourglass control". Int. J. Num.
!     Meth. Engg., Vol. 17, 679-706, 1981.
!
!     author: Otto-Ernst Bernhardi
!
      implicit none
!
      integer i,j,k,l
      real*8 shp(4,20),xl(3,20),xsj,vol
      real*8 x1,x2,x3,x4,x5,x6,x7,x8
      real*8 y1,y2,y3,y4,y5,y6,y7,y8
      real*8 z1,z2,z3,z4,z5,z6,z7,z8
      real*8 gb(8,4),gs(8,4),s0,a
!
c      data gb / 1.0d0, 1.0d0,-1.0d0,-1.0d0,-1.0d0,-1.0d0, 1.0d0, 1.0d0,
c     &          1.0d0,-1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0, 1.0d0,-1.0d0,
c     &          1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0,
c     &         -1.0d0, 1.0d0,-1.0d0, 1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0 /
      gb = reshape((
     &      / 1.0d0, 1.0d0,-1.0d0,-1.0d0,-1.0d0,-1.0d0, 1.0d0, 1.0d0,
     &        1.0d0,-1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0, 1.0d0,-1.0d0,
     &        1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0,
     &       -1.0d0, 1.0d0,-1.0d0, 1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0 /),
     &        (/8,4/))
!
!     shape functions and their global derivatives
!
!     shape functions
!
      shp(4, 1)=1.0d0/8.0d0
      shp(4, 2)=1.0d0/8.0d0
      shp(4, 3)=1.0d0/8.0d0
      shp(4, 4)=1.0d0/8.0d0
      shp(4, 5)=1.0d0/8.0d0
      shp(4, 6)=1.0d0/8.0d0
      shp(4, 7)=1.0d0/8.0d0
      shp(4, 8)=1.0d0/8.0d0
!
!     extract node point coordinates of element
      x1=xl(1,1)
      x2=xl(1,2)
      x3=xl(1,3)
      x4=xl(1,4)
      x5=xl(1,5)
      x6=xl(1,6)
      x7=xl(1,7)
      x8=xl(1,8)
      y1=xl(2,1)
      y2=xl(2,2)
      y3=xl(2,3)
      y4=xl(2,4)
      y5=xl(2,5)
      y6=xl(2,6)
      y7=xl(2,7)
      y8=xl(2,8)
      z1=xl(3,1)
      z2=xl(3,2)
      z3=xl(3,3)
      z4=xl(3,4)
      z5=xl(3,5)
      z6=xl(3,6)
      z7=xl(3,7)
      z8=xl(3,8)
!
!     Average displacement derivative operator matrix,
!     using eqn 16 in the reference above.
!     Generated using maxima/wxmaxima and the following input lines.
!     Note that shp array must be divided by the element volume.
!     h1:(1-r)*(1-s)*(1-t)/8;
!     h2:(1+r)*(1-s)*(1-t)/8;
!     h3:(1+r)*(1+s)*(1-t)/8;
!     h4:(1-r)*(1+s)*(1-t)/8;
!     h5:(1-r)*(1-s)*(1+t)/8;
!     h6:(1+r)*(1-s)*(1+t)/8;
!     h7:(1+r)*(1+s)*(1+t)/8;
!     h8:(1-r)*(1+s)*(1+t)/8;
!     H:matrix([h1,h2,h3,h4,h5,h6,h7,h8])$
!     Bx:diff(H,r);
!     By:diff(H,s);
!     Bz:diff(H,t);
!     B:matrix([Bx[1,1],Bx[1,2],Bx[1,3],Bx[1,4],Bx[1,5],Bx[1,6],Bx[1,7],Bx[1,8]],
!              [By[1,1],By[1,2],By[1,3],By[1,4],By[1,5],By[1,6],By[1,7],By[1,8]],
!              [Bz[1,1],Bz[1,2],Bz[1,3],Bz[1,4],Bz[1,5],Bz[1,6],Bz[1,7],Bz[1,8]]);
!     x:matrix([x1,x2,x3,x4,x5,x6,x7,x8],
!              [y1,y2,y3,y4,y5,y6,y7,y8],
!              [z1,z2,z3,z4,z5,z6,z7,z8]);
!     xr:B.transpose(x);
!     det:determinant(xr);
!     vol:factor(integrate(integrate(integrate(det,r,-1,1),s,-1,1),t,-1,1));
!     shp: matrix( [0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0])$
!     shp[1][1]:ratsimp(diff(vol,x1));
!     shp[1][2]:ratsimp(diff(vol,x2));
!     shp[1][3]:ratsimp(diff(vol,x3));
!     ... (more lines)
!     shp[3][8]:ratsimp(diff(vol,z8));
!     fortran(shp);
!
!     local derivatives of the shape functions: xi-derivative
!
      shp(1,1) = ((y5-y4)*z8+(y2-y5)*z6+(-y8+y6-y4+y2)*z5+(y8+y5-y3-y2)
     1   *z4+(y4-y2)*z3+(-y6-y5+y4+y3)*z2)/1.2d+1
      shp(1,2) = -((y6-y3)*z7+(-y7+y5-y3+y1)*z6+(y1-y6)*z5+(y3-y1)*z4+(
     1   y7+y6-y4-y1)*z3+(-y6-y5+y4+y3)*z1)/1.2d+1
      shp(1,3) = -((y7-y4)*z8+(-y8+y6-y4+y2)*z7+(y2-y7)*z6+(y8+y7-y2-y1
     1   )*z4+(-y7-y6+y4+y1)*z2+(y4-y2)*z1)/1.2d+1
      shp(1,4) = -((y7-y5+y3-y1)*z8+(y3-y8)*z7+(y8-y1)*z5+(-y8-y7+y2+y1
     1   )*z3+(y1-y3)*z2+(y8+y5-y3-y2)*z1)/1.2d+1
      shp(1,5) = ((y7+y6-y4-y1)*z8+(y6-y8)*z7+(-y8-y7+y2+y1)*z6+(y8-y1)
     1   *z4+(y1-y6)*z2+(y8-y6+y4-y2)*z1)/1.2d+1
      shp(1,6) = ((y7-y5)*z8+(-y8-y5+y3+y2)*z7+(y8+y7-y2-y1)*z5+(y2-y7)
     1   *z3+(-y7+y5-y3+y1)*z2+(y5-y2)*z1)/1.2d+1
      shp(1,7) = -((y6+y5-y4-y3)*z8+(-y8-y5+y3+y2)*z6+(y6-y8)*z5+(y8-y3
     1   )*z4+(y8-y6+y4-y2)*z3+(y3-y6)*z2)/1.2d+1
      shp(1,8) = ((y6+y5-y4-y3)*z7+(y5-y7)*z6+(-y7-y6+y4+y1)*z5+(y7-y5+
     1   y3-y1)*z4+(y7-y4)*z3+(y4-y5)*z1)/1.2d+1
!
!     local derivatives of the shape functions: eta-derivative
!
      shp(2,1) = -((x5-x4)*z8+(x2-x5)*z6+(-x8+x6-x4+x2)*z5+(x8+x5-x3-x2
     1   )*z4+(x4-x2)*z3+(-x6-x5+x4+x3)*z2)/1.2d+1
      shp(2,2) = ((x6-x3)*z7+(-x7+x5-x3+x1)*z6+(x1-x6)*z5+(x3-x1)*z4+(x
     1   7+x6-x4-x1)*z3+(-x6-x5+x4+x3)*z1)/1.2d+1
      shp(2,3) = ((x7-x4)*z8+(-x8+x6-x4+x2)*z7+(x2-x7)*z6+(x8+x7-x2-x1)
     1   *z4+(-x7-x6+x4+x1)*z2+(x4-x2)*z1)/1.2d+1
      shp(2,4) = ((x7-x5+x3-x1)*z8+(x3-x8)*z7+(x8-x1)*z5+(-x8-x7+x2+x1)
     1   *z3+(x1-x3)*z2+(x8+x5-x3-x2)*z1)/1.2d+1
      shp(2,5) = -((x7+x6-x4-x1)*z8+(x6-x8)*z7+(-x8-x7+x2+x1)*z6+(x8-x1
     1   )*z4+(x1-x6)*z2+(x8-x6+x4-x2)*z1)/1.2d+1
      shp(2,6) = -((x7-x5)*z8+(-x8-x5+x3+x2)*z7+(x8+x7-x2-x1)*z5+(x2-x7
     1   )*z3+(-x7+x5-x3+x1)*z2+(x5-x2)*z1)/1.2d+1
      shp(2,7) = ((x6+x5-x4-x3)*z8+(-x8-x5+x3+x2)*z6+(x6-x8)*z5+(x8-x3)
     1   *z4+(x8-x6+x4-x2)*z3+(x3-x6)*z2)/1.2d+1
      shp(2,8) = -((x6+x5-x4-x3)*z7+(x5-x7)*z6+(-x7-x6+x4+x1)*z5+(x7-x5
     1   +x3-x1)*z4+(x7-x4)*z3+(x4-x5)*z1)/1.2d+1
!
!     local derivatives of the shape functions: zeta-derivative
!
      shp(3,1) = ((x5-x4)*y8+(x2-x5)*y6+(-x8+x6-x4+x2)*y5+(x8+x5-x3-x2)
     1   *y4+(x4-x2)*y3+(-x6-x5+x4+x3)*y2)/1.2d+1
      shp(3,2) = -((x6-x3)*y7+(-x7+x5-x3+x1)*y6+(x1-x6)*y5+(x3-x1)*y4+(
     1   x7+x6-x4-x1)*y3+(-x6-x5+x4+x3)*y1)/1.2d+1
      shp(3,3) = -((x7-x4)*y8+(-x8+x6-x4+x2)*y7+(x2-x7)*y6+(x8+x7-x2-x1
     1   )*y4+(-x7-x6+x4+x1)*y2+(x4-x2)*y1)/1.2d+1
      shp(3,4) = -((x7-x5+x3-x1)*y8+(x3-x8)*y7+(x8-x1)*y5+(-x8-x7+x2+x1
     1   )*y3+(x1-x3)*y2+(x8+x5-x3-x2)*y1)/1.2d+1
      shp(3,5) = ((x7+x6-x4-x1)*y8+(x6-x8)*y7+(-x8-x7+x2+x1)*y6+(x8-x1)
     1   *y4+(x1-x6)*y2+(x8-x6+x4-x2)*y1)/1.2d+1
      shp(3,6) = ((x7-x5)*y8+(-x8-x5+x3+x2)*y7+(x8+x7-x2-x1)*y5+(x2-x7)
     1   *y3+(-x7+x5-x3+x1)*y2+(x5-x2)*y1)/1.2d+1
      shp(3,7) = -((x6+x5-x4-x3)*y8+(-x8-x5+x3+x2)*y6+(x6-x8)*y5+(x8-x3
     1   )*y4+(x8-x6+x4-x2)*y3+(x3-x6)*y2)/1.2d+1
      shp(3,8) = ((x6+x5-x4-x3)*y7+(x5-x7)*y6+(-x7-x6+x4+x1)*y5+(x7-x5+
     1   x3-x1)*y4+(x7-x4)*y3+(x4-x5)*y1)/1.2d+1
!
!     computation of element volume (eqn 15)
!
      vol=0.0d0
      do k=1,8
        vol=vol+xl(1,k)*shp(1,k)
      enddo
!
!     computation of the average jacobian determinant
!
      xsj=vol/8.0d0
!
!     hourglass control vectors(see appendix 2 from the reference above).
!     divide shp array by element volume
      a=0.0d0
      do i=1,8
         do j=1,3
            a=a+shp(j,i)*shp(j,i)
            s0=shp(j,i)/vol
            shp(j,i)=s0
            do k=1,4
               gs(i,k)=gb(i,k)
               do l=1,8
                  gs(i,k)=gs(i,k)-s0*xl(j,l)*gb(l,k)
               enddo
            enddo
         enddo
      enddo
c
c     calculate hourglass control stiffness factor a
c     (to be used in hgstiffness() and hgforce())
c
c     in ABAQUS, a 0.005 is used as default value.
      a=0.0005*a/vol
c
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
      program tester

          use DiscreteSampler

          type(binID) myBinID
          real*8 mValue
          integer i
          real*8 jacobian
          integer bid

          call DS_register_dimension('Dimension1',20)
          mValue=1.0d1
          do i=1,20
            mValue=mValue+1.0d1
            call DS_add_entry('Dimension1',i,mValue)
          enddo
          do i=1,20
            mValue=mValue+1.0d1
            call DS_add_entry('Dimension1',i,mValue)
            if (i.eq.10) then
            call DS_add_entry('Dimension1',i,1.0d0)
            endif
          enddo
          call DS_update_grid()
          call DS_write_grid('grids.dsg')
          call DS_clear()
          call DS_register_dimension('Dimension2',30)
          do i=1,30
            mValue=mValue+1.0d1
            call DS_add_entry('Dimension2',i,mValue)
          enddo
          do i=1,30
            mValue=mValue+1.0d1
            call DS_add_entry('Dimension2',i,mValue)
            if (i.eq.10) then
            call DS_add_entry('Dimension2',i,1.0d0)
            endif
          enddo
          call DS_update_grid()
          call DS_write_grid('grids.dsg')
          call DS_clear()
          write(*,*) 'Before grid load'
          call DS_print_global_info()
          call DS_load_grid('grids.dsg')
          write(*,*) 'after grid load'
          call DS_print_global_info()

          call DS_get_point('Dimension1',
     &           0.4d0, bid, jacobian, 'norm')
          write(*,*) 'I got point ',bid,' with jac',jacobian

          call DS_get_point('Dimension1',
     &           0.1d0, bid, jacobian, 'norm')
          write(*,*) 'I got point ',bid,' with jac',jacobian

          call DS_get_point('Dimension2',
     &           0.7d0, bid, jacobian,'variance')
          write(*,*) 'I got point ',bid,' with jac',jacobian

      end program tester


text:
      subroutine setprob
      implicit double precision (a-h,o-z)
      common /comrp/ grav
      common /comsrc/ ndim
      common/cdisc/ x0,y0,alf,beta,r0,idisc
      common /comic/ hin,hout
c
c     # Set the material parameters for the acoustic equations
c     # Passed to the Riemann solver rp1.f in a common block
c
      open(unit=7,file='setprob.data',status='old',form='formatted')

c     # ndim = space dimensions (2 = cylindrical symmetry, 3 = spherical)
      read(7,*) ndim

c     # gravitational constant:
      read(7,*) grav

c     # data for radial dam-break problem:
      read(7,*) x0,y0,r0
      read(7,*) hin,hout
c
      return
      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw