text:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C Program BSL.FOR
C
C Purpose: Interactive data appraisal for 2-d images.
C
C Author: G.R.Mant & J.Bordas, Daresbury lab., Warrington
C
C Limitations: Maximum image size 512*512 ( MAXDIM in common /SIZES/ )
C
C Non-ANSI code:
C 1: All routines called by BSL main program, including itself, use
C    the INCLUDE statement for common block insertion.
C 2: Wherever information is requested from the terminal the '$' format
C    control character is used. This leaves the cursor at the end of the
C    line.
C
C Implementation notes:
C 1: Subroutine ERRMSG uses ESCAPE sequences specific to VT100
C    emulator terminals to switch into inverse video mode.
C 2: Several subroutines are taken from the OTOKO Program library.
C 3: Subroutines requiring changes for unix implementation:
C
C Updates:
C 09/06/89 GM Initail implementation of UNIX version of BSL started
C 30/05/96 GM Corrected bug in .VER & .HOR, changed askno to askyes
C             Corrected bug in .SUM (for 2nd picture crash).
C
      IMPLICIT NONE
      INCLUDE 'COMMON.FOR'
      COMMON /MYGRAF/ OPENGR
      LOGICAL         OPENGR
C
C CALLS 43:  ADD    , ADDCON  , ADDFIL , ADDNRM  , AVERAG , BACK
C            CHANGE , CONPLOT , CUT    , DISPLAY , DIVCON , DIVFIL
C            DIVIDE , DIVNRM  , DUPE   , EXPON   , FFT    , GETINS
C            HORINT , IFT     , INTEG  , INTERP  , INVPOL , LOGTHM
C            MASK   , MAXVAL  , MIRROR , MULCON  , MULFIL , MULNRM
C            MULT   , PACK    , POLAR  , POWER   , PRTVAL , RADII
C            REMAP  , REPLACE , ROTATE , SHIFT   , SUBNRM , SUM
C            VERINT , WINDO   , ZERO   , SECTOR  , IMPRINT, CONVOLVE
C            MERGE  , GAUSSIAN, INSERT , ROTINT  , CIRSCN , RADSCN
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C LOCAL VARIABLES:
C
      INTEGER IOPT
C
C IOPT   : INSTRUCTION SELECTED
C
C-----------------------------------------------------------------------
C
C========GET INSTRUCTION FROM USER OR <CTRL-Z> TO END
C
10    CALL GETINS (ITERM,IPRINT,IOPT)
      IF (IOPT.EQ.0) GOTO 999
C
C========PROCESS INSTRUCTION
C
      GOTO (100,110,120,130,140,150,160,170,180,190,
     &      200,210,220,230,240,250,260,270,280,290,
     &      300,310,320,330,340,350,360,370,380,390,
     &      400,410,420,430,440,450,460,470,480,490,
     &      500,510,520,530,540,550,190,570,580,590,
     &      600,610,620,630,640,650,660,670,680) IOPT
C
      CALL ERRMSG ('Error: In software Please report!')
      GOTO 999
C
100   CALL ADDCON
      GOTO 10
C
110   CALL ADD
      GOTO 10
C
120   CALL ADDNRM
      GOTO 10
C
130   CALL ADDFIL
      GOTO 10
C
140   CALL AVERAG
      GOTO 10
C
150   CALL CONPLOT
      GOTO 10
C
160   CALL CUT
      GOTO 10
C
170   CALL DIVCON
      GOTO 10
C
180   CALL DIVNRM
      GOTO 10
C
190   CALL DISPLAY
      GOTO 10
C
200   CALL DIVIDE
      GOTO 10
C
210   CALL DIVFIL
      GOTO 10
C
220   CALL DUPE
      GOTO 10
C
230   CALL EXPON
      GOTO 10
C
240   CALL FFT
      GOTO 10
C
250   CALL HORINT
      GOTO 10
C
260   CALL IFT
      GOTO 10
C
270   CALL INTEG
      GOTO 10
C
280   CALL INVPOL
      GOTO 10
C
290   CALL INTERP
      GOTO 10
C
300   CALL LOGTHM
      GOTO 10
C
310   CALL MAXVAL
      GOTO 10
C
320   CALL MIRROR
      GOTO 10
C
330   CALL MULFIL
      GOTO 10
C
340   CALL MASK
      GOTO 10
C
350   CALL MULCON
      GOTO 10
C
360   CALL MULT
      GOTO 10
C
370   CALL MULNRM
      GOTO 10
C
380   CALL PACK
      GOTO 10
C
390   CALL POLAR
      GOTO 10
C
400   CALL POWER
      GOTO 10
C
410   CALL PRTVAL
      GOTO 10
C
420   CALL RADII
      GOTO 10
C
430   CALL DUPE
      GOTO 10
C
440   CALL REPLACE
      GOTO 10
C
450   CALL REMAP
      GOTO 10
C
460   CALL ROTATE
      GOTO 10
C
470   CALL SHIFT
      GOTO 10
C
480   CALL SUM
      GOTO 10
C
490   CALL SUBNRM
      GOTO 10
C
500   CALL VERINT
      GOTO 10
C
510   CALL WINDO
      GOTO 10
C
520   CALL ZERO
      GOTO 10
C
530   CALL BACK
      GOTO 10
C
540   CALL CHANGE
      GOTO 10
C
550   CALL GLITCH
      GOTO 10
C
 570  CALL IMPRINT
      GOTO 10
C
 580  CALL SECTOR
      GOTO 10
C
 590  CALL GAUSSIAN
      GOTO 10
C
 600  CALL CONVOLVE
      GOTO 10
C
 610  CALL MERGE
      GOTO 10
C
 620  CALL INSERT
      GOTO 10
C
 630  CALL ROTINT
      GOTO 10
C
 670  CONTINUE
 640  CALL CIRSCN
      GOTO 10
C
 680  CONTINUE
 650  CALL RADSCN
      GOTO 10
C
 660  CALL SURFACE
      GOTO 10
C
999   IF (OPENGR) CALL GREND
      STOP
      END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      BLOCK DATA BLKDAT
C
      INCLUDE 'COMMON.FOR'
      COMMON /MYGRAF/ OPENGR
      LOGICAL         OPENGR
C
C-----------------------------------------------------------------------
C
C========INITIALISE FORTRAN I/O CHANNEL UNITS
C
      DATA JUNIT/-1/, KUNIT/-1/, LUNIT/-1/
      DATA ITERM/5/, IPRINT/6/, IUNIT/1/
C
      DATA MAXDIM/1048576/
C
      DATA OPENGR/.FALSE./
      END

ext: f


text:
MODULE test_mod
    IMPLICIT NONE

    DOUBLE PRECISION, POINTER :: x_comp_sp(:)
    DOUBLE PRECISION, POINTER :: concen(:,:)

    INTEGER :: nx = 5

END MODULE test_mod

SUBROUTINE use_array(value1,value2)
    USE test_mod
    IMPLICIT NONE

    DOUBLE PRECISION value1, value2
    INTEGER i,j

    ALLOCATE(x_comp_sp(nx), concen(nx,nx))

    do i = 1,nx
        x_comp_sp(i) = value1
        do j = 1,nx
            concen(i,j) = value2
        end do
    end do

end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
subroutine setaux1d(ixyz,mbc,mx1,maxm,x1lower,x2val,x3val,dx1,dx2,dx3,t,maux,aux,iout)

    ! set auxiliary variables in current row with direction specified by ixyz
    !
    ! ixyz = 1 => row in x direction, e2 is y direction, e3 is z direction
    ! ixyz = 2 => row in y direction, e2 is z direction, e3 is x direction
    ! ixyz = 3 => row in z direction, e2 is x direction, e3 is y direction
    !
    ! Auxiliary variables:
    !       1 rho
    !       2 lambda
    !       3 mu
    !       4 cp
    !       5 cs
    !       6 slip (obtained from 2d aux in step3.f)
    !       7 nx at lower wall in e1 direction
    !       8 ny at lower wall in e1 direction
    !       9 nz at lower wall in e1 direction
    !       10 area ratio of lower wall in e1 direction
    !       11 nx at lower wall in e2 direction
    !       12 ny at lower wall in e2 direction
    !       13 nz at lower wall in e2 direction
    !       14 area ratio of lower wall in e2 direction
    !       15 nx at lower wall in e3 direction
    !       16 ny at lower wall in e3 direction
    !       17 nz at lower wall in e3 direction
    !       18 area ratio of lower wall in e3 direction

    implicit none
    integer, intent(in) :: ixyz,mbc,mx1,maxm,maux,iout
    real(kind=8), intent(in) :: x1lower,x2val,x3val,dx1,dx2,dx3,t
    real(kind=8), intent(out) ::  aux(maux,1-mbc:maxm+mbc,3)

    real(kind=8) :: x1cell, cp, cs, rho_cell, lambda_cell, mu_cell
    real(kind=8) :: x1ccorn(4), x2ccorn(4), x3ccorn(4)
    real(kind=8) :: x1pcorn(4), x2pcorn(4), x3pcorn(4), mag
    real(kind=8) :: nx1, nx2, nx3
    integer :: i

    lambda_cell = 60.d9  ! Pa
    mu_cell = 30.d9      ! Pa
    rho_cell = 2500.d0   ! kg/m**3
    cp = dsqrt((lambda_cell + 2.d0*mu_cell)/rho_cell)
    cs = dsqrt(mu_cell/rho_cell)

    do i=1-mbc,mx1 + mbc
      x1cell = x1lower + (i-0.5d0)*dx1
      aux(1,i,iout) = rho_cell
      aux(2,i,iout) = lambda_cell
      aux(3,i,iout) = mu_cell
      aux(4,i,iout) = cp
      aux(5,i,iout) = cs

      ! compute mapping info for lower face in e1 direction
      x1ccorn(1) = x1cell - 0.5d0*dx1 ! p1 dot e1
      x2ccorn(1) = x2val - 0.5d0*dx2  ! p1 dot e2
      x3ccorn(1) = x3val - 0.5d0*dx3  ! p1 dot e3

      x1ccorn(2) = x1cell - 0.5d0*dx1 ! p2 dot e1
      x2ccorn(2) = x2val - 0.5d0*dx2  ! p2 dot e2
      x3ccorn(2) = x3val + 0.5d0*dx3  ! p2 dot e3

      x1ccorn(3) = x1cell - 0.5d0*dx1 ! p3 dot e1
      x2ccorn(3) = x2val + 0.5d0*dx2  ! p3 dot e2
      x3ccorn(3) = x3val + 0.5d0*dx3  ! p3 dot e3

      x1ccorn(4) = x1cell - 0.5d0*dx1 ! p4 dot e1
      x2ccorn(4) = x2val + 0.5d0*dx2  ! p4 dot e2
      x3ccorn(4) = x3val - 0.5d0*dx3  ! p4 dot e3
      if (ixyz == 1) then
        call mapc2p(x1ccorn(1), x2ccorn(1), x3ccorn(1), x1pcorn(1), x2pcorn(1), x3pcorn(1))
        call mapc2p(x1ccorn(2), x2ccorn(2), x3ccorn(2), x1pcorn(2), x2pcorn(2), x3pcorn(2))
        call mapc2p(x1ccorn(3), x2ccorn(3), x3ccorn(3), x1pcorn(3), x2pcorn(3), x3pcorn(3))
        call mapc2p(x1ccorn(4), x2ccorn(4), x3ccorn(4), x1pcorn(4), x2pcorn(4), x3pcorn(4))
      else if (ixyz == 2) then
        call mapc2p(x3ccorn(1), x1ccorn(1), x2ccorn(1), x3pcorn(1), x1pcorn(1), x2pcorn(1))
        call mapc2p(x3ccorn(2), x1ccorn(2), x2ccorn(2), x3pcorn(2), x1pcorn(2), x2pcorn(2))
        call mapc2p(x3ccorn(3), x1ccorn(3), x2ccorn(3), x3pcorn(3), x1pcorn(3), x2pcorn(3))
        call mapc2p(x3ccorn(4), x1ccorn(4), x2ccorn(4), x3pcorn(4), x1pcorn(4), x2pcorn(4))
      else if (ixyz == 3) then
        call mapc2p(x2ccorn(1), x3ccorn(1), x1ccorn(1), x2pcorn(1), x3pcorn(1), x1pcorn(1))
        call mapc2p(x2ccorn(2), x3ccorn(2), x1ccorn(2), x2pcorn(2), x3pcorn(2), x1pcorn(2))
        call mapc2p(x2ccorn(3), x3ccorn(3), x1ccorn(3), x2pcorn(3), x3pcorn(3), x1pcorn(3))
        call mapc2p(x2ccorn(4), x3ccorn(4), x1ccorn(4), x2pcorn(4), x3pcorn(4), x1pcorn(4))
      end if

      nx1 = (x2pcorn(3) - x2pcorn(1))*(x3pcorn(2) - x3pcorn(4)) - (x2pcorn(2) - x2pcorn(4))*(x3pcorn(3) - x3pcorn(1))
      nx2 = (x1pcorn(2) - x1pcorn(4))*(x3pcorn(3) - x3pcorn(1)) - (x1pcorn(3) - x1pcorn(1))*(x3pcorn(2) - x3pcorn(4))
      nx3 = (x1pcorn(3) - x1pcorn(1))*(x2pcorn(2) - x2pcorn(4)) - (x1pcorn(2) - x1pcorn(4))*(x2pcorn(3) - x2pcorn(1))
      mag = dsqrt(nx1*nx1 + nx2*nx2 + nx3*nx3)
      nx1 = nx1/mag
      nx2 = nx2/mag
      nx3 = nx3/mag

      if (ixyz == 1) then
        aux(7,i,iout) = nx1
        aux(8,i,iout) = nx2
        aux(9,i,iout) = nx3
      else if (ixyz == 2) then
        aux(7,i,iout) = nx3
        aux(8,i,iout) = nx1
        aux(9,i,iout) = nx2
      else if (ixyz == 3) then
        aux(7,i,iout) = nx2
        aux(8,i,iout) = nx3
        aux(9,i,iout) = nx1
      end if

      aux(10,i,iout) = 0.5d0*mag/(dx2*dx3)

      ! compute mapping info for lower face in e2 direction
      x1ccorn(1) = x1cell - 0.5d0*dx1 ! p1 dot e1
      x2ccorn(1) = x2val - 0.5d0*dx2  ! p1 dot e2
      x3ccorn(1) = x3val - 0.5d0*dx3  ! p1 dot e3

      x1ccorn(2) = x1cell + 0.5d0*dx1 ! p2 dot e1
      x2ccorn(2) = x2val - 0.5d0*dx2  ! p2 dot e2
      x3ccorn(2) = x3val - 0.5d0*dx3  ! p2 dot e3

      x1ccorn(3) = x1cell + 0.5d0*dx1 ! p3 dot e1
      x2ccorn(3) = x2val - 0.5d0*dx2  ! p3 dot e2
      x3ccorn(3) = x3val + 0.5d0*dx3  ! p3 dot e3

      x1ccorn(4) = x1cell - 0.5d0*dx1 ! p4 dot e1
      x2ccorn(4) = x2val - 0.5d0*dx2  ! p4 dot e2
      x3ccorn(4) = x3val + 0.5d0*dx3  ! p4 dot e3
      if (ixyz == 1) then
        call mapc2p(x1ccorn(1), x2ccorn(1), x3ccorn(1), x1pcorn(1), x2pcorn(1), x3pcorn(1))
        call mapc2p(x1ccorn(2), x2ccorn(2), x3ccorn(2), x1pcorn(2), x2pcorn(2), x3pcorn(2))
        call mapc2p(x1ccorn(3), x2ccorn(3), x3ccorn(3), x1pcorn(3), x2pcorn(3), x3pcorn(3))
        call mapc2p(x1ccorn(4), x2ccorn(4), x3ccorn(4), x1pcorn(4), x2pcorn(4), x3pcorn(4))
      else if (ixyz == 2) then
        call mapc2p(x3ccorn(1), x1ccorn(1), x2ccorn(1), x3pcorn(1), x1pcorn(1), x2pcorn(1))
        call mapc2p(x3ccorn(2), x1ccorn(2), x2ccorn(2), x3pcorn(2), x1pcorn(2), x2pcorn(2))
        call mapc2p(x3ccorn(3), x1ccorn(3), x2ccorn(3), x3pcorn(3), x1pcorn(3), x2pcorn(3))
        call mapc2p(x3ccorn(4), x1ccorn(4), x2ccorn(4), x3pcorn(4), x1pcorn(4), x2pcorn(4))
      else if (ixyz == 3) then
        call mapc2p(x2ccorn(1), x3ccorn(1), x1ccorn(1), x2pcorn(1), x3pcorn(1), x1pcorn(1))
        call mapc2p(x2ccorn(2), x3ccorn(2), x1ccorn(2), x2pcorn(2), x3pcorn(2), x1pcorn(2))
        call mapc2p(x2ccorn(3), x3ccorn(3), x1ccorn(3), x2pcorn(3), x3pcorn(3), x1pcorn(3))
        call mapc2p(x2ccorn(4), x3ccorn(4), x1ccorn(4), x2pcorn(4), x3pcorn(4), x1pcorn(4))
      end if

      nx1 = (x2pcorn(3) - x2pcorn(1))*(x3pcorn(2) - x3pcorn(4)) - (x2pcorn(2) - x2pcorn(4))*(x3pcorn(3) - x3pcorn(1))
      nx2 = (x1pcorn(2) - x1pcorn(4))*(x3pcorn(3) - x3pcorn(1)) - (x1pcorn(3) - x1pcorn(1))*(x3pcorn(2) - x3pcorn(4))
      nx3 = (x1pcorn(3) - x1pcorn(1))*(x2pcorn(2) - x2pcorn(4)) - (x1pcorn(2) - x1pcorn(4))*(x2pcorn(3) - x2pcorn(1))
      mag = dsqrt(nx1*nx1 + nx2*nx2 + nx3*nx3)
      nx1 = nx1/mag
      nx2 = nx2/mag
      nx3 = nx3/mag

      if (ixyz == 1) then
        aux(11,i,iout) = nx1
        aux(12,i,iout) = nx2
        aux(13,i,iout) = nx3
      else if (ixyz == 2) then
        aux(11,i,iout) = nx3
        aux(12,i,iout) = nx1
        aux(13,i,iout) = nx2
      else if (ixyz == 3) then
        aux(11,i,iout) = nx2
        aux(12,i,iout) = nx3
        aux(13,i,iout) = nx1
      end if

      aux(14,i,iout) = 0.5d0*mag/(dx1*dx3)

      ! compute mapping info for lower face in e3 direction
      x1ccorn(1) = x1cell - 0.5d0*dx1 ! p1 dot e1
      x2ccorn(1) = x2val - 0.5d0*dx2  ! p1 dot e2
      x3ccorn(1) = x3val - 0.5d0*dx3  ! p1 dot e3

      x1ccorn(2) = x1cell - 0.5d0*dx1 ! p2 dot e1
      x2ccorn(2) = x2val + 0.5d0*dx2  ! p2 dot e2
      x3ccorn(2) = x3val - 0.5d0*dx3  ! p2 dot e3

      x1ccorn(3) = x1cell + 0.5d0*dx1 ! p3 dot e1
      x2ccorn(3) = x2val + 0.5d0*dx2  ! p3 dot e2
      x3ccorn(3) = x3val - 0.5d0*dx3  ! p3 dot e3

      x1ccorn(4) = x1cell + 0.5d0*dx1 ! p4 dot e1
      x2ccorn(4) = x2val - 0.5d0*dx2  ! p4 dot e2
      x3ccorn(4) = x3val - 0.5d0*dx3  ! p4 dot e3
      if (ixyz == 1) then
        call mapc2p(x1ccorn(1), x2ccorn(1), x3ccorn(1), x1pcorn(1), x2pcorn(1), x3pcorn(1))
        call mapc2p(x1ccorn(2), x2ccorn(2), x3ccorn(2), x1pcorn(2), x2pcorn(2), x3pcorn(2))
        call mapc2p(x1ccorn(3), x2ccorn(3), x3ccorn(3), x1pcorn(3), x2pcorn(3), x3pcorn(3))
        call mapc2p(x1ccorn(4), x2ccorn(4), x3ccorn(4), x1pcorn(4), x2pcorn(4), x3pcorn(4))
      else if (ixyz == 2) then
        call mapc2p(x3ccorn(1), x1ccorn(1), x2ccorn(1), x3pcorn(1), x1pcorn(1), x2pcorn(1))
        call mapc2p(x3ccorn(2), x1ccorn(2), x2ccorn(2), x3pcorn(2), x1pcorn(2), x2pcorn(2))
        call mapc2p(x3ccorn(3), x1ccorn(3), x2ccorn(3), x3pcorn(3), x1pcorn(3), x2pcorn(3))
        call mapc2p(x3ccorn(4), x1ccorn(4), x2ccorn(4), x3pcorn(4), x1pcorn(4), x2pcorn(4))
      else if (ixyz == 3) then
        call mapc2p(x2ccorn(1), x3ccorn(1), x1ccorn(1), x2pcorn(1), x3pcorn(1), x1pcorn(1))
        call mapc2p(x2ccorn(2), x3ccorn(2), x1ccorn(2), x2pcorn(2), x3pcorn(2), x1pcorn(2))
        call mapc2p(x2ccorn(3), x3ccorn(3), x1ccorn(3), x2pcorn(3), x3pcorn(3), x1pcorn(3))
        call mapc2p(x2ccorn(4), x3ccorn(4), x1ccorn(4), x2pcorn(4), x3pcorn(4), x1pcorn(4))
      end if

      nx1 = (x2pcorn(3) - x2pcorn(1))*(x3pcorn(2) - x3pcorn(4)) - (x2pcorn(2) - x2pcorn(4))*(x3pcorn(3) - x3pcorn(1))
      nx2 = (x1pcorn(2) - x1pcorn(4))*(x3pcorn(3) - x3pcorn(1)) - (x1pcorn(3) - x1pcorn(1))*(x3pcorn(2) - x3pcorn(4))
      nx3 = (x1pcorn(3) - x1pcorn(1))*(x2pcorn(2) - x2pcorn(4)) - (x1pcorn(2) - x1pcorn(4))*(x2pcorn(3) - x2pcorn(1))
      mag = dsqrt(nx1*nx1 + nx2*nx2 + nx3*nx3)
      nx1 = nx1/mag
      nx2 = nx2/mag
      nx3 = nx3/mag

      if (ixyz == 1) then
        aux(15,i,iout) = nx1
        aux(16,i,iout) = nx2
        aux(17,i,iout) = nx3
      else if (ixyz == 2) then
        aux(15,i,iout) = nx3
        aux(16,i,iout) = nx1
        aux(17,i,iout) = nx2
      else if (ixyz == 3) then
        aux(15,i,iout) = nx2
        aux(16,i,iout) = nx3
        aux(17,i,iout) = nx1
      end if

      aux(18,i,iout) = 0.5d0*mag/(dx1*dx2)

    end do

end subroutine setaux1d

max_stars_repo_name: rjleveque/seismic

max_issues_repo_name: clawpack/seismic

max_forks_repo_name: rjleveque/seismic


text:
    !-------------------------------------------------------------------------------------------------------------
    !
    !> \file        ModuleGEMSolver.f90
    !> \brief       Fortran module for input/output of the non-linear solver.
    !> \author      M.H.A. Piro
    !
    !
    ! Pertinent variables:
    ! ====================
    !
    !
    !> \param iterLast          The last global iteration that the phase assemblage was adjusted.
    !> \param iterHistory       An integer matrix representing all of the indices of phases that contribute
    !!                           to the equilibrium phase assemblage at each stage in the iteration history.
    !> \param dGEMFunctionNorm  A double real scalar representing the norm of the functional vector in the
    !!                          GEMSolver.
    !> \param dGEMFunctionNormLast A double real scalar representing the norm of the functional vector in the
    !!                              GEMSolver from the last iteration.
    !> \param dSumMolFractionSoln  A double real vector representing the sum of mole fractions in each solution
    !!                              phase.
    !> \param dUpdateVar        A double real vector representing the direction vector that updates the function
    !                            vector.
    !> \param dPartialExcessGibbs  A double real vector representing the partial molar excess Gibbs energy of
    !!                              mixing of each species in the system.
    !> \param dEffStoichSolnPhase   A double real matrix representing the effective stoichiometry of each
    !!                               solution phase.
    !> \param lDebugMode        A logical variable used for debugging purposes.  When it is TRUE, a number
    !!                           of print statements are applied.
    !> \param lRevertSystem     A logical variable identifying whether the system should be reverted (TRUE)
    !!                           or not (FALSE).
    !> \param lConverged        A logical variable identifying whether the system has converged (TRUE) or not
    !!                           (FALSE).
    !> \param lSolnPhases       A logical vector indicating whether a particular solution phase is currently
    !!                           assumed to be stable (true) or not (false).
    !> \param lMiscibility      A logical vector indicating whether a particular solution phase has a
    !!                           miscibility gap (true) or not (false).
    !
    !
    ! CONSIDER REMOVING THE FOLLOWING VARIABLES:
    !
    ! iterLastCon           The last iteration that a pure condensed phase was either added to or removed from
    !                        the estimated phase assemblage.
    ! iterLastSoln          The last iteration that a pure condensed phase was either added to or removed from
    !                        the estimated phase assemblage.
    ! iConPhaseLast         The species index of the last pure condensed phase that was either added to, or
    !                        removed from, the estimated phase assemblage.
    ! iSolnPhaseLast        The species index of the last pure condensed phase that was either added to, or
    !                        removed from, the estimated phase assemblage.
    !
    !
    !-------------------------------------------------------------------------------------------------------------

module ModuleGEMSolver

    implicit none

    SAVE

    integer                              ::  iterLast,      iterStep, iterRevert, iterGlobal
    integer                              ::  iterLastCon,   iterLastSoln,         iterSwap,   iterLastMiscGapCheck
    integer                              ::  iConPhaseLast, iSolnPhaseLast,       iSolnSwap,  iPureConSwap
    integer,                 parameter   ::  iterGlobalMax = 3000
    integer, dimension(:,:), allocatable ::  iterHistory

    real(8)                              ::  dGEMFunctionNorm,    dGEMFunctionNormLast
    real(8), dimension(:),   allocatable ::  dSumMolFractionSoln, dMolesPhaseLast, dUpdateVar, dDrivingForceSoln
    real(8), dimension(:),   allocatable ::  dPartialExcessGibbs, dPartialExcessGibbsLast
    real(8), dimension(:,:), allocatable ::  dEffStoichSolnPhase

    logical                              ::  lDebugMode, lRevertSystem, lConverged
    logical, dimension(:),   allocatable ::  lSolnPhases, lMiscibility

end module ModuleGEMSolver

max_stars_repo_name: nrc-fuels/thermochimica

max_issues_repo_name: nrc-fuels/thermochimica

max_forks_repo_name: nrc-fuels/thermochimica


text:
	program ungrasp
C************************************************************************
C									*
C UNGRASP -                                                             *
C									*
C	This program will take the information from a GRASP surface     *
C file and convert it to a form usable by the render program in         *
C Raster3D.  Note that the 2nd + 3rd records list the file contents,    *
C but I haven't added code here to interpret all of them.  The program  *
C can currently handle potentials and vertexcolors.                     *
C									*
C 30 Sep 1994	Ethan A Merritt - initial version                       *
C 24 Jan 1995	EAM - option to output grid rather than surface		*
C 16 Feb 1996	EAM - fixed scale conversion of potential to color info *
C		      coordinated with Raster3d Version 2.2             *
C		      Recognize GRASP Version 1.2 files (format2) also	*
C 18 Mar 1998	EAM - updated code to use explicit vertex colors in	*
C		      Raster3D V2.4d					*
C  7 May 1998	EAM - ditto for calculated potentials			*
C									*
C************************************************************************
C
	parameter (MAXPTS = 200000)
C
	real*4	vertex(3,MAXPTS)
	real*4	normal(3,MAXPTS)
	real*4	access(3,MAXPTS)
	real*4	potent(MAXPTS)
	real*4	curvature(MAXPTS)
	real*4	distances(MAXPTS)
	real*4	gprop1(MAXPTS)
	real*4	gprop2(MAXPTS)
	real*4	discol(MAXPTS)
	integer*2	triangl2(3,MAXPTS)
	integer*4	triangle(3,MAXPTS)
C
	integer		i,j,k,ic
	real*4		xlo,xhi,ylo,yhi,zlo,zhi
	logical		clip
	logical		grid
	logical		pots, potcol
	logical		cols, vercol
	character*80	line, line3
	character*64	infile, outfile
	integer*4	nvert, ntriangles, gridsize
	real*4		lattice, center(3)
	integer		filetype
C
	integer		indexc(256)
	real*4		coltab(256,3)
	real*4		rgb(3,3)
C

	type 1,'Input file: '
1	format($,a)
	accept '(a)', infile
	open (unit=1, file=infile, status='old', form='unformatted')

c
c	Can we read this file?
	read (1) line
	type *,'header line 1:',line
	if (line(1:8).eq."format=1") then
	    filetype = 1
	else if (line(1:8).eq."format=2") then
	    filetype = 2
	else
	    type *,'Sorry, unrecognized file format'
	    call exit
	endif
c
	read (1) line
	type *,'header line 2:',line
c
	read (1) line3
	type *,'header line 3:',line3

	read (1) line
C	type *,'header line 4:',line
	read (line,*) nvert, ntriangles, gridsize, lattice
	type *, 'vertices, triangles, gridsize, lattice spacing'
	type *,  nvert, ntriangles, gridsize, lattice

	read (1) line
C	type *,'header line 5:',line
	read (line,*) center
	type *, 'center'
	type *, center(1), center(2), center(3)
C
C	Interpret line2 to find contents of file
C
	read (1) (vertex(1,i),vertex(2,i),vertex(3,i), i=1,nvert)
	read (1) (access(1,i),access(2,i),access(3,i), i=1,nvert)
	read (1) (normal(1,i),normal(2,i),normal(3,i), i=1,nvert)
	if (filetype.eq.1) then
	    read (1) (triangl2(1,i),triangl2(2,i),triangl2(3,i),
     &               i=1,ntriangles)
	    do i=1,ntriangles
		triangle(1,i) = triangl2(1,i)
		triangle(2,i) = triangl2(2,i)
		triangle(3,i) = triangl2(3,i)
	    enddo
	endif
	if (filetype.eq.2) then
	    read (1) (triangle(1,i),triangle(2,i),triangle(3,i),
     &               i=1,ntriangles)
	endif
C
C	Interpret line3 to find contents of file
C
	pots = .false.
	cols = .false.
	l = 1
   11	continue
	if (line3(l:l).eq." ") then
		goto 12
	else if (line3(l:l).eq.",") then
		l = l+1
	else if (line3(l:l+9).eq."potentials") then
		pots = .true.
		l = l+10
	else if (line3(l:l+10).eq."vertexcolor") then
		cols = .true.
		l = l+11
	else
		type *,'Don''t know how to handle ', line3(l:80)
		goto 12
	endif
	goto 11
   12	continue

	if (pots) then
	    potmin =  999.
	    potmax = -999.
	    type *,'This file contains potentials'
	    read (1) (potent(i), i=1,nvert)
	    do i=1,nvert
		if (potent(i).gt.potmax) potmax=potent(i)
		if (potent(i).lt.potmin) potmin=potent(i)
	    enddo
	    type *,'Potential range:',potmin,potmax
	    type 1,'Color by potential? '
	    accept '(a80)',line
	    if (line(1:1).eq.'y' .or. line(1:1).eq.'Y') then
	    	potcol =  .true.
		type 1,'What did GRASP say the max potential was? '
		accept *,potxxx
		potmax =  potxxx
		potmin = -potxxx
	    else
	    	potcol = .false.
		potxxx = 1.0
	    endif
	    nclip = 0
	    do i=1,nvert
		if (potent(i) .gt.  potxxx) then
		    potent(i) =  potxxx
		    nclip = nclip + 1
		endif
		if (potent(i) .lt. -potxxx) then
		    potent(i) = -potxxx
		    nclip = nclip + 1
		endif
	    enddo
	    type *,nclip,' truncated color values'
	endif

	if (cols) then
	    type *,'This file contains discrete vertex colors'
	    read (1) (discol(i), i=1,nvert)
	    do j=1,256
		indexc(j) = 0
	    enddo
	    do i=1,nvert
		j = discol(i)
		indexc(j) = indexc(j) + 1
	    enddo
	    type 1,'Use vertex colours? '
	    accept '(a80)',line
	    if (line(1:1).eq.'y' .or. line(1:1).eq.'Y') then
	    	vercol =  .true.
	    	do i=1,256
		    if (indexc(i).gt.0) then
			type '($,a,i4,a)','RGB components for color',i,': '
			accept *, coltab(i,1),coltab(i,2),coltab(i,3)
		    endif
	    	enddo
	    endif
	endif

C
C DEBUG - print out first 10 of each
C
  	goto 100
	do i=1,10
	    type 2, i,	(vertex(j,i),j=1,3)
    2	    format(' Vertex',i4,':',t20,3f8.2)
	enddo
	do i=1,10
	    type 4, i,	(normal(j,i),j=1,3)
    4	    format(' Normal',i4,':',t20,3f8.2)
	enddo
	do i=1,10
	    type 6, i, (triangle(j,i),j=1,3)
    6	    format(' Triangle',i4,':',t20,3i8)
	enddo
C
	do i=1,10
	    i1 = triangle(1,i)
	    i2 = triangle(2,i)
	    i3 = triangle(3,i)
	    type 8, i, (vertex(j,i1),j=1,3), (vertex(k,i2),k=1,3),
     *	 	    (vertex(l,i3),l=1,3)
    8	    format(' Triangle',i4,':',t20,3f8.2,2(/,t20,3f8.2))
	enddo
C
	if (cols) type *, 'Vertex colours',(discol(i),i=1,10)
  100	continue
C
C Allow windowing of triangles output, either to select a
C specific subset of the original file contents, or to
C reduce the processing overhead in Raster3D by omitting
C triangles which will be out of the figure anyway.
C
	type 1,'Clipping bounds (xlo,xhi,ylo,yhi,zlo,zhi:'
	accept '(a80)',line
	read (line,*,err=101) xlo,xhi,ylo,yhi,zlo,zhi
  101	continue
	if (xhi.gt.xlo .and. yhi.gt.ylo .and. zhi.gt.zlo) then
	    clip = .TRUE.
	    type 102,xlo,xhi,ylo,yhi,zlo,zhi
  102	    format('Clipping limits will be: ', 6f10.3)
	else
	    clip = .FALSE.
	    type *,'No clipping limits'
	end if
C
C Allow two output modes
C  1)	Surface made up of triangles with explicit normals
C  2)	Grid/mesh made up of thin cylinders along edges of the triangles
C
	type 1,'Output grid rather than surface?'
	accept '(a80)',line
	if (line(1:1).eq.'y' .or. line(1:1).eq.'Y') then
	    grid =  .true.
	else
	    grid = .false.
C
C Grasp seems a bit chaotic about the direction of the normal vectors
C it puts out.  Allow inversion here.
C
	type 1,'Invert normals?'
	accept '(a80)',line
	if (line(1:1).eq.'y' .or. line(1:1).eq.'Y') then
	    invert = -1
	else
	    invert =  1
	end if
	end if
C
C Write out a Raster3D format file consisting of triangles
C and surface normals (obviously needs a new version of render!)
C
	type 1,'Output file: '
	accept '(a)', outfile
	open (unit=2, file=outfile, form='FORMATTED',
     *	      carriagecontrol='LIST', status='UNKNOWN')
C
C Assume constant colour for surface, unless discrete vertex colors are
C in the GRASP file.
C
	red   = 0.855
	green = 0.855
	blue  = 0.855
	radius = 0.018

	do i=1,ntriangles
	    i1 = triangle(1,i)
	    i2 = triangle(2,i)
	    i3 = triangle(3,i)

	    if (clip) then
		if (vertex(1,i1).lt.xlo .or. vertex(1,i1).gt.xhi) goto 200
		if (vertex(2,i1).lt.ylo .or. vertex(2,i1).gt.yhi) goto 200
		if (vertex(3,i1).lt.zlo .or. vertex(3,i1).gt.zhi) goto 200
	    end if

	    if (potcol) then
		call potrgb( potent(i1), potmin, potmax, rgb(1,1) )
		call potrgb( potent(i2), potmin, potmax, rgb(1,2) )
		call potrgb( potent(i3), potmin, potmax, rgb(1,3) )
	    endif

	    if (vercol) then
		j1 = discol(i1)
		j2 = discol(i2)
		j3 = discol(i3)
		red   = (coltab(j1,1) + coltab(j2,1) + coltab(j3,1)) / 3.
		green = (coltab(j1,2) + coltab(j2,2) + coltab(j3,2)) / 3.
		blue  = (coltab(j1,3) + coltab(j2,3) + coltab(j3,3)) / 3.
	    endif

C	    /* Grid of lines along triangle edges */
	    if (grid) then
	    	write (2,10)	(vertex(j,i1),j=1,3), radius,
     *				(vertex(k,i2),k=1,3), radius,
     *				red, green, blue
	    	write (2,10)	(vertex(k,i2),k=1,3), radius,
     *				(vertex(l,i3),l=1,3), radius,
     *				red, green, blue
	    	write (2,10)	(vertex(l,i3),l=1,3), radius,
     *				(vertex(j,i1),j=1,3), radius,
     *				red, green, blue
   10	    	format('3',/, 3f8.3, f8.3, 3f8.3, f8.3, 1x, 3f6.2)

C	    /* Surface made up of triangles with explicit normals */
	    else
	    	write (2,14)	(vertex(j,i1),j=1,3),
     *				(vertex(k,i2),k=1,3),
     *	 	    		(vertex(l,i3),l=1,3),
     *				red, green, blue
   14	    	format('1',/, 3f8.3, 3f8.3, 3f8.3, 1x, 3f6.2)
	    	write (2,16)	(invert * normal(j,i1),j=1,3),
     *				(invert * normal(k,i2),k=1,3),
     *	 	    		(invert * normal(l,i3),l=1,3)
   16	    	format('7',/, 3f8.3, 3f8.3, 3f8.3, 1x, 3f6.2)
		if (vercol) then
		write (2,17) (coltab(j1,ic),ic=1,3),
     *                       (coltab(j2,ic),ic=1,3),
     *                       (coltab(j3,ic),ic=1,3)
		else if (potcol) then
		write (2,17) (rgb(k,1),k=1,3),
     *			     (rgb(k,2),k=1,3),
     *			     (rgb(k,3),k=1,3)
		endif
   17		format('17',/, 3f8.3, 3f8.3, 3f8.3, 1x, 3f6.2)

	    end if

  200	continue
	enddo

C
	end

	subroutine potrgb( potential, potmin, potmax, rgb )
	real*4	potential, potmin, potmax
	real*4	rgb(3)
c
	if (potential.lt.0) then
	    rgb(1) = 0.8 + 0.1*(potential/potmin)
	    rgb(2) = 0.8 - 0.8*(potential/potmin)
	    rgb(3) = 0.8 - 0.8*(potential/potmin)
	else
	    rgb(3) = 0.8 + 0.1*(potential/potmin)
	    rgb(2) = 0.8 - 0.8*(potential/potmax)
	    rgb(1) = 0.8 - 0.8*(potential/potmax)
	endif
	rgb(1) = rgb(1)*rgb(1)
	rgb(2) = rgb(2)*rgb(2)
	rgb(3) = rgb(3)*rgb(3)
	return
	end


text:
!
!
! Copyright (c) 2011 - 2015
!   University of Houston System and UT-Battelle, LLC.
! Copyright (c) 2009 - 2015
!   Silicon Graphics International Corp.  SHMEM is copyrighted
!   by Silicon Graphics International Corp. (SGI) The OpenSHMEM API
!   (shmem) is released by Open Source Software Solutions, Inc., under an
!   agreement with Silicon Graphics International Corp. (SGI).
!
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
! o Redistributions of source code must retain the above copyright notice,
!   this list of conditions and the following disclaimer.
!
! o Redistributions in binary form must reproduce the above copyright
!   notice, this list of conditions and the following disclaimer in the
!   documentation and/or other materials provided with the distribution.
!
! o Neither the name of the University of Houston System, UT-Battelle, LLC
!   nor the names of its contributors may be used to endorse or promote
!   products derived from this software without specific prior written
!   permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!

program test_shmem_lock
  implicit none
  include 'shmem.fh'

  integer*8, save :: L
  integer  , save :: x

  integer         :: me, npes;
  integer         :: slp;
  integer         :: ret_val;
  integer         :: new_val;

  ! Function definitions
  integer                   :: shmem_my_pe, shmem_n_pes
  integer                   :: shmem_test_lock

  call shmem_init()
  me   = shmem_my_pe()
  npes = shmem_n_pes()
  L = 0
  x = 0
  ret_val = -1

  if(npes .gt. 1) then

    slp = 1
    call shmem_barrier_all()

    call shmem_set_lock(L)

    call shmem_integer_get(new_val, x, 1, 0)
    new_val = new_val + 1

    call shmem_integer_put(x, new_val, 1, 0)  ! increment x on PE 0
    call shmem_quiet()

    call shmem_clear_lock(L)

    call shmem_barrier_all()

    if(me .eq. 0) then
      if(x .eq. npes) then
        write(*,*) "Test for set, and clear lock: Passed"
      else
        write(*,*) "Test for set, and clear lock: Failed"
      end if
      x=0
    end if

    call shmem_barrier_all()

    ret_val = 1
    do while(ret_val .eq. 1)
      ret_val = shmem_test_lock(L)
    end do

    call shmem_integer_get(new_val, x, 1, 0)
    new_val = new_val + 1

    call shmem_integer_put(x, new_val, 1, 0)  ! increment x on PE 0

    call shmem_quiet()

    call shmem_clear_lock(L)

    call shmem_barrier_all()

    if(me .eq. 0) then
      if(x .eq. npes) then
        write(*,*) "Test for test lock: Passed"
      else
        write(*,*) "Test for test lock: Failed"
      end if

    end if

    call shmem_barrier_all()
  else
    write(*,*) "Number of PEs must be > 1 to test locks, test skipped"
  end if

  call shmem_finalize()

end program test_shmem_lock

max_stars_repo_name: openshmem-org/tests-uh

max_issues_repo_name: openshmem-org/tests-uh

max_forks_repo_name: openshmem-org/tests-uh


text:
module mod_d2q9_type
  !> Summary: Define the specifics of the D2Q9 lattice structure

  use mod_kinds, only: ik, rk
  use mod_lattice_type, only: base_2d_lattice_t

  implicit none
  ! private
  ! public :: d2q9_lattice_t

  type, extends(base_2d_lattice_t) :: d2q9_lattice_t
  contains
    procedure :: initialize => init_d2q9
  end type d2q9_lattice_t

contains
  subroutine init_d2q9(self)
    class(d2q9_lattice_t), intent(inout) :: self

    self%n_weights = 9
    self%name = 'D2Q9'
    self%dimension = 2

    allocate(self%f_eq(0:self%n_weights - 1))
    self%f_eq = 0.0_rk

    allocate(self%w_i(0:self%n_weights - 1))
    self%w_i(0) = 4.0_rk / 9.0_rk
    self%w_i(1:4) = 1.0_rk / 9.0_rk
    self%w_i(5:8) = 1.0_rk / 36.0_rk

    allocate(self%v_hat(0:self%n_weights - 1, self%dimension))
    self%v_hat(0, :) = [0, 0]
    self%v_hat(1, :) = [1, 0]
    self%v_hat(2, :) = [0, 1]
    self%v_hat(3, :) = [-1, 0]
    self%v_hat(4, :) = [0, -1]
    self%v_hat(5, :) = [1, 1]
    self%v_hat(6, :) = [-1, 1]
    self%v_hat(7, :) = [-1, -1]
    self%v_hat(8, :) = [1, -1]

    allocate(self%neighbor_indices(1:self%n_weights - 1))
    self%neighbor_indices(1) = 3
    self%neighbor_indices(2) = 4
    self%neighbor_indices(3) = 1
    self%neighbor_indices(4) = 2
    self%neighbor_indices(5) = 7
    self%neighbor_indices(6) = 8
    self%neighbor_indices(7) = 5
    self%neighbor_indices(8) = 6

  end subroutine init_d2q9

end module mod_d2q9_type

max_stars_repo_name: smillerc/lbm

max_issues_repo_name: smillerc/lbm

max_forks_repo_name: smillerc/lbm


text:
subroutine step1(maxm,meqn,maux,mbc,mx,qold,aux,dx,dt,cflgrid,fm,fp,rp1)
!
!     clawpack routine ...  modified for AMRCLAW
!
!     Take one time step, updating q.
!     On entry, qold gives
!        initial data for this step
!        and is unchanged in this version.
!
!     fm, fp are fluxes to left and right of single cell edge
!     See the flux1 documentation for more information.
!
!     Converted to f90 2012-1-04 (KTM)
!     Modified from step2 to work with 1d 2016-6-6 (BND)

    use amr_module

    implicit none

    external rp1

    ! Arguments
    integer, intent(in) :: maxm,meqn,maux,mbc,mx
    real(kind=8), intent(in) :: dx,dt
    real(kind=8), intent(inout) :: cflgrid
    real(kind=8), intent(inout) :: qold(meqn, 1-mbc:mx+mbc)
    real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc)
    real(kind=8), intent(inout) :: fm(meqn, 1-mbc:mx+mbc)
    real(kind=8), intent(inout) :: fp(meqn,1-mbc:mx+mbc)

    ! Local storage for flux accumulation
    real(kind=8) :: faddm(meqn,1-mbc:maxm+mbc)
    real(kind=8) :: faddp(meqn,1-mbc:maxm+mbc)

    ! Scratch storage for Sweeps and Riemann problems
    real(kind=8) ::  q1d(meqn,1-mbc:maxm+mbc)
    real(kind=8) :: aux1(maux,1-mbc:maxm+mbc)
    real(kind=8) :: dtdx1d(1-mbc:maxm+mbc)
    real(kind=8) :: dtdy1d(1-mbc:maxm+mbc)

    real(kind=8) ::  wave(meqn, mwaves, 1-mbc:maxm+mbc)
    real(kind=8) ::     s(mwaves, 1-mbc:maxm + mbc)
    real(kind=8) ::  amdq(meqn,1-mbc:maxm + mbc)
    real(kind=8) ::  apdq(meqn,1-mbc:maxm + mbc)
    real(kind=8) ::  cqxx(meqn,1-mbc:maxm + mbc)
    real(kind=8) :: bmadq(meqn,1-mbc:maxm + mbc)
    real(kind=8) :: bpadq(meqn,1-mbc:maxm + mbc)

    ! Looping scalar storage
    integer :: i,j,thread_num
    real(kind=8) :: dtdx,cfl1d

    ! Common block storage
    integer :: icom
    real(kind=8) :: dtcom,dxcom,tcom
    common /comxt/ dtcom,dxcom,tcom,icom

    ! Store mesh parameters in common block
    dxcom = dx
    dtcom = dt

    cflgrid = 0.d0
    dtdx = dt/dx

    fm = 0.d0
    fp = 0.d0

    ! ============================================================================
    ! Perform X-Sweep
        ! Copy old q into 1d slice
        q1d(:,1-mbc:mx+mbc) = qold(:,1-mbc:mx+mbc)

        ! Set dtdx slice if a capacity array exists
        if (mcapa > 0)  then
            dtdx1d(1-mbc:mx+mbc) = dtdx / aux(mcapa,1-mbc:mx+mbc)
        else
            dtdx1d = dtdx
        endif

        ! Copy aux array into slices
        if (maux > 0) then
            aux1(:,1-mbc:mx+mbc) = aux(:,1-mbc:mx+mbc)
        endif

        ! Compute modification fadd to fluxes along this slice:
        call flux1(maxm,meqn,maux,mbc,mx,q1d,dtdx1d,aux1, &
                   faddm,faddp,cfl1d,wave,s, &
                   amdq,apdq,cqxx,bmadq,bpadq,rp1)

        cflgrid = max(cflgrid,cfl1d)

        ! Update fluxes
        fm(:,1:mx+1) = fm(:,1:mx+1) + faddm(:,1:mx+1)
        fp(:,1:mx+1) = fp(:,1:mx+1) + faddp(:,1:mx+1)

end subroutine step1

max_stars_repo_name: navravi/amrclaw

max_issues_repo_name: navravi/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw


text:
module m_inputFluidProperty
	implicit none; private

	type, public :: inputFluidPropertyObj
    real(8) :: molecularWeight,temperature,gamma,molarConcentration,waterViscosity,gasViscosity,waterSaturateRate,gasSaturateRate, &
		& rhow0,p0,kw
	contains
		procedure :: InitialSet
	end type

	contains
	subroutine InitialSet(inputFluidProperty,number)
	use stdio
	implicit none
	class(inputFluidPropertyObj), intent(inout) :: inputFluidProperty
	integer, intent(in), optional :: number
  integer :: i

	type(fileObj) :: InputFile
	character(ncardWidth) :: readingCard

	inputFluidProperty%gamma = 0.8
	!per unit salt concentration Viscosity increase rate
	inputFluidProperty%molarConcentration = 58.44d-3
	inputFluidProperty%molecularWeight = 44.01d-3
	inputFluidProperty%rhow0 = 1d3
	inputFluidProperty%p0 = 101325
	inputFluidProperty%kw = 2.2d9

	InputFile%filename="./input/inputFluidProperty.env"
	if(present(number)) then
		InputFile%nunit=number*10
		open(unit=InputFile%nunit,file=InputFile%filename,status="old",action="read")
	else
		open(newunit=InputFile%nunit,file=InputFile%filename,status="old",action="read")
	endif

	do
		call InputFile%readline(readingCard)
		if(InputFile%nerr/=0) exit
		if(trim(readingCard)=="#temperature") then
			call InputFile%readline(readingCard)
			read(readingCard,*) inputFluidProperty%temperature
		endif
		if(trim(readingCard)=="#waterViscosity") then
			call InputFile%readline(readingCard)
			read(readingCard,*) inputFluidProperty%waterViscosity
		endif
		if(trim(readingCard)=="#gasViscosity") then
			call InputFile%readline(readingCard)
			read(readingCard,*) inputFluidProperty%gasViscosity
		endif
		if(trim(readingCard)=="#waterSaturateRate") then
			call InputFile%readline(readingCard)
			read(readingCard,*) inputFluidProperty%waterSaturateRate
		endif
		if(trim(readingCard)=="#gasSaturateRate") then
			call InputFile%readline(readingCard)
			read(readingCard,*) inputFluidProperty%gasSaturateRate
		endif
	enddo
	close(InputFile%nunit)

	return
	end subroutine
end module

max_stars_repo_name: ShunsukeKobayashi/fuzzy-broccoli

max_issues_repo_name: ShunsukeKobayashi/fuzzy-broccoli

max_forks_repo_name: ShunsukeKobayashi/fuzzy-broccoli


text:
      subroutine wmeas

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine reads in wind speed data from file and assigns the
!!    data to HRUs

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    hru_sub(:)  |none          |subbasin in which HRU is located
!!    id1         |julian date   |first day of simulation in current year
!!    ifirstw     |none          |wind speed data search code
!!                               |0 first day of wind speed data located in
!!                               |  file
!!                               |1 first day of wind speed data not
!!                               |  located in file
!!    iwgage(:)   |none          |HRU wind speed gage data code (gage # for
!!                               |wind speed data used in HRU)
!!    iyr         |none          |beginning year of simulation
!!    mrg         |none          |maximum number of rainfall/temp gages
!!    nhru        |none          |number of HRUs in watershed
!!    nwtot       |none          |total number of wind speed records
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ifirstw     |none          |wind speed data search code
!!                               |0 first day of wind speed data located in
!!                               |  file
!!                               |1 first day of wind speed data not
!!                               |  located in file
!!    u10(:)      |m/s           |wind speed for the day in HRU
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    idap        |julian date   |julian date of measured weather data
!!    inum3sprev  |none          |subbasin number of previous HRU
!!    iyp         |none          |last 2 digits of year measured weather data
!!    k           |none          |counter
!!    l           |none          |counter
!!    u10bsb      |m/s           |generated wind speed for subbasin
!!    wndmeas(:)  |m/s           |wind speed read from file
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    SWAT: wndgen

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

      integer :: k, iyp, idap, l, inum3sprev
      real :: u10bsb
      real, dimension (mrg) :: wndmeas

      !! initialize variables for the day
      wndmeas = 0.

        !! read wind speed data from file
          if (ifirstw == 0) then
            read (139,5200) (wndmeas(l), l = 1, nwtot)
          else
            ifirstw = 0
            do
              iyp = 0
              idap = 0
              read (139,5300) iyp, idap, (wndmeas(l), l = 1, nwtot)
              if (iyp + idap <= 0) exit
              if (iyp == iyr .and. idap == id1) exit
            end do
          end if

        !! assign wind speed data to HRUs
        inum3sprev = 0
        do k = 1, nhru
          u10(k) = wndmeas(iwgage(hru_sub(k)))
          !! generate values to replace missing data
          if (u10(k) <  -97.) then
            !! use same generated data for all HRUs in a subbasin
            if (hru_sub(k) == inum3sprev .and. hru_sub(k) /= 0) then
              u10(k) = u10bsb
            else
              call wndgen(k)
              !! set subbasin generated values
              inum3sprev = 0
              u10bsb = 0.
              inum3sprev = hru_sub(k)
              u10bsb = u10(k)
            end if
          end if
        end do

      return
! 5200 format (7x,300f8.3)
! 5300 format (i4,i3,300f8.3)
 5200 format (7x,1800f8.3)
 5300 format (i4,i3,1800f8.3)
      end


text:
! ----------------------------------------------------------------------
      subroutine nescoil (loopcount)
!................................................................
      USE stel_kinds
      USE Vvacuum2, ONLY: iota_edge, phip_edge, curpol
      USE Vvacuum3
      USE Vprecal1, ONLY: np
      USE Vbiopre
      USE Vculine1
      USE Vprecal4
      use SvdCtrl, ONLY : noaccuracy
      use LoopCtrl
      use Vbnorm, ONLY : extension
      use NumParams, ONLY : inesc
      use safe_open_mod
      implicit none
!................................................................
C   D u m m y   A r g u m e n t s
c................................................................
      integer :: loopcount
c................................................................
C   L o c a l   V a r i a b l e s
c................................................................
      integer :: istat
      real(rprec) :: unit, file, t1, t2, tbeg, tend
      character*10 :: date, time

c................................................................
c
c      The NESCOIL code is based on a method described in
c         [1]  P.Merkel,Nucl.Fusion 27 (1987) 867.
c         [2]  P.Merkel, In: Theory of Fusion Plasmas,Eds. A.Bondeson
c              E.Sindoni,and F.Troyon, Varenna,Italy,EUR 11336 EN,25-46.
c
c      Introducing two nested toroidally closed surfaces a current
c      on the outer surface is determined such that the normal
c      component of the magnetic field on the inner surface is
c      minimized.
c
c      The current carrying surface is mapped onto the unit square:
c      0< u <1, 0< v <1.
c
c
c         ^
c       1.|
c         |                  cup - poloidal current per period
c         |                  cut - toroidal current
c         |
c   ^     |
c   |   u |
c   |     |
c         |
c  cup    |
c         |
c         |
c         |
c       0 ----------------------------->
c         0                v           1.
c                     <----- cut
c
c
c
c................................................................
c     PARAMETERS READ IN FROM INPUT FILE:
c
c     1. DIMENSIONS:
c        nu    - number of poloidal meshpoints for coil surface
c        nv    - number of toroidal meshpoints for coil surface
c
c        nu1   - number of poloidal meshpoints for plasma surface
c        nv1   - number of toroidal meshpoints for plasma surface
c
c        mf    - number of poloidal modes for current potential
c        nf    - number of toroidal modes for current potential
c
c        md    - number of poloidal modes for plasma and coil surface shape
c        nd    - number of toroidal modes for plasma and coil surface shape
c
c        npol  - number of segments of a modular or helical filament
c        ntor  - number of filaments per period
c
c        np    - number of field periods
c
c     2. NESCOIL CONTROL PARAMETERS:
c        cut   - Net toroidal current +-1 or 0
c        cup   - Net poloidal current +-1 or 0
c
c        ibex  - An  external field can be added by supplying
c                a subroutine bexter with IBEX = 1
c
c     3. INFORMATION FROM VMEC PLASMA:
c        iota_edge - iota at plasma edge
c        phip_edge - phi_prime at plasma edge
c        curpol - Total poloidal current Amps per field period from VMEC
c
c     4. SVD CONTROL PARAMETERS:
c        1. MSTRT = Method + svdscan start if >1
c                 >=0:Berr, <0:Xerr, unless MSTEP <=0: Least square
c        2. MSTEP = Method + svdscan stepsize:
c                 <=0: LeastSquare, =0: use old f04abe (now SOLVER), no svd
c        3. MKEEP = svd/scan control: =0: svdscan, else keep |nkeep| wgts
c                 <0: write all weights to output
c        4. MDSPW = 2+exponent of dsur multiplying bfn,ben:
c                 <0: post-calculate Xerr svdscan and write to output
c        5. CURWT = Weight for surface current minimization
c                   Works ONLY in LSQ branch
c        6. TRGWT : Not implemented yet (PMV)

c................................................................
c        Secondary parameters:
c
c        md,nd - dimension of arrays: md>=mf,nd>=nf
c        nmax  - dimension of array:  nmax> 3*(nu+2)*nv
c
c................................................................
c      Set loop counter iloop for internal use
      iloop = loopcount

c     Strategy based on iloop:
c      1. iloop = 0 or  1 : single run or first call, allocate
c      2. iloop = 0 or -1 : single run or last  call, deallocate
c................................................................
c     Read the input file only in the first loop
c     Read input data from file specified on command line:
c     xnesopt input
c     NOTE: SINCE ALL DIMENSIONS ARE ALLOCATABLE NOW,
c           YOU MUST FIRST READ INPUT FILE (TO GET DIMENSIONS)
c           BEFORE DOING ANYTHING ELSE
      if (iloop == 0 .or. iloop == 1) call read_nescoil_input
c     This ends first-loop tasks, rest is done in each loop
c................................................................
c
      write (inesc, '(a,i3,a)') '-----  Begin Nescoil run num ', iloop,
     1    ' -----'

c     Now that all inputs are here, begin nescoil run
      call date_and_time(date,time)
      write (inesc, 100) date(5:6),date(7:8),date(1:4),
     1  time(1:2),time(3:4),time(5:6)
 100  format('DATE = ',a2,'-',a2,'-',a4,' ',' TIME = ',2(a2,':'),a2)
      call second0(tbeg)   !This one for timing the whole run
c
c     setup constants, arrays, and allocations based on nwire
c     uncomment for NT      USE NUMERICAL_LIBRARIES
c
      nw = ntor*np*(npol+1)
      istat=0
      if (.not.allocated(vx))
     1 allocate (vx(nw), vy(nw), vz(nw), dx(nw), dy(nw), dz(nw),
     1           xw(nw), yw(nw), zw(nw), curre(nw),
     2           cok(0:np - 1), sik(0:np - 1), stat=istat)
      curre = 0; xw = 0; yw = 0; zw = 0                                !!MUST be initialized (SPH)
      if (istat .ne. 0) stop 'allocation error nw in NESCOIL'

c     Do all precalculations
      call precal
      call second0(t2)
      write (inesc,"('PRECAL took ',g12.3,' sec')") t2-tbeg

c     Compute quantities on plasma surface
      write (inesc, 10) '----- Calling Surface_Plasma -----'
      call second0(t1)
      call surface_plas
      call second0(t2)
      write (inesc,"('Time in Surface_Plasma: ',g12.3,' sec')") t2-t1

c     Compute quantities on coil surface
      write (inesc, 10) '----- Calling Surface_Coil -----'
      call second0(t1)
      call surface_coil
      call second0(t2)
      write (inesc,"('Time in Surface_Coil: ',g12.3,' sec')") t2-t1

c................................................................
c     Solve boundary value problem
      write (inesc, 10) '----- Calling Solver -----'
      call second0(t1)
      call solver_nescoil
      call second0(t2)
      write (inesc,"('Time in Solver: ',g12.3,' sec')") t2-t1

c................................................................
c     Post-process solution for various answers
      write (inesc, 10) '----- Calling Surfcur_Diag -----'
      call second0(t1)
      call surfcur_diag
      call second0(t2)
      write (inesc,"('Time in Surfcur_Diag: ',g12.3,' sec')") t2-t1

      if(noaccuracy .eq. 0) then
         write (inesc, 10) '----- Calling Accuracy -----'
         call second0(t1)
         call accuracy
         call second0(t2)
         write (inesc,"('ACCURACY took ',g12.3,' sec')") t2-t1
      endif

      call second0(tend)
      write (inesc,"('ONE NESCOIL RUN took ',g12.3,' sec')")
     1    tend-tbeg

 10   format(a)
c................................................................
c     Tasks to be done in last loop only
      if (iloop .le. 0) call nescoil_cleanup
c................................................................

      end subroutine nescoil

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
C    @(#)p_ctlpow.f	20.5 8/20/98
C****************************************************************
C
C   File: p_ctlpow.f
C   Purpose: IPF shell program to process miscellaneous commands
C
C   Author: Walt Powell  Date: 29 November 1992
C                        Modified:
C   Called by: pf_proc, p_commnd
C
C****************************************************************
C
	integer function p_ctlpow (in_buffer, out_buffer)

        character in_buffer*(*)
        character out_buffer*(*)

        include 'ipfinc/parametr.inc'
        include 'ipfinc/blank.inc'
        include 'ipfinc/lfiles.inc'
        include 'ipfinc/filnam.inc'
        include 'ipfinc/jobctl.inc'
        include 'ipfinc/dtaiop.inc'
        include 'ipfinc/changr.inc'
        include 'ipfinc/delete.inc'
        include 'ipfinc/epridc.inc'
        include 'ipfinc/errorsw.inc'
        include 'ipfinc/pfstates.inc'
        include 'ipfinc/prt.inc'
        include 'ipfinc/errmsg.inc'
        include 'ipfinc/errorx.inc'

        character null * 1, linefeed * 1, text * 80
        integer o2, apdoutbuf

        null = char(0)
        linefeed = char(10)
        p_ctlpow = 0     ! default return SUCCESS state
        numerr = 0       ! reinitialize error count
        out_buffer(1:1) = null

        buf = inrcd
        call ctlpow
        buf = inrcd
        call loadarcv
c
c       Append error messages to buffer
c
        j = 1
        length = 1
        o2 = index (out_buffer,null)
        do while (j .le. numerr .and. length .gt. 0)
           length = apdoutbuf(o2, errm(j), out_buffer(o2:))
           o2 = o2 + length
           j = j + 1
        enddo
c
c 	Append summary
c
        write (text, 340) 'p_ctlpow.f', p_ctlpow, ostates
  340   format ('/', a, ' return status: ', i2, ' IPF state: ', i2)
        length = apdoutbuf(o2, text, out_buffer(o2:))
        o2 = o2 + length
c
c       Reset error flag
c
        call setercnt (0, ' ')
        return
	end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
SUBMODULE(Times) testTimesSubmodule
#include "UnitTest.h"
USE UnitTest
USE IntrType
USE Times

IMPLICIT NONE

TYPE(TimerType) :: testTimer

INTEGER :: idum1,idum2,idum3,ioerr
CHARACTER(LEN=1) :: adum1,adum2
CHARACTER(LEN=5) :: adum3
CHARACTER(LEN=2) :: adum4
CHARACTER(LEN=MAXLEN_DATE_STRING) :: adate
CHARACTER(LEN=MAXLEN_CLOCK_STRING) :: aclock
REAL(SRK) :: totalElapsed
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
SUBROUTINE runTestTimes()

  CREATE_TEST('Timers')

  REGISTER_SUBTEST('Basic Timers',testTimers)
  REGISTER_SUBTEST('Parent Timers',testParentTimers)

  FINALIZE_TEST()

ENDSUBROUTINE runTestTimes
!
!-------------------------------------------------------------------------------
SUBROUTINE testTimers()

  COMPONENT_TEST('getDate()')
  !
  !Test getDate()
  adate=getDate()
  idum1=0
  idum2=0
  idum3=0
  READ(adate,'(i2,a1,i2,a1,i4)',iostat=ioerr) idum1,adum1,idum2,adum2,idum3
  ASSERT_EQ(ioerr,0,'ioerr')
  ASSERT_GT(idum1,0,'month')
  ASSERT_LT(idum1,13,'month')
  ASSERT_GT(idum2,0,'day')
  ASSERT_LT(idum2,32,'day')
  ASSERT_GT(idum3,0,'year')
  ASSERT_EQ(adum1,adum2,'adum')
  ASSERT_EQ(adum1,'/','adum1')
  INFO(0) 'getDate() = '//getDate()
  ASSERT_EQ(getDate(),getDate(1),'getDate(1)')

  idum1=0
  idum2=0
  adate=getDate(2)
  READ(adate,'(a5,i2,a2,i4)',iostat=ioerr) adum3,idum1,adum4,idum2
  ASSERT_EQ(ioerr,0,'ioerr')
  ASSERT_GT(idum1,0,'day')
  ASSERT_LT(idum1,32,'day')
  ASSERT_GT(idum2,0,'year')
  ASSERT_EQ(adum4(1:1),',','month')

  !Test getTimeFromDate
  COMPONENT_TEST('getTimeFromDate')
  ASSERT_APPROXEQA(getTimeFromDate('12/01/1990','12/02/1990'),1.0_SRK,'check Defaults')
  ASSERT_APPROXEQA(getTimeFromDate('12/1/1990','12/2/1990','HOUR'),24.0_SRK,'check HOUR and MM/D/YYYY fmt')
  ASSERT_APPROXEQA(getTimeFromDate('1/1/1990','1/2/1990','MIN'),1440.0_SRK,'check MIN and M/D/YYYY fmt')
  ASSERT_APPROXEQA(getTimeFromDate('1/10/1990','1/11/1990','SEC'),86400.0_SRK,'check SEC and M/DD/YYYY fmt')
  ASSERT_APPROXEQA(getTimeFromDate('1990/12/01','1994/12/01','DAY'),1461.0_SRK,'check DAY, YYYY/MM/DD fmt and leapyear calls')
  ASSERT_APPROXEQA(getTimeFromDate('1990/12/1','1998/12/1','DAY'),2922.0_SRK,'check DAY, YYYY/MM/D fmt and leapyear calls')
  ASSERT_APPROXEQA(getTimeFromDate('1990/1/01','2002/1/01','DAY'),4383.0_SRK,'check DAY, YYYY/M/DD fmt and leapyear calls')
  ASSERT_APPROXEQA(getTimeFromDate('1890/1/1','1906/1/1','DAY'),5843.0_SRK,'check DAY, YYYY/M/D fmt and leapyear calls')
  ASSERT_APPROXEQA(getTimeFromDate('09/05/1997','10/05/1997','SEC'),2592000.0_SRK,'check SEC 09/05/1997 and 10/05/1997')
  ASSERT_APPROXEQA(getTimeFromDate('9/20/2009','10/19/2009','HOUR'),696.0_SRK,'check HOUR 9/20/2009 and 10/19/2009')

  !Test getClockTime
  COMPONENT_TEST('getClockTime()')
  idum1=0
  idum2=0
  idum3=0
  aclock=getClockTime()
  READ(aclock,'(i2,a1,i2,a1,i2)',iostat=ioerr) idum1,adum1,idum2,adum2,idum3
  ASSERT_EQ(ioerr,0,'ioerr')
  ASSERT_GT(idum1,-1,'hour')
  ASSERT_LT(idum1,24,'hour')
  ASSERT_GT(idum2,-1,'minute')
  ASSERT_LT(idum2,60,'minute')
  ASSERT_GT(idum3,-1,'minute')
  ASSERT_LT(idum3,60,'minute')
  ASSERT_EQ(adum1,':',':')
  ASSERT_EQ(adum1,adum2,':')

  COMPONENT_TEST('HI-RES TIMER')
  ASSERT_EQ(LEN_TRIM(testTimer%getTimerName()),0,'%getTimerName()')
  CALL testTimer%setTimerName('myName')
  ASSERT_EQ(TRIM(testTimer%getTimerName()),'myName','%setTimerName()')
  testTimer%elapsedtime=0.0001_SRK
  ASSERT_EQ(testTimer%getTimeReal(),0.0001_SRK,'%getTimeReal()')
  ASSERT_EQ(testTimer%getTimeChar(),'  100.000 microsec','%getTimeChar() (us)')
  ASSERT_EQ(testTimer%getTimeHHMMSS(),'00:00.00   ','%getTimeHHMMSS() (us)')
  testTimer%elapsedtime=0.999_SRK
  ASSERT_EQ(testTimer%getTimeChar(),'  999.000 ms      ','%getTimeChar() (ms)')
  ASSERT_EQ(testTimer%getTimeHHMMSS(),'00:01.00   ','%getTimeHHMMSS() (ms)')
  testTimer%elapsedtime=100.637_SRK
  ASSERT_EQ(testTimer%getTimeChar(),'  100.637 s       ','%getTimeChar() (s)')
  ASSERT_EQ(testTimer%getTimeHHMMSS(),'01:40.64   ','%getTimeHHMMSS() (s)')
  testTimer%elapsedtime=100000.6_SRK
  ASSERT_EQ(testTimer%getTimeChar(),'27:46:40.60 hh:mm:ss','%getTimeChar() (hr)')
  ASSERT_EQ(testTimer%getTimeHHMMSS(),'27:46:40.60','%getTimeHHMMSS() (hr)')
  testTimer%elapsedtime=3719.6_SRK
  ASSERT_EQ(testTimer%getTimeHHMMSS(),'01:01:59.60','%getTimeHHMMSS() (round)')
  CALL testTimer%ResetTimer()
  ASSERT_EQ(LEN_TRIM(testTimer%getTimername()),0,'name')
  ASSERT_EQ(testTimer%elapsedtime,0._SRK,'elapsedtime')
  ASSERT_EQ(testTimer%getTimeChar(),'    0.000 microsec','time char')
  CALL testTimer%tic()
  CALL sleep(1)
  CALL testTimer%toc()
  totalElapsed=testTimer%elapsedtime
  ASSERT_GT(totalElapsed,0.0_SRK,'tic/toc')
  CALL testTimer%tic()
  CALL sleep(1)
  CALL testTimer%toc()
  ASSERT_GT(testTimer%elapsedtime,totalElapsed,'tic/toc')
  totalElapsed=testTimer%elapsedtime
  CALL testTimer%tic()
  CALL sleep(1)
  CALL testTimer%toc()
  ASSERT_GT(testTimer%elapsedtime,totalElapsed,'tic/toc')
  totalElapsed=testTimer%elapsedtime
  ASSERT_LT(totalElapsed,60.0_SRK,'tic/toc too slow')
  ASSERT(testTimer%getTimerHiResMode(),'%getTimerHiResMode()')

  COMPONENT_TEST('LO-RES TIMER')
  CALL testTimer%setTimerHiResMode(.FALSE.)
  ASSERT(.NOT.testTimer%getTimerHiResMode(),'%getTimerHiResMode()')
  adate=testTimer%getDate()
  idum1=0
  idum2=0
  idum3=0
  READ(adate,'(i2,a1,i2,a1,i4)',iostat=ioerr) idum1,adum1,idum2,adum2,idum3
  ASSERT_EQ(ioerr,0,'ioerr')
  ASSERT_GT(idum1,0,'month')
  ASSERT_LT(idum1,13,'month')
  ASSERT_GT(idum2,0,'day')
  ASSERT_LT(idum2,32,'day')
  ASSERT_GT(idum3,0,'year')
  ASSERT_EQ(adum1,adum2,'adum')
  ASSERT_EQ(adum1,'/','adum1')
  ASSERT_EQ(testTimer%getDate(),testTimer%getDate(1),'%getDate(1)')
  idum1=0
  idum2=0
  adate=testTimer%getDate(2)
  READ(adate,'(a5,i2,a2,i4)',iostat=ioerr) adum3,idum1,adum4,idum2
  ASSERT_EQ(ioerr,0,'ioerr')
  ASSERT_GT(idum1,0,'day')
  ASSERT_LT(idum1,32,'day')
  ASSERT(idum2 > 0,'year')
  ASSERT_EQ(adum4(1:1),',','month')
  idum1=0
  idum2=0
  idum3=0
  aclock=testTimer%getClockTime()
  READ(aclock,'(i2,a1,i2,a1,i2)',iostat=ioerr) idum1,adum1,idum2,adum2,idum3
  ASSERT_EQ(ioerr,0,'ioerr')
  ASSERT_GT(idum1,-1,'hour')
  ASSERT_LT(idum1,24,'hour')
  ASSERT_GT(idum2,-1,'minute')
  ASSERT_LT(idum2,60,'minute')
  ASSERT_GT(idum3,-1,'minute')
  ASSERT_LT(idum3,60,'minute')
  ASSERT_EQ(adum1,':',':')
  ASSERT_EQ(adum1,adum2,':')
  testTimer%elapsedtime=0.0001_SRK

  COMPONENT_TEST('testTimer%getRemainingTime')
  CALL testTimer%tic()
  CALL sleep(1)
  CALL testTimer%toc()
  totalElapsed=testTimer%elapsedtime
  ASSERT_GT(totalElapsed,0.0_SRK,'tic/toc')
  CALL testTimer%tic()
  CALL sleep(1)
  CALL testTimer%toc()
  ASSERT_GT(testTimer%elapsedtime,totalElapsed,'tic/toc')
  totalElapsed=testTimer%elapsedtime
  CALL testTimer%tic()
  CALL sleep(1)
  CALL testTimer%toc()
  ASSERT_GT(testTimer%elapsedtime,totalElapsed,'tic/toc')
  totalElapsed=testTimer%elapsedtime
  ASSERT_LT(totalElapsed,60.0_SRK,'tic/toc too slow')

ENDSUBROUTINE testTimers
!
!-------------------------------------------------------------------------------
SUBROUTINE testParentTimers()
  TYPE(ParentTimerType) :: parentTimer
  TYPE(TimerType),POINTER :: tmpptr
  CLASS(TimerType),POINTER :: timerptr

  COMPONENT_TEST('clear')
  CALL parentTimer%clear()
  ALLOCATE(parentTimer%timers(1))
  !Doing this tmpptr thing because of what seems like a compiler bug.  Allocating
  !timers(1)%t directly causes a segfault unless I have a WRITE statement first.
  !It only happens in this test, not elsewhere
  tmpptr => NULL()
  ALLOCATE(tmpptr)
  parentTimer%timers(1)%t => tmpptr
  parentTimer%name='testName'
  CALL parentTimer%clear()
  ASSERT(.NOT.ALLOCATED(parentTimer%timers),'ALLOCATED %timers')
  ASSERT_EQ(parentTimer%getTimerName(),'','%name')

  COMPONENT_TEST('addTimer')
  CALL parentTimer%addTimer('test1')
  ASSERT_EQ(SIZE(parentTimer%timers),1,'SIZE %timers')
  SELECTTYPE(t => parentTimer%timers(1)%t)
  TYPE IS(TimerType)
    ASSERT(.TRUE.,'timer(1) type')
  CLASS DEFAULT
    ASSERT(.FALSE.,'timer(1) type')
  ENDSELECT
  CALL parentTimer%addTimer('test2')
  ASSERT_EQ(SIZE(parentTimer%timers),2,'SIZE %timers')
  SELECTTYPE(t => parentTimer%timers(2)%t)
  TYPE IS(TimerType)
    ASSERT(.TRUE.,'timer(2) type')
  CLASS DEFAULT
    ASSERT(.FALSE.,'timer(2) type')
  ENDSELECT
  CALL parentTimer%addTimer('test3 -> subtest1')
  ASSERT_EQ(SIZE(parentTimer%timers),3,'SIZE %timers')
  SELECTTYPE(t => parentTimer%timers(3)%t)
  TYPE IS(ParentTimerType)
    ASSERT(.TRUE.,'timer(3) type')
    ASSERT_EQ(SIZE(t%timers),1,'SIZE %timers(3)%timers')
  CLASS DEFAULT
    ASSERT(.FALSE.,'timer(3) type')
  ENDSELECT
  CALL parentTimer%addTimer('test4 -> subtest2')
  ASSERT_EQ(SIZE(parentTimer%timers),4,'SIZE %timers')
  SELECTTYPE(t => parentTimer%timers(4)%t)
  TYPE IS(ParentTimerType)
    ASSERT(.TRUE.,'timer(4) type')
    ASSERT_EQ(SIZE(t%timers),1,'SIZE %timers(4)%timers')
  CLASS DEFAULT
    ASSERT(.FALSE.,'timer(4) type')
  ENDSELECT
  CALL parentTimer%addTimer('subtest1')
  ASSERT_EQ(SIZE(parentTimer%timers),5,'SIZE %timers')
  SELECTTYPE(t => parentTimer%timers(1)%t)
  TYPE IS(TimerType)
    ASSERT(.TRUE.,'timer(5) type')
  CLASS DEFAULT
    ASSERT(.FALSE.,'timer(5) type')
  ENDSELECT
  CALL parentTimer%addTimer('test4 -> subtest1')
  ASSERT_EQ(SIZE(parentTimer%timers),5,'SIZE %timers')
  SELECTTYPE(t => parentTimer%timers(4)%t)
  TYPE IS(ParentTimerType)
    ASSERT(.TRUE.,'timer(4) type')
    ASSERT_EQ(SIZE(t%timers),2,'SIZE %timers(4)%timers')
  CLASS DEFAULT
    ASSERT(.FALSE.,'timer(4) type')
  ENDSELECT
  CALL parentTimer%addTimer('test4 -> test1')
  ASSERT_EQ(SIZE(parentTimer%timers),5,'SIZE %timers')
  SELECTTYPE(t => parentTimer%timers(4)%t)
  TYPE IS(ParentTimerType)
    ASSERT(.TRUE.,'timer(4) type')
    ASSERT_EQ(SIZE(t%timers),3,'SIZE %timers(4)%timers')
  CLASS DEFAULT
    ASSERT(.FALSE.,'timer(4) type')
  ENDSELECT
  CALL parentTimer%addTimer('test2')
  ASSERT_EQ(SIZE(parentTimer%timers),5,'SIZE %timers')
  CALL parentTimer%addTimer('test4 -> test1')
  ASSERT_EQ(SIZE(parentTimer%timers),5,'SIZE %timers')
  SELECTTYPE(t => parentTimer%timers(4)%t)
  TYPE IS(ParentTimerType)
    ASSERT(.TRUE.,'timer(4) type')
    ASSERT_EQ(SIZE(t%timers),3,'SIZE %timers(4)%timers')
  CLASS DEFAULT
    ASSERT(.FALSE.,'timer(4) type')
  ENDSELECT

  COMPONENT_TEST('getTimer')
  timerptr => parentTimer%getTimer('test1')
  ASSERT(ASSOCIATED(timerptr),'ASSOCIATED %getTimer("test1")')
  ASSERT_EQ(CHAR(timerptr%name),'test1','getTimer("test1")')
  timerptr => parentTimer%getTimer('test2')
  ASSERT(ASSOCIATED(timerptr),'ASSOCIATED %getTimer("test2")')
  ASSERT_EQ(CHAR(timerptr%name),'test2','getTimer("test2")')
  timerptr => parentTimer%getTimer('test3 -> subtest1')
  ASSERT(ASSOCIATED(timerptr),'ASSOCIATED %getTimer("test3 -> subtest1")')
  ASSERT_EQ(CHAR(timerptr%name),'subtest1','getTimer("test3 -> subtest1")')
  timerptr => parentTimer%getTimer('test4 -> subtest2')
  ASSERT(ASSOCIATED(timerptr),'ASSOCIATED %getTimer("test4 -> subtest1")')
  ASSERT_EQ(CHAR(timerptr%name),'subtest2','getTimer("test4 -> subtest2")')
  timerptr => parentTimer%getTimer('subtest1')
  ASSERT(ASSOCIATED(timerptr),'ASSOCIATED %getTimer("subtest1")')
  ASSERT_EQ(CHAR(timerptr%name),'subtest1','getTimer("subtest1")')
  timerptr => parentTimer%getTimer('test4 -> subtest1')
  ASSERT(ASSOCIATED(timerptr),'ASSOCIATED %getTimer("test4 -> subtest1")')
  ASSERT_EQ(CHAR(timerptr%name),'subtest1','getTimer("test4 -> subtest1")')
  timerptr => parentTimer%getTimer('test4 -> test1')
  ASSERT(ASSOCIATED(timerptr),'ASSOCIATED %getTimer("test4 -> test1")')
  ASSERT_EQ(CHAR(timerptr%name),'test1','getTimer("test4 -> test1")')
  timerptr => parentTimer%getTimer('subtest2')
  ASSERT(.NOT.ASSOCIATED(timerptr),'ASSOCIATED %getTimer("subtest2")')
  timerptr => parentTimer%getTimer('test4 -> subtest1')
  ASSERT(ASSOCIATED(timerptr),'ASSOCIATED %getTimer("test4 -> subtest1")')
  ASSERT_EQ(CHAR(timerptr%name),'subtest1','getTimer("test4 -> subtest1")')
  timerptr => parentTimer%getTimer('test')
  ASSERT(.NOT.ASSOCIATED(timerptr),'ASSOCIATED %getTimer("test")')

ENDSUBROUTINE testParentTimers
!
!
!-------------------------------------------------------------------------------
ENDSUBMODULE testTimesSubmodule
PROGRAM testTimes
  USE Times
  CALL runTestTimes()
ENDPROGRAM testTimes

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: picmc/Futility