text:
cBHEADER**********************************************************************
c Copyright (c) 2008,  Lawrence Livermore National Security, LLC.
c Produced at the Lawrence Livermore National Laboratory.
c This file is part of HYPRE.  See file COPYRIGHT for details.
c
c HYPRE is free software; you can redistribute it and/or modify it under the
c terms of the GNU Lesser General Public License (as published by the Free
c Software Foundation) version 2.1 dated February 1999.
c
c $Revision: 1.3 $
cEHEADER**********************************************************************

c*****************************************************************************
c HYPRE_IJMatrix Fortran interface
c*****************************************************************************

c--------------------------------------------------------------------------
c HYPRE_IJMatrixCreate
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixcreate(fcomm, filower, fiupper,
     1                                 fjlower, fjupper, fmatrix)

      integer ierr
      integer*8 fmatrix
      integer*8 fcomm
      integer filower
      integer fiupper
      integer fjlower
      integer fjupper

      call HYPRE_IJMatrixCreate(fcomm, filower, fiupper, fjlower,
     1                          fjupper, fmatrix, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixcreate error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixDestroy
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixdestroy(fmatrix)

      integer ierr
      integer*8 fmatrix

      call HYPRE_IJMatrixDestroy(fmatrix, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixdestroy error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixInitialize
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixinitialize(fmatrix)

      integer ierr
      integer*8 fmatrix

      call HYPRE_IJMatrixInitialize(fmatrix, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixinitialize error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixAssemble
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixassemble(fmatrix)

      integer ierr
      integer*8 fmatrix

      call HYPRE_IJMatrixAssemble(fmatrix, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixassemble error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixSetRowSizes
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixsetrowsizes(fmatrix, fizes)

      integer ierr
      integer*8 fmatrix
      integer fsizes

      call HYPRE_IJMatrixSetRowSizes(fmatrix, fsizes, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixsetrowsizes error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixSetDiagOffdSizes
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixsetdiagoffdsizes(fmatrix, fdiag_sizes,
     1                                           foffd_sizes)

      integer ierr
      integer*8 fmatrix
      integer fdiag_sizes
      integer foffd_sizes

      call HYPRE_IJMatrixSetDiagOffdSizes(fmatrix, fdiag_sizes,
     1                                    foffd_sizes, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixsetdiagoffdsizes error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixSetValues
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixsetvalues(fmatrix, fnrows, fncols,
     1                                    frows, fcols, fvalues)

      integer ierr
      integer*8 fmatrix
      integer fnrows
      integer fncols
      integer frows
      integer fcols
      double precision fvalues

      call HYPRE_IJMatrixSetValues(fmatrix, fnrows, fncols, frows,
     1                             fcols, fvalues, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixsetvalues error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixAddToValues
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixaddtovalues(fmatrix, fnrows, fncols,
     1                                      frows, fcols, fvalues)

      integer ierr
      integer*8 fmatrix
      integer fnrows
      integer fncols
      integer frows
      integer fcols
      double precision fvalues

      call HYPRE_IJMatrixAddToValues(fmatrix, fnrows, fncols, frows,
     1                               fcols, fvalues, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixaddtovalues error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixSetObjectType
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixsetobjecttype(fmatrix, ftype)

      integer ierr
      integer*8 fmatrix
      integer ftype

      call HYPRE_IJMatrixSetObjectType(fmatrix, ftype, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixsetobjecttype error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixGetObjectType
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixgetobjecttype(fmatrix, ftype)

      integer ierr
      integer*8 fmatrix
      integer ftype

      call HYPRE_IJMatrixGetObjectType(fmatrix, ftype, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixgetobjecttype error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixGetObject
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixgetobject(fmatrix, fobject)

      integer ierr
      integer*8 fmatrix
      integer*8 fobject

      call HYPRE_IJMatrixGetObject(fmatrix, fobject, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixgetobject error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixRead
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixread(ffilename, fcomm, fobject_type,
     1                               fmatrix)

      integer ierr
      integer*8 fmatrix
      integer*8 fcomm
      integer fobject_type
      character*(*) ffilename

      call HYPRE_IJMatrixRead(ffilename, fcomm, fobject_type, fmatrix,
     1                        ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixread error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJMatrixPrint
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixprint(fmatrix, ffilename)

      integer ierr
      integer*8 fmatrix
      character*(*) ffilename

      call HYPRE_IJMatrixPrint(fmatrix, ffilename, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixprint error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c hypre_IJMatrixSetObject
c--------------------------------------------------------------------------
      subroutine fhypre_ijmatrixsetobject(fmatrix, fobject)

      integer ierr
      integer*8 fmatrix
      integer*8 fobject

      call hypre_IJMatrixSetObject(fmatrix, fobject, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijmatrixsetobject error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorCreate
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorcreate(fcomm, fjlower, fjupper, fvector)

      integer ierr
      integer*8 fvector
      integer fcomm
      integer fjlower
      integer fjupper

      call HYPRE_IJVectorCreate(fcomm, fjlower, fjupper, fvector, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorcreate error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorDestroy
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectordestroy(fvector)

      integer ierr
      integer*8 fvector

      call HYPRE_IJVectorDestroy(fvector, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectordestroy error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorInitialize
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorinitialize(fvector)

      integer ierr
      integer*8 fvector

      call HYPRE_IJVectorInitialize(fvector, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorinitialize error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorSetValues
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorsetvalues(fvector, fnum_values,
     1                                    findices, fvalues)

      integer ierr
      integer*8 fvector
      integer fnum_values
      integer findices
      double precision fvalues

      call HYPRE_IJVectorSetValues(fvector, fnum_values, findices,
     1                             fvalues, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorsetvalues error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorAddToValues
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectoraddtovalues(fvector, fnum_values,
     1                                      findices, fvalues)

      integer ierr
      integer*8 fvector
      integer fnum_values
      integer findices
      double precision fvalues

      call HYPRE_IJVectorAddToValues(fvector, fnum_values, findices,
     1                               fvalues, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectoraddtovalues error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorAssemble
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorassemble(fvector)

      integer ierr
      integer*8 fvector

      call HYPRE_IJVectorAssemble(fvector , ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorassemble error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorGetValues
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorgetvalues(fvector, fnum_values,
     1                                    findices, fvalues)

      integer ierr
      integer*8 fvector
      integer fnum_values
      integer findices
      double precision fvalues

      call HYPRE_IJVectorGetValues(fvector, fnum_values, findices,
     1                             fvalues, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorgetvalues error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorSetObjectType
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorsetobjecttype(fvector, ftype)

      integer ierr
      integer*8 fvector
      integer ftype

      call HYPRE_IJVectorSetObjectType(fvector, ftype, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorsetobjecttype error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorGetObjectType
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorgetobjecttype(fvector, ftype)

      integer ierr
      integer*8 fvector
      integer ftype

      call HYPRE_IJVectorGetObjectType(fvector, ftype, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorgetobjecttype error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorGetObject
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorgetobject(fvector, fobject)

      integer ierr
      integer*8 fvector
      integer*8 fobject

      call HYPRE_IJVectorGetObject(fvector, fobject, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorgetobject error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorRead
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorread(ffilename, fcomm, fobject_type,
     1                               fvector)

      integer ierr
      integer*8 fvector
      integer*8 fcomm
      integer fobject_type
      character*(*) ffilename

      call HYPRE_IJVectorRead(ffilename, fcomm, fobject_type, fvector,
     1                        ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorread error = ', ierr
      endif

      return
      end

c--------------------------------------------------------------------------
c HYPRE_IJVectorPrint
c--------------------------------------------------------------------------
      subroutine fhypre_ijvectorprint(fvector, ffilename)

      integer ierr
      integer*8 fvector
      character*(*) ffilename

      call HYPRE_IJVectorPrint(fvector, ffilename, ierr)

      if(ierr .ne. 0) then
         print *, 'fhypre_ijvectorprint error = ', ierr
      endif

      return
      end

max_issues_repo_name: leyel/BoSSS

max_forks_repo_name: leyel/BoSSS


text:
module riemann_mod

use constants_mod

enum, bind(c)

    !! "NUNK" ... FLASH speak
    enumerator :: B5_NUNK_VARS_BEGIN = 0

    enumerator :: DENS_VAR
    enumerator :: MOMX_VAR
    enumerator :: MOMY_VAR
    enumerator :: MOMZ_VAR
    enumerator :: ENER_VAR

    enumerator :: MAGX_VAR
    enumerator :: MAGY_VAR
    enumerator :: MAGZ_VAR

    enumerator :: VELX_VAR
    enumerator :: VELY_VAR
    enumerator :: VELZ_VAR

    enumerator :: PRES_VAR
    enumerator :: EINT_VAR

    enumerator :: GAMC_VAR

    enumerator :: B5_NUNK_VARS_END

endenum

integer, parameter :: B5_NUNK_VARS = B5_NUNK_VARS_END - 1

enum, bind(c)

    enumerator :: B5_FLUX_VARS_BEGIN = 0

    enumerator :: DENS_FLUX
    enumerator :: XMOM_FLUX
    enumerator :: YMOM_FLUX
    enumerator :: ZMOM_FLUX
    enumerator :: ENER_FLUX

    enumerator :: MAGX_FLUX
    enumerator :: MAGY_FLUX
    enumerator :: MAGZ_FLUX

    enumerator :: VELX_FLUX
    enumerator :: VELY_FLUX
    enumerator :: VELZ_FLUX

    enumerator :: PRES_FLUX
    enumerator :: EINT_FLUX

    enumerator :: GAMC_FLUX

    enumerator :: B5_FLUX_VARS_END

endenum

integer, parameter :: B5_FLUX_VARS = B5_FLUX_VARS_END - 1

contains

!! 5-wave Riemann solver (Bouchut, 2008)
pure function xriemann(U_L,U_R) result(F)

    use equations_mod, only: pressure
    use setup_config_mod, only: kappa

    real(dp), intent(in) :: U_L(N_VARS),U_R(N_VARS)
    real(dp)             :: F(N_VARS)

    real(dp) :: speed, vint

    real(dp) :: lm,lp,lam,lap,lbm,lbp

    real(dp) :: cal, car, cbl, cbr, cbovrhol, cbovrhor, caovrhol, caovrhor, aql, aqr, aql2, aqr2
    real(dp) :: pixl, pixr, Bx2l, Bx2r, asl2, asr2, asl, asr, BxBtl,BxBtr, Xl, Xr
    real(dp) :: xsmr, xsml, B2l, B2r, alpha, Etotl, Etotr, ptotl, ptotr

    real(dp), DIMENSION(N_DIMS) :: pil, pir, ustar,  Bstarl, Bstarr, pinstar, ctabl, ctabr, ustarmid
    real(dp) :: pi2l, pi2r,pistar2, rhostarl, rhostarr, eintstarl, eintstarr
    real(dp) :: B2starl, B2starr, u2star, pistar,Etotstarl, Etotstarr
    real(dp) :: dBn

    real(dp) :: Um(B5_NUNK_VARS), Up(B5_NUNK_VARS), Flux(B5_FLUX_VARS)

    integer, parameter :: ivn  = VELX_VAR
    integer, parameter :: ibn  = MAGX_VAR
    integer, parameter :: idir = 1

    Up(DENS_VAR:ENER_VAR) = U_L
    Um(DENS_VAR:ENER_VAR) = U_R

    Up(MAGX_VAR:MAGZ_VAR) = 0.0_dp
    Um(MAGX_VAR:MAGZ_VAR) = 0.0_dp

    Up(VELX_VAR:VELZ_VAR) = Up(MOMX_VAR:MOMZ_VAR)/Up(DENS_VAR)
    Um(VELX_VAR:VELZ_VAR) = Um(MOMX_VAR:MOMZ_VAR)/Um(DENS_VAR)

    Up(PRES_VAR) = pressure(Up(DENS_VAR:ENER_VAR))
    Um(PRES_VAR) = pressure(Um(DENS_VAR:ENER_VAR))

    Up(EINT_VAR) = Up(PRES_VAR)/(kappa-1)/Up(DENS_VAR)
    Um(EINT_VAR) = Um(PRES_VAR)/(kappa-1)/Um(DENS_VAR)

    Up(GAMC_VAR) = kappa
    Um(GAMC_VAR) = kappa

    Flux  = 0.0_dp
    speed = 0.0_dp
    vint  = 0.0_dp

    alpha = 0.5_dp*(max(Up(GAMC_VAR),Um(GAMC_VAR))+1.0_dp) + 1e-4_dp

    Bx2l = Up(ibn)**2
    Bx2r = Um(ibn)**2

    !! Total pressure and energy
    B2l   = Up(MAGX_VAR)**2 + Up(MAGY_VAR)**2 + Up(MAGZ_VAR)**2
    ptotl = Up(PRES_VAR) + 0.5_dp*B2l
    pixl  = ptotl - Bx2l
    Etotl = Up(DENS_VAR)*(Up(EINT_VAR) + 0.5_dp*(Up(VELX_VAR)**2 + Up(VELY_VAR)**2 + Up(VELZ_VAR)**2)) + 0.5_dp*B2l

    B2r   = Um(MAGX_VAR)**2 + Um(MAGY_VAR)**2 + Um(MAGZ_VAR)**2
    ptotr = Um(PRES_VAR) + 0.5_dp*B2r
    pixr  = ptotr - Bx2r
    Etotr = Um(DENS_VAR)*(Um(EINT_VAR) + 0.5_dp*(Um(VELX_VAR)**2 + Um(VELY_VAR)**2 + Um(VELZ_VAR)**2)) + 0.5_dp*B2r

    pil         = 0.0_dp
    pil(idir)   = ptotl
    pil         = pil - Up(ibn)*Up(MAGX_VAR:MAGZ_VAR)
    pi2l        = dot_product(pil,pil)

    pir         = 0.0_dp
    pir(idir)   = ptotr
    pir         = pir - Um(ibn)*Um(MAGX_VAR:MAGZ_VAR)
    pi2r        = dot_product(pir,pir)

    !! Sound soeed
    asl2        = Up(GAMC_VAR)*Up(PRES_VAR)/Up(DENS_VAR)
    asr2        = Um(GAMC_VAR)*Um(PRES_VAR)/Um(DENS_VAR)
    asl2        = max(0.0_dp,asl2)
    asr2        = max(0.0_dp,asr2)
    asl         = sqrt(asl2)
    asr         = sqrt(asr2)

    BxBtl       = abs(Up(ibn))*sqrt(max(0.0_dp,B2l-Bx2l))
    BxBtr       = abs(Um(ibn))*sqrt(max(0.0_dp,B2r-Bx2r))
    aql2        = abs(asl2 + (B2l + BxBtl)/Up(DENS_VAR))
    aqr2        = abs(asr2 + (B2r + BxBtr)/Um(DENS_VAR))
    aql         = sqrt(aql2)
    aqr         = sqrt(aqr2)

    Xl          = max(Up(ivn)-Um(ivn),0.0_dp) + max(pixr-pixl,0.0_dp)/(Up(DENS_VAR)*aql + Um(DENS_VAR)*aqr)
    Xr          = max(Up(ivn)-Um(ivn),0.0_dp) + max(pixl-pixr,0.0_dp)/(Um(DENS_VAR)*aqr + Up(DENS_VAR)*aql)

    xsml        = 1.0_dp - Xl/(aql+alpha*Xl)
    xsmr        = 1.0_dp - Xr/(aqr+alpha*Xr)

    caovrhol    = sqrt((Bx2l + BxBtl)/Up(DENS_VAR)/xsml)
    caovrhor    = sqrt((Bx2r + BxBtr)/Um(DENS_VAR)/xsmr)
    cbovrhol    = aql/sqrt(xsml) + alpha*Xl
    cbovrhor    = aqr/sqrt(xsmr) + alpha*Xr

    !! Avoid ill-posedness.
    caovrhol    = max(caovrhol, 1.0e-12_dp*cbovrhol)
    caovrhor    = max(caovrhor, 1.0e-12_dp*cbovrhor)

    cal         = Up(DENS_VAR)*caovrhol
    car         = Um(DENS_VAR)*caovrhor
    cbl         = Up(DENS_VAR)*cbovrhol
    cbr         = Um(DENS_VAR)*cbovrhor

    !! Velocities
    ustarmid        = (cal*Up(VELX_VAR:VELZ_VAR) + car*Um(VELX_VAR:VELZ_VAR) + pil-pir) / (cal+car)
    ustar(idir)     = (cbl*Up(ivn) + cbr*Um(ivn) + pixl - pixr) / (cbl+cbr)
    ustarmid(idir)  = ustar(idir)

    lbm         = Up(ivn) - cbovrhol
    lbp         = Um(ivn) + cbovrhor

    lam         = Up(ivn) - cal/Up(DENS_VAR)
    lap         = Um(ivn) + car/Um(DENS_VAR)

    lm          = min(lam,lbm)
    lp          = max(lap,lbp)

    !! "Output"
    speed       = max(lp,-lm)
    vint        = ustar(idir)

    !! Computation of states and fluxes:

    if (lm .ge. 0.0_dp) then

        Flux(DENS_FLUX)           = Up(DENS_VAR)*Up(ivn)
        Flux(XMOM_FLUX:ZMOM_FLUX) = Flux(DENS_FLUX)*Up(VELX_VAR:VELZ_VAR) + pil
        Flux(MAGX_FLUX:MAGZ_FLUX) = Up(ivn) * Up(MAGX_VAR:MAGZ_VAR) - Up(ibn)*Up(VELX_VAR:VELZ_VAR)
        Flux(ENER_FLUX)           = Etotl*Up(ivn) + dot_product(pil,Up(VELX_VAR:VELZ_VAR))

    else if (ustar(idir) .gt. 0.0_dp) then

        if (lbm .gt. 0.0_dp) then

            ustar           = (cal*Up(VELX_VAR:VELZ_VAR) + car*Um(VELX_VAR:VELZ_VAR) +pil - pir) / (cal+car)
            ustar(idir)     = Up(ivn)
            pinstar         = (car*pil + cal*pir - cal*car*(Um(VELX_VAR:VELZ_VAR) - Up(VELX_VAR:VELZ_VAR))) / (cal+car)
            pinstar(idir)   = pixl

            rhostarl        = Up(DENS_VAR)

            Bstarl          = Up(MAGX_VAR:MAGZ_VAR) + rhostarl*Up(ibn)*(pil-pinstar)/(cal**2)
            Bstarl(idir)    = Up(ibn)
            u2star          = dot_product(ustar,ustar)
            B2starl         = dot_product(Bstarl,Bstarl)
            pistar2         = dot_product(pinstar,pinstar)
            eintstarl       = Up(EINT_VAR) + 0.5_dp*(B2l/Up(DENS_VAR) - B2starl/rhostarl) - (pi2l - pistar2)/(2.0_dp*cal**2)
            Etotstarl       = rhostarl*(0.5_dp*u2star + eintstarl) + 0.5_dp*B2starl

            Flux(DENS_FLUX)             = rhostarl*ustar(idir)
            Flux(XMOM_FLUX:ZMOM_FLUX)   = rhostarl*ustar(idir)*ustar + pinstar
            Flux(MAGX_FLUX:MAGZ_FLUX)   = ustar(idir)*Bstarl - Up(ibn)*ustar
            Flux(ENER_FLUX)             = Etotstarl*ustar(idir) + dot_product(pinstar(1:3),ustar(1:3))

        else if (lam .lt. 0.0_dp) then

            ctabl       = (/cal, cal, cal/)
            ctabl(idir) = cbl
            ctabr       = (/car, car, car/)
            ctabr(idir) = cbr

            ustar           = (ctabl*Up(VELX_VAR:VELZ_VAR) + ctabr*Um(VELX_VAR:VELZ_VAR) + pil - pir) / (ctabl+ctabr)
            !ustar(idir)    = (cbl*Up(ivn,i-1)+cbr*Um(ivn,i)+pixl-pixr) / (cbl+cbr)
            pinstar         = (ctabr*pil+ctabl*pir-ctabl*ctabr*(Um(VELX_VAR:VELZ_VAR) - Up(VELX_VAR:VELZ_VAR))) / (ctabl+ctabr)
            !pinstar(idir)  = (cbr*pixl+cbl*pixr-cbl*cbr*(Um(ivn,i)-Up(ivn,i-1))) / (cnl+cnr)

            rhostarl        = 1.0_dp/(1.0_dp/Up(DENS_VAR) + (pil(idir) - pinstar(idir))/(cbl**2))
            Bstarl          = rhostarl*(Up(MAGX_VAR:MAGZ_VAR)/Up(DENS_VAR) + Up(ibn)*(pil-pinstar)/(ctabl**2))
            u2star          = dot_product(ustar,ustar)
            B2starl         = dot_product(Bstarl,Bstarl)
            pi2l            = dot_product(pil/ctabl,pil/ctabl)
            pistar2         = dot_product(pinstar/ctabl,pinstar/ctabl)
            eintstarl       = Up(EINT_VAR) + 0.5_dp*(B2l/Up(DENS_VAR) - B2starl/rhostarl) - 0.5_dp*(pi2l-pistar2)
            Etotstarl       = rhostarl*(0.5_dp*u2star + eintstarl) + 0.5_dp*B2starl

            Flux(DENS_FLUX)             = rhostarl*ustar(idir)
            Flux(XMOM_FLUX:ZMOM_FLUX)   = rhostarl*ustar(idir)*ustar + pinstar
            Flux(MAGX_FLUX:MAGZ_FLUX)   = ustar(idir)*Bstarl - Up(ibn)*ustar
            Flux(ENER_FLUX)             = Etotstarl*ustar(idir)+dot_product(pinstar(1:3),ustar(1:3))

        else !lam > 0, lbm < 0

            ustar           = Up(VELX_VAR:VELZ_VAR)
            ustar(idir)     = vint
            pinstar         = pil
            pinstar(idir)   = (cbr*pixl+cbl*pixr-cbl*cbr*(Um(ivn)-Up(ivn))) / (cbl+cbr)

            rhostarl        = 1.0_dp/(1.0_dp/Up(DENS_VAR) + (pil(idir)-pinstar(idir))/(cbl**2))
            Bstarl          = rhostarl*(Up(MAGX_VAR:MAGZ_VAR)/Up(DENS_VAR))
            Bstarl(idir)    = Up(ibn)
            u2star          = dot_product(ustar,ustar)
            B2starl         = dot_product(Bstarl,Bstarl)
            !pi2l           = dot_product(pil,pil)
            !pistar2        = dot_product(pinstar,pinstar)
            eintstarl       = Up(EINT_VAR) + 0.5_dp*(B2l/Up(DENS_VAR) - B2starl/rhostarl)-(pixl**2-pinstar(idir)**2)/(2.0_dp*cbl**2)
            Etotstarl       = rhostarl*(u2star/2. + eintstarl) + B2starl/2

            Flux(DENS_FLUX)             = rhostarl*ustar(idir)
            Flux(XMOM_FLUX:ZMOM_FLUX)   = rhostarl*ustar(idir)*ustar + pinstar
            Flux(MAGX_FLUX:MAGZ_FLUX)   = ustar(idir)*Bstarl - Up(ibn)*ustar
            Flux(ENER_FLUX)             = Etotstarl*ustar(idir) + dot_product(pinstar(1:3),ustar(1:3))

        endif ! Done with left wave fan

    else if (lp.gt.0.0_dp) then

        if (lbp .lt. 0.0_dp ) then !lap>0

            ustar           = (cal*Up(VELX_VAR:VELZ_VAR) + car*Um(VELX_VAR:VELZ_VAR) + pil-pir) / (cal+car)
            ustar(idir)     = Um(ivn)
            pinstar         = (car*pil + cal*pir - cal*car*(Um(VELX_VAR:VELZ_VAR) - Up(VELX_VAR:VELZ_VAR))) / (cal+car)
            pinstar(idir)   = pixr

            rhostarr        = Um(DENS_VAR)

            Bstarr          = Um(MAGX_VAR:MAGZ_VAR) + rhostarr*Um(ibn)*(pir-pinstar)/(car**2)
            Bstarr(idir)    = Um(ibn)
            u2star          = dot_product(ustar,ustar)
            B2starr         = dot_product(Bstarr,Bstarr)
            pistar2         = dot_product(pinstar,pinstar)
            eintstarr       = Um(EINT_VAR) + (B2r/Um(DENS_VAR) - B2starr/rhostarr)/2 - (pi2r-pistar2)/(2*car**2)
            Etotstarr       = rhostarr*(u2star/2 + eintstarr) + B2starr/2

            Flux(DENS_FLUX)             = rhostarr*ustar(idir)
            Flux(XMOM_FLUX:ZMOM_FLUX)   = rhostarr*ustar(idir)*ustar + pinstar
            Flux(MAGX_FLUX:MAGZ_FLUX)   = ustar(idir)*Bstarr - Um(ibn)*ustar
            Flux(ENER_FLUX)             = Etotstarr*ustar(idir) + dot_product(pinstar(1:3),ustar(1:3))

        else if (lap .gt. 0) then !lbp>0, lap>0

            ctabl           = (/cal, cal, cal/)
            ctabl(idir)     = cbl
            ctabr           = (/car, car,car/)
            ctabr(idir)     = cbr

            ustar           = (ctabl*Up(VELX_VAR:VELZ_VAR) + ctabr*Um(VELX_VAR:VELZ_VAR) + pil-pir) / (ctabl+ctabr)
            !ustar(idir)    = (cbl*Up(ivn,i-1)+cbr*Um(ivn,i)+pixl-pixr) / (cbl+cbr)
            pinstar         = (ctabr*pil+ctabl*pir-ctabl*ctabr*(Um(VELX_VAR:VELZ_VAR) - Up(VELX_VAR:VELZ_VAR))) / (ctabl+ctabr)
            !pinstar(idir)  = (cbr*pixl+cbl*pixr-cbl*cbr*(Um(ivn,i)-Up(ivn,i-1))) / (cnl+cnr)

            rhostarr        = 1.0_dp/(1.0_dp/Um(DENS_VAR) + (pir(idir) - pinstar(idir))/(cbr**2))
            Bstarr          = rhostarr*(Um(MAGX_VAR:MAGZ_VAR)/Um(DENS_VAR) + Um(ibn)*(pir-pinstar)/(ctabr**2))
            u2star          = dot_product(ustar,ustar)
            B2starr         = dot_product(Bstarr,Bstarr)
            pi2r            = dot_product(pir/ctabr,pir/ctabr)
            pistar2         = dot_product(pinstar/ctabr,pinstar/ctabr)
            eintstarr       = Um(EINT_VAR) + (B2r/Um(DENS_VAR) - B2starr/rhostarr)/2 - (pi2r-pistar2)/2
            Etotstarr       = rhostarr*(u2star/2 + eintstarr) + B2starr/2

            Flux(DENS_FLUX)             = rhostarr*ustar(idir)
            Flux(XMOM_FLUX:ZMOM_FLUX)   = rhostarr*ustar(idir)*ustar + pinstar
            Flux(MAGX_FLUX:MAGZ_FLUX)   = ustar(idir)*Bstarr - Um(ibn)*ustar
            Flux(ENER_FLUX)             = Etotstarr*ustar(idir) + dot_product(pinstar(1:3),ustar(1:3))

        else !lbp>0, lap<=0

            ustar           = Um(VELX_VAR:VELZ_VAR)
            ustar(idir)     = vint
            pinstar         = pir
            pinstar(idir)   = (cbr*pixl + cbl*pixr - cbl*cbr*(Um(ivn)-Up(ivn))) / (cbl+cbr)

            rhostarr        = 1.0_dp/(1.0_dp/Um(DENS_VAR) + (pir(idir) - pinstar(idir))/(cbr**2))
            Bstarr          = rhostarr*(Um(MAGX_VAR:MAGZ_VAR)/Um(DENS_VAR))
            Bstarr(idir)    = Um(ibn)
            u2star          = dot_product(ustar,ustar)
            B2starr         = dot_product(Bstarr,Bstarr)
            !pi2l=dot_product(pil,pil)
            !pistar2=dot_product(pinstar,pinstar)
            eintstarr       = Um(EINT_VAR) + (B2r/Um(DENS_VAR) - B2starr/rhostarr)/2 - (pixr**2-pinstar(idir)**2)/(2*cbr**2)
            Etotstarr       = rhostarr*(u2star/2 + eintstarr) + B2starr/2

            Flux(DENS_FLUX)             = rhostarr*ustar(idir)
            Flux(XMOM_FLUX:ZMOM_FLUX)   = rhostarr*ustar(idir)*ustar + pinstar
            Flux(MAGX_FLUX:MAGZ_FLUX)   = ustar(idir)*Bstarr - Um(ibn)*ustar
            Flux(ENER_FLUX)             = Etotstarr*ustar(idir) + dot_product(pinstar(1:3),ustar(1:3))

        endif !! End of right inner states

    else ! Leftwards supermagnetosonic

        Flux(DENS_FLUX)             = Um(DENS_VAR)*Um(ivn)
        Flux(XMOM_FLUX:ZMOM_FLUX)   = Flux(DENS_FLUX)*Um(VELX_VAR:VELZ_VAR) + pir
        Flux(MAGX_FLUX:MAGZ_FLUX)   = Um(ivn)*Um(MAGX_VAR:MAGZ_VAR) - Um(ibn)*Um(VELX_VAR:VELZ_VAR)
        Flux(ENER_FLUX)             = Etotr*Um(ivn) + dot_product(pir,Um(VELX_VAR:VELZ_VAR))

    endif

    !! Nonconservative bit added to source
    !! dBn = Um(ibn)-Up(ibn)
    !! if (vint .le. 0.0_dp) then
    !!     Src(MAGX_VAR:MAGZ_VAR) = Src(MAGX_VAR:MAGZ_VAR) - ustarmid*dBn/dx
    !! else
    !!     Src(MAGX_VAR:MAGZ_VAR) = Src(MAGX_VAR:MAGZ_VAR) - ustarmid*dBn/dx
    !! endif

    F(DENS_VAR)             = Flux(DENS_FLUX)
    F(MOMX_VAR:MOMZ_VAR)    = Flux(XMOM_FLUX:ZMOM_FLUX)
    F(ENER_VAR)             = Flux(ENER_FLUX)

end function

pure function yriemann(uL,uR) result(sf)

    use equations_const_mod

    real(dp), intent(in) :: uL(N_VARS), uR(N_VARS)
    real(dp)             :: sf(N_VARS)

    real(dp) :: ruL(N_VARS), ruR(N_VARS), rf(N_VARS) !! rotated

    ruL(DENS_VAR) = uL(DENS_VAR)
    ruL(MOMX_VAR) = uL(MOMY_VAR)
    ruL(MOMY_VAR) = uL(MOMZ_VAR)
    ruL(MOMZ_VAR) = uL(MOMX_VAR)
    ruL(ENER_VAR) = uL(ENER_VAR)

    ruR(DENS_VAR) = uR(DENS_VAR)
    ruR(MOMX_VAR) = uR(MOMY_VAR)
    ruR(MOMY_VAR) = uR(MOMZ_VAR)
    ruR(MOMZ_VAR) = uR(MOMX_VAR)
    ruR(ENER_VAR) = uR(ENER_VAR)

    rf = xriemann(ruL,ruR)

    sf(DENS_VAR) = rf(DENS_VAR)
    sf(MOMX_VAR) = rf(MOMZ_VAR)
    sf(MOMY_VAR) = rf(MOMX_VAR)
    sf(MOMZ_VAR) = rf(MOMY_VAR)
    sf(ENER_VAR) = rf(ENER_VAR)

end function

pure function zriemann(uL,uR) result(sf)

    use equations_const_mod

    real(dp), intent(in) :: uL(N_VARS), uR(N_VARS)
    real(dp)             :: sf(N_VARS)

    real(dp) :: ruL(N_VARS), ruR(N_VARS), rf(N_VARS) !! rotated

    ruL(DENS_VAR) = uL(DENS_VAR)
    ruL(MOMX_VAR) = uL(MOMZ_VAR)
    ruL(MOMY_VAR) = uL(MOMX_VAR)
    ruL(MOMZ_VAR) = uL(MOMY_VAR)
    ruL(ENER_VAR) = uL(ENER_VAR)

    ruR(DENS_VAR) = uR(DENS_VAR)
    ruR(MOMX_VAR) = uR(MOMZ_VAR)
    ruR(MOMY_VAR) = uR(MOMX_VAR)
    ruR(MOMZ_VAR) = uR(MOMY_VAR)
    ruR(ENER_VAR) = uR(ENER_VAR)

    rf = xriemann(ruL,ruR)

    sf(DENS_VAR) = rf(DENS_VAR)
    sf(MOMX_VAR) = rf(MOMY_VAR)
    sf(MOMY_VAR) = rf(MOMZ_VAR)
    sf(MOMZ_VAR) = rf(MOMX_VAR)
    sf(ENER_VAR) = rf(ENER_VAR)

end function

end module

max_stars_repo_name: jmark/nemo

max_issues_repo_name: jmark/nemo

max_forks_repo_name: jmark/nemo


MODULE Phyllochronmod
    IMPLICIT NONE
CONTAINS

    SUBROUTINE model_Phyllochron(fixPhyll, &
        leafNumber, &
        lincr, &
        ldecr, &
        pdecr, &
        pincr, &
        ptq, &
        gAImean, &
        kl, &
        pTQhf, &
        B, &
        p, &
        choosePhyllUse, &
        areaSL, &
        areaSS, &
        lARmin, &
        lARmax, &
        sowingDensity, &
        lNeff, &
        phyllochron)
        IMPLICIT NONE
        REAL, INTENT(IN) :: fixPhyll
        REAL, INTENT(IN) :: leafNumber
        REAL, INTENT(IN) :: lincr
        REAL, INTENT(IN) :: ldecr
        REAL, INTENT(IN) :: pdecr
        REAL, INTENT(IN) :: pincr
        REAL, INTENT(IN) :: ptq
        REAL, INTENT(IN) :: gAImean
        REAL, INTENT(IN) :: kl
        REAL, INTENT(IN) :: pTQhf
        REAL, INTENT(IN) :: B
        REAL, INTENT(IN) :: p
        CHARACTER(65), INTENT(IN) :: choosePhyllUse
        REAL, INTENT(IN) :: areaSL
        REAL, INTENT(IN) :: areaSS
        REAL, INTENT(IN) :: lARmin
        REAL, INTENT(IN) :: lARmax
        REAL, INTENT(IN) :: sowingDensity
        REAL, INTENT(IN) :: lNeff
        REAL, INTENT(OUT) :: phyllochron
        REAL:: gaiLim
        REAL:: LAR
        !- Name: Phyllochron -Version: 1.0, -Time step: 1
        !- Description:
    !            * Title: Phyllochron Model
    !            * Author: Pierre Martre
    !            * Reference: Modeling development phase in the
    !                Wheat Simulation Model SiriusQuality.
    !                See documentation at http://www1.clermont.inra.fr/siriusquality/?page_id=427
    !            * Institution: INRA Montpellier
    !            * Abstract: Calculate different types of phyllochron
        !- inputs:
    !            * name: fixPhyll
    !                          ** description : Sowing date corrected Phyllochron
    !                          ** inputtype : variable
    !                          ** variablecategory : auxiliary
    !                          ** datatype : DOUBLE
    !                          ** default : 5.0
    !                          ** min : 0.0
    !                          ** max : 10000.0
    !                          ** unit : °C d leaf-1
    !                          ** uri : some url
    !            * name: leafNumber
    !                          ** description : Actual number of phytomers
    !                          ** inputtype : variable
    !                          ** variablecategory : state
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 25.0
    !                          ** unit : leaf
    !                          ** uri : some url
    !            * name: lincr
    !                          ** description : Leaf number above which the phyllochron is increased by Pincr
    !                          ** inputtype : parameter
    !                          ** parametercategory : species
    !                          ** datatype : DOUBLE
    !                          ** default : 8.0
    !                          ** min : 0.0
    !                          ** max : 30.0
    !                          ** unit : leaf
    !                          ** uri : some url
    !            * name: ldecr
    !                          ** description : Leaf number up to which the phyllochron is decreased by Pdecr
    !                          ** inputtype : parameter
    !                          ** parametercategory : species
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 100.0
    !                          ** unit : leaf
    !                          ** uri : some url
    !            * name: pdecr
    !                          ** description : Factor decreasing the phyllochron for leaf number less than Ldecr
    !                          ** inputtype : parameter
    !                          ** parametercategory : species
    !                          ** datatype : DOUBLE
    !                          ** default : 0.4
    !                          ** min : 0.0
    !                          ** max : 10.0
    !                          ** unit : -
    !                          ** uri : some url
    !            * name: pincr
    !                          ** description : Factor increasing the phyllochron for leaf number higher than Lincr
    !                          ** inputtype : parameter
    !                          ** parametercategory : species
    !                          ** datatype : DOUBLE
    !                          ** default : 1.5
    !                          ** min : 0.0
    !                          ** max : 10.0
    !                          ** unit : -
    !                          ** uri : some url
    !            * name: ptq
    !                          ** description : Photothermal quotient
    !                          ** inputtype : variable
    !                          ** variablecategory : state
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 10000.0
    !                          ** unit : MJ °C-1 d-1 m-2)
    !                          ** uri : some url
    !            * name: gAImean
    !                          ** description : Green Area Index
    !                          ** inputtype : variable
    !                          ** variablecategory : state
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 10000.0
    !                          ** unit : m2 m-2
    !                          ** uri : some url
    !            * name: kl
    !                          ** description : Exctinction Coefficient
    !                          ** inputtype : parameter
    !                          ** parametercategory : species
    !                          ** datatype : DOUBLE
    !                          ** default : 0.45
    !                          ** min : 0.0
    !                          ** max : 50.0
    !                          ** unit : -
    !                          ** uri : some url
    !            * name: pTQhf
    !                          ** description : Slope to intercept ratio for Phyllochron  parametrization with PhotoThermal Quotient
    !                          ** inputtype : parameter
    !                          ** parametercategory : genotypic
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : °C d leaf-1
    !                          ** uri : some url
    !            * name: B
    !                          ** description : Phyllochron at PTQ equal 1
    !                          ** inputtype : parameter
    !                          ** parametercategory : species
    !                          ** datatype : DOUBLE
    !                          ** default : 20.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : °C d leaf-1
    !                          ** uri : some url
    !            * name: p
    !                          ** description : Phyllochron (Varietal parameter)
    !                          ** inputtype : parameter
    !                          ** parametercategory : species
    !                          ** datatype : DOUBLE
    !                          ** default : 120.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : °C d leaf-1
    !                          ** uri : some url
    !            * name: choosePhyllUse
    !                          ** description : Switch to choose the type of phyllochron calculation to be used
    !                          ** inputtype : parameter
    !                          ** parametercategory : species
    !                          ** datatype : STRING
    !                          ** default : Default
    !                          ** min :
    !                          ** max :
    !                          ** unit : -
    !                          ** uri : some url
    !            * name: areaSL
    !                          ** description :  Area Leaf
    !                          ** inputtype : parameter
    !                          ** parametercategory : genotypic
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : cm2
    !                          ** uri : some url
    !            * name: areaSS
    !                          ** description : Area Sheath
    !                          ** inputtype : parameter
    !                          ** parametercategory : genotypic
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : cm2
    !                          ** uri : some url
    !            * name: lARmin
    !                          ** description : LAR minimum
    !                          ** inputtype : parameter
    !                          ** parametercategory : genotypic
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : leaf-1 °C
    !                          ** uri : some url
    !            * name: lARmax
    !                          ** description : LAR maximum
    !                          ** inputtype : parameter
    !                          ** parametercategory : genotypic
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : leaf-1 °C
    !                          ** uri : some url
    !            * name: sowingDensity
    !                          ** description : Sowing Density
    !                          ** inputtype : parameter
    !                          ** parametercategory : genotypic
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : plant m-2
    !                          ** uri : some url
    !            * name: lNeff
    !                          ** description : Leaf Number efficace
    !                          ** inputtype : parameter
    !                          ** parametercategory : genotypic
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 1000.0
    !                          ** unit : leaf
    !                          ** uri : some url
        !- outputs:
    !            * name: phyllochron
    !                          ** description :  the rate of leaf appearance
    !                          ** variablecategory : state
    !                          ** datatype : DOUBLE
    !                          ** min : 0
    !                          ** max : 1000
    !                          ** unit :  °C d leaf-1
    !                          ** uri : some url
        phyllochron = 0.0
        LAR = 0.0
        gaiLim = lNeff * ((areaSL + areaSS) / 10000.0) * sowingDensity
        IF(choosePhyllUse .EQ. 'Default') THEN
            IF(leafNumber .LT. ldecr) THEN
                phyllochron = fixPhyll * pdecr
            ELSE IF ( leafNumber .GE. ldecr .AND. leafNumber .LT. lincr) THEN
                phyllochron = fixPhyll
            ELSE
                phyllochron = fixPhyll * pincr
            END IF
        END IF
        IF(choosePhyllUse .EQ. 'PTQ') THEN
            IF(gAImean .GT. gaiLim) THEN
                LAR = (lARmin + ((lARmax - lARmin) * ptq / (pTQhf + ptq))) / (B *  &
                        gAImean)
            ELSE
                LAR = (lARmin + ((lARmax - lARmin) * ptq / (pTQhf + ptq))) / (B *  &
                        gaiLim)
            END IF
            phyllochron = 1.0 / LAR
        END IF
        IF(choosePhyllUse .EQ. 'Test') THEN
            IF(leafNumber .LT. ldecr) THEN
                phyllochron = p * pdecr
            ELSE IF ( leafNumber .GE. ldecr .AND. leafNumber .LT. lincr) THEN
                phyllochron = p
            ELSE
                phyllochron = p * pincr
            END IF
        END IF
    END SUBROUTINE model_Phyllochron

END MODULE


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

C=======================================================================
      SUBROUTINE SCAELB (A, USESEL, IELBST,
     &   VALMN, NUMMN, XYZMN, ISTMN, VALMX, NUMMX, XYZMX, ISTMX,
     &   VALMIN, NUMMIN, XYZMIN, ISTMIN, VALMAX, NUMMAX, XYZMAX, ISTMAX)
C=======================================================================

C   --*** SCAELB *** (BLOT) Scale element variable for selected blocks
C   --   Written by Amy Gilkey - revised 04/01/88
C   --
C   --SCAELB returns the minimum and maximum values for the selected
C   --element blocks.  The minimum and maximums for each block have already
C   --been calculated and stored.
C   --
C   --Parameters:
C   --   A - IN - the dynamic memory base array
C   --   USESEL - IN - use the element blocks selected array iff true,
C   --      else all selected
C   --   IELBST - IN - the element block status (>0 if selected)
C   --   VALMN, VALMX - IN - the minimum and maximum value for each
C   --      element block
C   --   NUMMN, NUMMX - IN - the element number of the minimum and
C   --      maximum value for each element block
C   --   XYZMN, XYZMX - IN - the coordinates of NUMMN, NUMMX
C   --   ISTMN, ISTMX - IN - the step number of the minimum and maximum
C   --      value for each element block
C   --   VALMIN, VALMAX - OUT - the minimum and maximum value
C   --      (with selected element block and birth/death)
C   --   NUMMIN, NUMMAX - OUT - the element number of the minimum and
C   --      maximum value
C   --   XYZMIN, XYZMAX - OUT - the coordinates of NUMMIN, NUMMAX
C   --   ISTMIN, ISTMAX - OUT - the step number of the minimum and maximum value
C   --
C   --Common Variables:
C   --   Uses NUMEL, NELBLK, NVAREL, NSTEPS of /DBNUMS/

      include 'dbnums.blk'

      DIMENSION A(*)
      LOGICAL USESEL
      INTEGER IELBST(*)
      REAL VALMN(0:NELBLK), VALMX(0:NELBLK)
      INTEGER NUMMN(0:NELBLK), NUMMX(0:NELBLK)
      REAL XYZMN(3,0:NELBLK), XYZMX(3,0:NELBLK)
      INTEGER ISTMN(0:NELBLK), ISTMX(0:NELBLK)
      REAL XYZMIN(3), XYZMAX(3)

      LOGICAL INIT

C   --Calculate min/max for element by element block

      IF (USESEL) THEN
         INIT = .TRUE.
         DO 100 IELB = 1, NELBLK
            IF ((.NOT. USESEL) .OR. (IELBST(IELB) .GT. 0)) THEN
               IF (ISTMN(IELB) .GT. 0) THEN
                  IF (INIT) THEN
                     IMIN = IELB
                     IMAX = IELB
                     INIT = .FALSE.
                  ELSE
                     IF (VALMN(IMIN) .GT. VALMN(IELB)) THEN
                        IMIN = IELB
                     END IF
                     IF (VALMX(IMAX) .LT. VALMX(IELB)) THEN
                        IMAX = IELB
                     END IF
                  END IF
               END IF
            END IF
  100    CONTINUE
      ELSE
         INIT = .FALSE.
         IMIN = 0
         IMAX = 0
      END IF

      IF (INIT) THEN
         VALMIN = 0.0
         NUMMIN = 0
         DO 110 I = 1, 3
            XYZMIN(I) = 0.0
  110    CONTINUE
         ISTMIN = 0
         VALMAX = 0.0
         NUMMAX = 0
         DO 120 I = 1, 3
            XYZMAX(I) = 0.0
  120    CONTINUE
         ISTMAX = 0
      ELSE
         VALMIN = VALMN(IMIN)
         NUMMIN = NUMMN(IMIN)
         DO 130 I = 1, 3
            XYZMIN(I) = XYZMN(I,IMIN)
  130    CONTINUE
         ISTMIN = ISTMN(IMIN)
         VALMAX = VALMX(IMAX)
         NUMMAX = NUMMX(IMAX)
         DO 140 I = 1, 3
            XYZMAX(I) = XYZMX(I,IMAX)
  140    CONTINUE
         ISTMAX = ISTMX(IMAX)
      END IF

      RETURN
      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
c { dg-do run }
      PROGRAM LABUG1

*  This program core dumps on mips-sgi-irix6.2 when compiled
*  with egcs-19981101, egcs-19981109 and egcs-19981122 snapshots
*  with -O2
*
*  Originally derived from LAPACK test suite.
*  Almost any change allows it to run.
*
*  David Billinghurst, (David.Billinghurst@riotinto.com.au)
*  25 November 1998
*
*     .. Parameters ..
      INTEGER   LDA, LDE
      PARAMETER ( LDA = 2500, LDE = 50  )
      COMPLEX   CZERO
      PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )

      INTEGER   I, J, M, N
      REAL      V
      COMPLEX   A(LDA),B(LDA),C(LDA),E(LDE,LDE),F(LDE,LDE)
      COMPLEX   Z

      N=2
      M=1
*
      do i = 1, m
         do j = 1, n
            e(i,j) = czero
            f(i,j) = czero
        end do
      end do
*
      DO J = 1, N
         DO I = 1, M
            V =  ABS( E(I,J) - F(I,J) )
         END DO
      END DO

      CALL SUB2(M,Z)

      END

      subroutine SUB2(I,A)
      integer i
      complex a
      end

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
        logical function f1()
        f1 = .TRUE.
        end

        logical*2 function f2()
        f2 = .TRUE.
        end

        character*5 function f3()
        f3 = 'abcde'
        end

        function f4()
        logical f4
        f4 = .TRUE.
        end

        integer*4 function f5()
        integer*4 f5
        f5 = 1234
        end

        character*6 function f6()
        character*6 f6
        f6 = 'abcdef'
        end

        character*(*) function f7()
        f7 = 'AB'
        end

max_stars_repo_name: camfort/fortran-vars

max_issues_repo_name: camfort/fortran-vars

max_forks_repo_name: camfort/fortran-vars