text:
module type_file_properties
#include <messenger.h>
    use mod_kinds,  only: rk, ik
    implicit none

    !>
    !!  @author Nathan A. Wukie
    !!  @date   2/3/2016
    !!
    !------------------------------------------------------------------------------------------------------------
    type, public :: file_properties_t

        logical     :: contains_grid           = .false.
        logical     :: contains_primary_fields = .false.
        logical     :: contains_adjoint_fields = .false.
        logical     :: adjoint_mode            = .false.

        integer(ik)                      :: ndomains = 0      !< Number of domains in the file
        character(len=1024)              :: time_integrator   !< Time integrator name
        character(len=1024), allocatable :: domain_names(:)
        character(len=1024), allocatable :: eqnset(:)         !< Equation set specified in the file.
        integer(ik),         allocatable :: order_c(:)        !< Coordinate order ( 1st, 2nd, 3rd, etc. )
        integer(ik),         allocatable :: order_s(:)        !< Solution order ( 1st, 2nd, 3rd, etc. )
        integer(ik),         allocatable :: nterms_c(:)       !< Number of terms in the coordinate expansion
        integer(ik),         allocatable :: nterms_s(:)       !< Number of terms in the solution expansion
        integer(ik)                      :: nfunctionals      !< Number of functionals (ie number adjoint vectors)
        integer(ik)                      :: istep             !< istep at which the solution was written out

    contains

        procedure   :: set_ndomains

    end type file_properties_t
    !************************************************************************************************************

contains

    !> Initialize the file_properties type by allocating the other components to the right sizes
    !! according to the number of domains detected.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/3/2016
    !!
    !------------------------------------------------------------------------------------------------------------
    subroutine set_ndomains(self,ndomains)
        class(file_properties_t),   intent(inout)   :: self
        integer(ik),                intent(in)      :: ndomains

        integer :: ierr

        ! Set number of domains in the file
        self%ndomains = ndomains

        ! Allocate storage for other data
        if ( allocated(self%order_c)  ) deallocate(self%order_c )
        if ( allocated(self%order_s)  ) deallocate(self%order_s )
        if ( allocated(self%nterms_c) ) deallocate(self%nterms_c)
        if ( allocated(self%nterms_s) ) deallocate(self%nterms_s)
        if ( allocated(self%eqnset)   ) deallocate(self%eqnset  )

        allocate(self%order_c(ndomains),     &
                 self%order_s(ndomains),     &
                 self%nterms_c(ndomains),    &
                 self%nterms_s(ndomains),    &
                 self%eqnset(ndomains), stat=ierr )

        if (ierr /= 0) call AllocationError

    end subroutine set_ndomains
    !***********************************************************************************************************

end module type_file_properties

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
subroutine mapc2p(xc, yc, zc, xp, yp, zp)

    use fault_module, only: center, theta, xcb, ycb

    implicit none

    real(kind=8), intent(in) :: xc, yc, zc
    real(kind=8), intent(out) :: xp, yp, zp

    ! Local variables
    real (kind=8) :: ls, alpha, xrot, zrot

    if (xc < xcb(1)) then
      ls = dsqrt((xc-xcb(1))**2 + (zc-center(3))**2)
    elseif (xc > xcb(2)) then
      ls = dsqrt((xc-xcb(2))**2 + (zc-center(3))**2)
    else
      ls = dabs(zc - center(3))
    end if

    alpha = ls/(-center(3))
    xrot = center(1) + dcos(theta)*(xc-center(1)) + dsin(theta)*(zc-center(3))
    zrot = center(3) - dsin(theta)*(xc-center(1)) + dcos(theta)*(zc-center(3))

    yp = yc
    if (alpha < 1.d0) then
      xp = (1.d0-alpha)*xrot + alpha*xc
      zp = (1.d0-alpha)*zrot + alpha*zc
    else
      xp = xc
      zp = zc
    end if

    return
end

max_stars_repo_name: rjleveque/seismic

max_issues_repo_name: clawpack/seismic

max_forks_repo_name: rjleveque/seismic


text:
c===========================================================================
c
c   This file is part of TISEAN
c
c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
c
c   TISEAN is free software; you can redistribute it and/or modify
c   it under the terms of the GNU General Public License as published by
c   the Free Software Foundation; either version 2 of the License, or
c   (at your option) any later version.
c
c   TISEAN is distributed in the hope that it will be useful,
c   but WITHOUT ANY WARRANTY; without even the implied warranty of
c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c   GNU General Public License for more details.
c
c   You should have received a copy of the GNU General Public License
c   along with TISEAN; if not, write to the Free Software
c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
c
c===========================================================================
c   Wiener filter (2): filter using periodogram in file
c   author T. Schreiber (1998)
c===========================================================================
      parameter(nx=1000000)
      dimension x(nx), a(nx)
      character*72 file, fout, ffin
      data h/1./, dh/0./
      data iverb/1/

      call whatido("Wiener filter (second part)",iverb)
      h=fcan("f",h)
      dh=fcan("w",dh)
      nmax=ican("l",nx)
      nexcl=ican("x",0)
      jcol=ican("c",0)
      isfin=igetfin(ffin,iverb)
      isout=igetout(fout,iverb)
      call nthstring(1,file)
      if(file.eq."-") stop "wiener2: cannot read stdin"
      call readfile(nmax,x,nexcl,jcol,file,iverb)
      call normal(nmax,x,sc,sd)
      nmaxp=nmore(nmax)

      if(dh.eq.0.) dh=h/nmaxp
      ibin=nmaxp*dh/(2*h)
      if(ibin.gt.0) write(istderr(),*)
     .   "wiener1: binning", 2*ibin+1," frequencies"
      if(isout.eq.1) call addsuff(fout,file,"_amp")
      if(fout.eq." ") fout="-"
      call infile(fout,iunit,iverb)
      read(iunit,*) dum, a(1)
      do 10 i=2+ibin,(nmaxp+1)/2-ibin,2*ibin+1
 10      read(iunit,*) dum, a(2*i-2)
      if(mod(nmaxp,2).eq.0)  read(iunit,*)
     .   dum, a(nmaxp)
      d=tospec(nmaxp,a,x,ibin)
      if(iv_io(iverb).eq.1) write(istderr(),*) "rms correction: ", d
      if(isfin.eq.1) call addsuff(ffin,file,"_wc")
      do 20 n=1,nmax
 20      x(n)=x(n)+sc
      call writefile(nmax,x,ffin,iverb)
      end

      subroutine usage()
c usage message

      call whatineed(
     .   "[-f# -w# -o outfile1 -O outfile -l# -x# -c# -V# -h] file")
      call ptext("to provide periodogram, first run:"//
     .   " wiener1 [-f# -w# -o outfile -l# -x# -c# -V# -h] file")
      call ptext("make sure -f# -w# are the same in both wiener calls")
      call popt("f","sampling rate (e.g. in Hz, default 1.)")
      call popt("w","frequency resolution (e.g. in Hz, default 1/N)")
      call popt("l","number of values to be read (all)")
      call popt("x","number of values to be skipped (0)")
      call popt("c","column to be read (1 or file,#)")
      call popt("o","output file of wiener1, just -o means file_amp")
      call popt("O","final output file name, just -O means file_wc")
      call pall()
      call ptext("Note: ""-"" not accepted as file")
      write(istderr(),'()')
      stop
      end

      function igetfin(fout,iverb)
c gets alternate output file name, default " "
c return 1 if fout must be determined from input file name
      character*(*) fout

      igetfin=0
      call stcan("O",fout," ")
      if(fout.ne." ") return
      igetfin=lopt("O",1)
      end

ext: f

max_stars_repo_name: homiak/jabaws

max_issues_repo_name: homiak/jabaws

max_forks_repo_name: homiak/jabaws


text: Phil Berghuis lives on Oeste Drive Oeste near the dorms, he enjoys his punk, his shisha and his veggie dinners.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


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

      SUBROUTINE CULL3 (DIRCOS, MASSLV, NIQM)
      DIMENSION DIRCOS(5,*), MASSLV(2,*)

      J = 0
      DO 10 I=1, NIQM
          IF (MASSLV(2,I) .NE. 0) THEN
              J = J + 1
              MASSLV(1, J) = MASSLV(1, I)
              MASSLV(2, J) = MASSLV(2, I)
              DIRCOS(1, J) = DIRCOS(1, I)
              DIRCOS(2, J) = DIRCOS(2, I)
              DIRCOS(3, J) = DIRCOS(3, I)
              DIRCOS(4, J) = DIRCOS(4, I)
              DIRCOS(5, J) = DIRCOS(5, I)
          END IF
   10 CONTINUE

      NIQM = J
      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 01/01                                 *
C  *    TEST TITLE : Testing the system state value        *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C system state value
      INTEGER    PPHCL,     PPHOP
      PARAMETER (PPHCL = 0, PPHOP = 1)

      INTEGER SYSTRA

      CALL INITGL ('01/01')

C  Use <inquire system state value> to determine:
C      systra = system state value

      CALL SETMSG ('3', '<Inquire system state value> should return ' //
     1             'closed as the system state value before PHIGS '   //
     2             'is opened.')

      SYSTRA = -6
      CALL PQSYS (SYSTRA)
      CALL IFPF (SYSTRA .EQ. PPHCL)

      CALL SETMSG ('1 3', 'After opening PHIGS <inquire system '      //
     1             'state value> should return the system state '     //
     2             'value as open.')

      CALL XPOPPH (ERRFIL, MEMUN)
      CALL PQSYS (SYSTRA)
      CALL IFPF (SYSTRA .EQ. PPHOP)

      CALL SETMSG ('2 3', 'After closing PHIGS <inquire system '      //
     1             'state value> should return the system state '     //
     2             'value as closed.')

      CALL PCLPH
      CALL PQSYS (SYSTRA)
      CALL IFPF (SYSTRA .EQ. PPHCL)

      CALL SETMSG ('1 3', 'After reopening PHIGS <inquire system '    //
     1             'state value> should return the system state '     //
     2             'value as open.')

      CALL XPOPPH (ERRFIL, MEMUN)
      CALL PQSYS (SYSTRA)
      CALL IFPF (SYSTRA .EQ. PPHOP)

666   CONTINUE
C close PHIGS
      CALL PCLPH
      CALL WINDUP
      END


text:
!############################################################################################
!
! Copyright 2021 Bharat Mahajan
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!    http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
!> \brief       Lorenz Test Main Program
!! \details     Main program to run Lorenz test.
!! \author      Bharat Mahajan
!! \date        09/17/2021
!
!############################################################################################

module LorenzDiffEq

    use FLINT

    use, intrinsic :: IEEE_ARITHMETIC

    implicit none

    ! user diff eq system

    type, extends(DiffEqSys) :: LorenzSys
        real(WP) :: sigma = 10.0_WP
        real(WP) :: rho = 28.0_WP
        real(WP) :: beta = 8.0_WP/3.0_WP
    contains
        procedure :: F => LorenzDE
    end type LorenzSys

    contains

    function LorenzDE(me, X, Y, Params)
        intrinsic :: size
        class(LorenzSys), intent(inout) :: me !< Differential Equation object
        real(WP), intent(in) :: X
        real(WP), intent(in), dimension(me%n) :: Y
        real(WP), intent(in), dimension(:), optional :: Params

        real(WP), dimension(size(Y)) :: LorenzDE

        LorenzDE(1) = me%sigma*(Y(2) - Y(1))
        LorenzDE(2) = Y(1)*(me%rho - Y(3)) - Y(2)
        LorenzDE(3) = Y(1)*Y(2) - me%beta*Y(3)
    end function LorenzDE

end module LorenzDiffEq

    program LorenzTest

    use iso_fortran_env, only: output_unit
    use FLINT
    use LorenzDiffEq

    implicit none

    ! Simulation Parameters

    ! Number of periodic orbits to run and loops for benchmark
    integer, parameter :: nloops = 100
    integer, parameter :: floops = 100000000

    ! solvers to run from 1 to 4
    integer, parameter :: nmethod = 4

    ! scalar tolerance
    real(WP), parameter :: rtol   = 1.0e-11_WP
    real(WP), parameter :: atol   = rtol*1.0e-3_WP

    ! max no. of steps
    integer, parameter :: MAX_STEPS = 100000
    logical, parameter :: CONST_STEPSZ = .FALSE.

    ! random dispersion of IC: this multiplies the random number
    real(WP), parameter :: randon = 0.000000000000_WP

    ! Interpolation points
    integer, parameter :: nIp = int(1000.0)

    ! results file
    character(len=20) :: fname = 'results_Lorenz.txt'
    character(len=20) :: spfname = 'states_ip_Lorenz.txt'

    type(ERK_class) :: erkvar
    type(LorenzSys) :: LorenzObj

    real(WP) :: x0, xf, xfval, rnum
    real(WP), dimension(3) :: y0, yf, y0r

    character(len=:), allocatable :: errstr

    real(WP) :: stepsz0, ipdx, stepsz
    real(WP), dimension(:), allocatable, target :: Xint, Xarr
    real(WP), dimension(:,:), allocatable, target :: Yint, Yarr
    integer :: stiffstatus, stifftestval
    integer :: nsteps, naccpt, nrejct, fcalls, steps, itr, nevents, ctr
    character(len=10) :: mname
    integer :: method

    real(WP) :: T0, Tf

    real(WP), dimension(3), parameter :: params = [10.0_WP, 28.0_WP, 8.0_WP/3.0_WP]

    y0 = [1.0_WP, 0.0_WP, 0.0_WP]

    x0 = 0.0_WP
    xF = 100.0_WP

    LorenzObj%n = 3
    LorenzObj%m = 0
    ! LorenzObj%G => SampleEvent

    stepsz0 = 0.0E-2_WP    ! let FLINT compute the initial step-size
    stifftestval = 1*0  ! check for stiffness and stop integration if stiff

    ! arrays for interpolated data
    ipdx = (xf-x0)/(nIp-1)
    Xarr = [(x0+ipdx*itr,itr=0,(nIp-1))]
    Xarr(nIp) = xf
    allocate(Yarr(3,size(Xarr)))

    ! create results file
    open(unit=17,file= fname, status = 'replace')
    write(17, *) '--- FLINT Stats ---'
    write(17, *) ' '
    write(17, *) 'A. Internal Step-size'
    write(17, '(6A12)') 'Method', 'Time(s)', 'FCalls', 'Accepted', 'Rejected', 'Total Steps'

    ! states using interpolation
    open(unit=19,file= spfname, status = 'replace')
    write(19, '(4A12)') 'X','Y1','Y2','Y3'

    ! Solution at internal step size

    do itr = 1, nmethod

        select case (itr)
        case (1)
            mname = 'DOP54'
            method = ERK_DOP54
        case (2)
            mname = 'DOP853'
            method = ERK_DOP853
        case (3)
            mname = 'Verner65E'
            method = ERK_Verner65E
        case (4)
            mname = 'Verner98R'
            method = ERK_Verner98R
        end select

        call erkvar%Init(LorenzObj, MAX_STEPS, Method=method, ATol=[atol], RTol=[rtol],&
            InterpOn=.FALSE.)
        if (erkvar%status == FLINT_SUCCESS) then

            y0r = y0
            call CPU_TIME(t0)
            do ctr = 1,nloops
                stiffstatus = stifftestval
                stepsz = stepsz0
                xfval = xf
                call RANDOM_NUMBER(rnum) ! randomize initial condition
                y0r(1) = y0r(1) + rnum*randon*y0r(1)

                call erkvar%Integrate(x0, y0r, xfval, yf, StepSz=stepsz, UseConstStepSz=CONST_STEPSZ, &
                    IntStepsOn=.FALSE.,StiffTest=stiffstatus)
            end do
            call CPU_TIME(tf)

            if (stiffstatus == -1) write(17, *) mname//': problem is stiff'

            if (erkvar%status == FLINT_SUCCESS .OR. erkvar%status == FLINT_EVENT_TERM &
                        .OR. erkvar%status == FLINT_ERROR_MAXSTEPS) then
                call erkvar%Info(stiffstatus, errstr, nAccept=naccpt, nReject=nrejct, nFCalls=fcalls, nSteps = steps)

                write(17, '(A12,1F12.3,4I12.1)') mname, (tf-t0), fcalls, naccpt, nrejct, steps
            else
                call erkvar%Info(stiffstatus, errstr)
                write(17,*) mname//': Integrate failed: ', erkvar%status, ':', errstr
            end if

            ! ! write states
            ! if (itr == 2) then
            !     do ctr = 1, size(Xint)
            !         write(18, '(7E18.9)') Xint(ctr), Yint(1:3,ctr)
            !     end do
            ! end if

            if (allocated(Xint)) deallocate(Xint)
            if (allocated(Yint)) deallocate(Yint)
        else
            write(17,*) mname//': Init failed: ', erkvar%status
        end if
    end do

     ! Solution at interpolated grid
     write(17, *) 'B. Interpolated Grid'
     write(17, '(6A12)') 'Method', 'Time(s)', 'FCalls', 'Accepted', 'Rejected', 'Total Steps'

     do itr = 1, nmethod

         select case (itr)
         case (1)
             mname = 'DOP54'
             method = ERK_DOP54
         case (2)
             mname = 'DOP853'
             method = ERK_DOP853
         case (3)
             mname = 'Verner65E'
             method = ERK_Verner65E
         case (4)
             mname = 'Verner98R'
             method = ERK_Verner98R
         end select

         call erkvar%Init(LorenzObj, MAX_STEPS, Method=method, ATol=[atol], RTol=[rtol],&
             InterpOn=.TRUE.)

         if (erkvar%status == FLINT_SUCCESS) then

             y0r = y0
             call CPU_TIME(t0)
             do ctr = 1,nloops
                 stiffstatus = stifftestval
                 stepsz = stepsz0
                 xfval = xf
                 call RANDOM_NUMBER(rnum) ! randomize initial condition
                 y0r(1) = y0r(1) + rnum*randon*y0r(1)
                 call erkvar%Integrate(x0, y0r, xfval, yf, StepSz=stepsz, UseConstStepSz=CONST_STEPSZ, &
                     StiffTest=stiffstatus, params=params)
             end do
             call CPU_TIME(tf)
             call erkvar%Interpolate(Xarr,Yarr,.TRUE.)

             if (stiffstatus == -1) write(17, *) mname//': problem is stiff'

             if (erkvar%status == FLINT_SUCCESS .OR. erkvar%status == FLINT_EVENT_TERM &
                         .OR. erkvar%status == FLINT_ERROR_MAXSTEPS) then
                 call erkvar%Info(stiffstatus, errstr, nAccept=naccpt, nReject=nrejct, nFCalls=fcalls, nSteps=steps)

                 write(17, '(A12,1F12.3,4I12.1)') mname, (tf-t0), fcalls, naccpt, nrejct, steps

             else
                 call erkvar%Info(stiffstatus, errstr)
                 write(17,*) mname//': Integrate failed: ', erkvar%status, ':',errstr
             end if

             ! write states
             if (itr == 2) then
                 do ctr = 1, size(Xarr)
                     write(19, '(7E18.9)') Xarr(ctr), Yarr(1:3,ctr)
                 end do
             end if

         else
             write(17,*) mname//': Init failed: ', erkvar%status
         end if
     end do

    ! function test
    y0r = y0*0.0_WP
    call CPU_TIME(t0)
    do ctr = 1,floops
        !y0 = y0 + 0.0001_WP*y0
        yf = LorenzObj%F(x0, y0, params)
        y0r = y0r + yf
    end do
    call CPU_TIME(tf)
    write(17, *) ' '
    write(17, '(I12,A13,F12.2,A5)') floops, ' Func calls: ', (tf-t0)*1.0e3_WP, ' ms'

    contains

    end program LorenzTest

max_stars_repo_name: princemahajan/FLINT

max_issues_repo_name: princemahajan/FLINT

max_forks_repo_name: princemahajan/FLINT


text:
!-----------------------------------------------------------------------------
!> Calculate jacobian for Michelsen state matrices.
!>
!> \todo Need trace-component functionallity.
!>
!-----------------------------------------------------------------------------
module state_functions
  !
  !
  use thermopack_var, only: nc
  use thermopack_constants
  use eos
  use numconstants, only: almost_zero, machine_prec
  implicit none
  private
  save

  public :: getStateFuncMatrix, dhdt_twoPhase, dvdp_twoPhase, dvdt_twoPhase
  public :: getStateFunc, dhdp_twoPhase, dnvdX, getSVDerivativesTwoPhase
  public :: getUVDerivativesTwoPhase
  public :: getJouleThompsonCoeff, dpdt_twoPhase

contains

  !-----------------------------------------------------------------------------
  !> Calculate jacobian matrix for Michelsen state functions
  !> Implemented for two phase PH, PS, PT and UV flash
  !> Included TV and PV for convenience.
  !>
  !> \author MHA, 2012-03-20
  !-----------------------------------------------------------------------------
  subroutine getStateFuncMatrix(t,p,Z,betaV,betaL,X,Y,spec,vspec,M,RHS,phase,simpleMatrix)
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Liquid phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    character(len=2), intent(in) :: spec !< Character specification
    real, dimension(2), intent(in) :: vspec !< Specification value
    real, dimension(:,:), intent(out) :: M !< Jacobian matrix
    real, dimension(:), intent(out) :: RHS !< Right hand side
    integer, dimension(2), optional, intent(in) :: phase !< Phase integer
    logical, optional, intent(in) :: simpleMatrix !< Deactivate composition diff.
    ! Locals
    real, dimension(nc) :: FUGL,FUGV,FUGLT,FUGVT,FUGLP,FUGVP
    real, dimension(nc,nc) :: FUGLX,FUGVX
    real :: hv,hl,dhvdt,dhldt,cpmix,dhdp,dhvdp,dhldp
    real :: sv,sl,dsldt,dsvdt,smix,hmix,vv,vl,vmix,dvdt,dvdp
    real :: dvldp,dvldt,dvvdp,dvvdt
    real :: denom
    integer :: i, phaseX, phaseY
    logical :: deactivateCompDiff

    phaseX = LIQPH
    phaseY = VAPPH
    if (present(phase)) then
      phaseX = phase(1)
      phaseY = phase(2)
    endif

    if (spec == 'TP' .OR. spec == 'PT') then
      call thermo(t,p,X,phaseX,FUGL,lnFUGX=FUGLX) ! ln phi liquid
      call thermo(t,p,Y,phaseY,FUGV,lnFUGX=FUGVX) ! ln phi vapour
    else if (spec == 'SP' .OR. spec == 'PS' .OR. spec == 'HP' .OR. spec == 'PH') then
      call thermo(t,p,X,phaseX,FUGL,lnFUGX=FUGLX,lnfugt=FUGLT) ! ln phi liquid
      call thermo(t,p,Y,phaseY,FUGV,lnFUGX=FUGVX,lnfugt=FUGVT) ! ln phi vapour
      if (spec == 'SP' .OR. spec == 'PS') then
        call entropy(t,p,X,phaseX,sl,dsdt=dsldt)
        call entropy(t,p,Y,phaseY,sv,dsdt=dsvdt)
        smix = betaV*sv + betaL*sl
        RHS(1+nc) = (vspec(1) - smix)/Rgas
        cpmix = T*(betaV*dsvdt + betaL*dsldt)
      else if (spec == 'HP' .OR. spec == 'PH') then
        call enthalpy(t,p,X,phaseX,hl,dhdt=dhldt)
        call enthalpy(t,p,Y,phaseY,hv,dhdt=dhvdt)
        hmix = betaV*hv + betaL*hl
        RHS(1+nc) = (vspec(1) - hmix)/(Rgas*T)
        cpmix = betaV*dhvdt + betaL*dhldt
      endif
      M(1+nc,1+nc) = -cpmix/Rgas
      M(1:nc,1+nc) = T*(FUGVT - FUGLT)
      M(1+nc,1:nc) = T*(FUGVT - FUGLT)
    else if (spec == 'UV') then
      call thermo(t,p,X,phaseX,FUGL,lnFUGX=FUGLX,lnfugt=FUGLT,lnfugp=FUGLP) ! ln phi liquid
      call thermo(t,p,Y,phaseY,FUGV,lnFUGX=FUGVX,lnfugt=FUGVT,lnfugp=FUGVP) ! ln phi vapour

      ! Temperature terms
      call enthalpy(t,p,X,phaseX,hl,dhdt=dhldt)
      call enthalpy(t,p,Y,phaseY,hv,dhdt=dhvdt)
      hmix = betaV*hv + betaL*hl
      RHS(1+nc) = (vspec(1) + P*vspec(2) - hmix)/(Rgas*T)
      cpmix = betaV*dhvdt + betaL*dhldt
      M(1+nc,1+nc) = -cpmix/Rgas
      M(1:nc,1+nc) = T*(FUGVT - FUGLT)
      M(1+nc,1:nc) = T*(FUGVT - FUGLT)

      ! Pressure terms
      call specificVolume(t,p,X,phaseX,vl,dvdt=dvldt,dvdp=dvldp)
      call specificVolume(t,p,Y,phaseY,vv,dvdt=dvvdt,dvdp=dvvdp)
      vmix = betaV*vv + betaL*vl
      RHS(2+nc) = P*(vmix-vspec(2))/(Rgas*T)
      dVdp = betaV*dvvdp + betaL*dvldp
      dVdt = betaV*dvvdt + betaL*dvldt
      M(2+nc,2+nc) = p*p*dVdp/(Rgas*t)
      M(1+nc,2+nc) = p*dVdt/Rgas
      M(2+nc,1+nc) = p*dVdt/Rgas
      M(1:nc,2+nc) = p*(FUGVP - FUGLP)
      M(2+nc,1:nc) = p*(FUGVP - FUGLP)
    else if (spec == 'TV' .OR. spec == 'VT') then
      call thermo(t,p,X,phaseX,FUGL,lnFUGX=FUGLX,lnfugp=FUGLP) ! ln phi liquid
      call thermo(t,p,Y,phaseY,FUGV,lnFUGX=FUGVX,lnfugp=FUGVP) ! ln phi vapour
      call specificVolume(t,p,X,phaseX,vl,dvdp=dvldp)
      call specificVolume(t,p,Y,phaseY,vv,dvdp=dvvdp)
      dVdp = betaV*dvvdp + betaL*dvldp
      M(1+nc,1+nc) = p*p*dVdp/(Rgas*t)
      M(1:nc,1+nc) = p*(FUGVP - FUGLP)
      M(1+nc,1:nc) = p*(FUGVP - FUGLP)
      vmix = betaV*vv + betaL*vl
      RHS(1+nc) = P*(vmix-vspec(1))/(Rgas*T)
    else if (spec == 'PV' .OR. spec == 'VP') then
      ! NOT SYMETRIC!
      call thermo(t,p,X,phaseX,FUGL,lnFUGX=FUGLX,lnfugt=FUGLT,lnfugp=FUGLP) ! ln phi liquid
      call thermo(t,p,Y,phaseY,FUGV,lnFUGX=FUGVX,lnfugt=FUGVT,lnfugp=FUGVP) ! ln phi vapour
      call specificVolume(t,p,X,phaseX,vl,dvdt=dvldt)
      call specificVolume(t,p,Y,phaseY,vv,dvdt=dvvdt)
      dVdt = betaV*dvvdt + betaL*dvldt
      M(1+nc,1+nc) = p*dVdt/Rgas
      M(1:nc,1+nc) = T*(FUGVT - FUGLT)
      M(1+nc,1:nc) = p*(FUGVP - FUGLP)
      vmix = betaV*vv + betaL*vl
      RHS(1+nc) = P*(vmix-vspec(1))/(Rgas*T)

    else if (spec == 'SV') then
      call thermo(t,p,X,phaseX,FUGL,lnFUGX=FUGLX,lnfugt=FUGLT,lnfugp=FUGLP) ! ln phi liquid
      call thermo(t,p,Y,phaseY,FUGV,lnFUGX=FUGVX,lnfugt=FUGVT,lnfugp=FUGVP) ! ln phi vapour

      ! Entropy terms
      call entropy(t,p,X,phaseX,sl,dsdt=dsldt)
      call entropy(t,p,Y,phaseY,sv,dsdt=dsvdt)
      smix = betaV*sv + betaL*sl
      RHS(1+nc) = (vspec(1) - smix)/Rgas
      cpmix = T*(betaV*dsvdt + betaL*dsldt)
      M(1+nc,1+nc) = -cpmix/Rgas
      M(1:nc,1+nc) = T*(FUGVT - FUGLT)
      M(1+nc,1:nc) = T*(FUGVT - FUGLT)

      ! Pressure terms
      call specificVolume(t,p,X,phaseX,vl,dvdt=dvldt,dvdp=dvldp)
      call specificVolume(t,p,Y,phaseY,vv,dvdt=dvvdt,dvdp=dvvdp)
      vmix = betaV*vv + betaL*vl
      RHS(2+nc) = P*(vmix-vspec(2))/(Rgas*T)
      dVdp = betaV*dvvdp + betaL*dvldp
      dVdt = betaV*dvvdt + betaL*dvldt
      M(2+nc,2+nc) = p*p*dVdp/(Rgas*t)
      M(1+nc,2+nc) = p*dVdt/Rgas
      M(2+nc,1+nc) = p*dVdt/Rgas
      M(1:nc,2+nc) = p*(FUGVP - FUGLP)
      M(2+nc,1:nc) = p*(FUGVP - FUGLP)
    else if (spec == 'TH' .OR. spec == 'HT') then
      call thermo(t,p,X,phaseX,FUGL,lnFUGX=FUGLX,lnfugt=FUGLT,lnfugp=FUGLP) ! ln phi liquid
      call thermo(t,p,Y,phaseY,FUGV,lnFUGX=FUGVX,lnfugt=FUGVT,lnfugp=FUGVP) ! ln phi vapour
      call enthalpy(t,p,X,phaseX,hl,dhdp=dhldp)
      call enthalpy(t,p,Y,phaseY,hv,dhdp=dhvdp)
      hmix = betaV*hv + betaL*hl
      RHS(1+nc) = (vspec(1) - hmix)/(Rgas*T)
      dhdp = betaV*dhvdp + betaL*dhldp
      M(1+nc,1+nc) = -P*dhdp/(Rgas*T)
      M(1:nc,1+nc) = P*(FUGVP - FUGLP)
      M(1+nc,1:nc) = T*(FUGVT - FUGLT)
    endif

    ! If simpleMatrix==.true., set M(1:nc,1:nc) to the identity matrix. A NR
    ! iteration will then not change compositions.
    if (present(simpleMatrix)) then
      deactivateCompDiff = simpleMatrix
    else
      deactivateCompDiff = .false.
    end if
    if (deactivateCompDiff) then
      M(1:nc,1:nc) = 0.0
      do i=1,nc
        M(i,i) = 1.0
        if (Z(i) > zLimit) then
          RHS(i) = log(Y(i)) - log(X(i)) + FUGV(i) - FUGL(i)
        else
          RHS(i) = 0.0
        endif
      end do

      return
    end if

    if (betaV*betaL < almost_zero) then
      M(1:nc,1:nc) = 0.0
    else
      M(1:nc,1:nc) = 1.0/(betaV*betaL)*(-1.0 + betaV*FUGLX + betaL*FUGVX)
    endif

    ! Add to the diagonal
    do i=1,nc
      if (Z(i) > zLimit) then
        RHS(i) = log(Y(i)) - log(X(i)) + FUGV(i) - FUGL(i)
        denom = X(i)*Y(i)*betaV*betaL
        if(denom > almost_zero) M(i,i) = M(i,i) + Z(i)/denom
      else
        RHS(i) = 0.0
        M(:,i) = 0.0
        M(i,:) = 0.0
        M(i,i) = 1.0
      endif
    enddo

  end subroutine getStateFuncMatrix

  !-----------------------------------------------------------------------------
  !> Calculate RHS function values for Michelsen state functions
  !> Values corrensponds to subroutine getStateFuncMatirx
  !>
  !> \author MH, 2014-10-23
  !-----------------------------------------------------------------------------
  subroutine getStateFunc(t,p,Z,betaV,betaL,X,Y,spec,vspec,RHS,phase)
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Liquid phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    character(len=2), intent(in) :: spec !< Character specification
    real, dimension(2), intent(in) :: vspec !< Specification value
    real, dimension(:), intent(out) :: RHS !< Right hand side
    integer, dimension(2), optional, intent(in) :: phase !< Phase integer
    ! Locals
    real, dimension(nc) :: FUGL,FUGV
    real :: sv,sl,hv,hl,smix,hmix,vv,vl,vmix
    integer :: i, phaseX, phaseY

    phaseX = LIQPH
    phaseY = VAPPH
    if (present(phase)) then
      phaseX = phase(1)
      phaseY = phase(2)
    endif

    if (spec == 'TP' .OR. spec == 'PT') then
      !
    else if (spec == 'SP' .OR. spec == 'PS') then
      call entropy(t,p,X,phaseX,sl)
      call entropy(t,p,Y,phaseY,sv)
      smix = betaV*sv + betaL*sl
      RHS(1+nc) = (vspec(1) - smix)/Rgas
    else if (spec == 'HP' .OR. spec == 'PH') then
      call enthalpy(t,p,X,phaseX,hl)
      call enthalpy(t,p,Y,phaseY,hv)
      hmix = betaV*hv + betaL*hl
      RHS(1+nc) = (vspec(1) - hmix)/(Rgas*T)
    else if (spec == 'UV') then
      ! Temperature terms
      call enthalpy(t,p,X,phaseX,hl)
      call enthalpy(t,p,Y,phaseY,hv)
      hmix = betaV*hv + betaL*hl
      RHS(1+nc) = (vspec(1) + P*vspec(2) - hmix)/(Rgas*T)

      ! Pressure terms
      call specificVolume(t,p,X,phaseX,vl)
      call specificVolume(t,p,Y,phaseY,vv)
      vmix = betaV*vv + betaL*vl
      RHS(2+nc) = P*(vmix-vspec(2))/(Rgas*T)
    else if (spec == 'TV' .OR. spec == 'VT') then
      call specificVolume(t,p,X,phaseX,vl)
      call specificVolume(t,p,Y,phaseY,vv)
      vmix = betaV*vv + betaL*vl
      RHS(1+nc) = P*(vmix-vspec(1))/(Rgas*T)
    else if (spec == 'PV' .OR. spec == 'VP') then
      ! NOT SYMETRIC!
      call specificVolume(t,p,X,phaseX,vl)
      call specificVolume(t,p,Y,phaseY,vv)
      vmix = betaV*vv + betaL*vl
      RHS(1+nc) = P*(vmix-vspec(1))/(Rgas*T)

    else if (spec == 'SV') then
      ! Entropy terms
      call entropy(t,p,X,phaseX,sl)
      call entropy(t,p,Y,phaseY,sv)
      smix = betaV*sv + betaL*sl
      RHS(1+nc) = (vspec(1) - smix)/Rgas

      ! Pressure terms
      call specificVolume(t,p,X,phaseX,vl)
      call specificVolume(t,p,Y,phaseY,vv)
      vmix = betaV*vv + betaL*vl
      RHS(2+nc) = P*(vmix-vspec(2))/(Rgas*T)

    else if (spec == 'TH' .OR. spec == 'HT') then
      call enthalpy(t,p,X,phaseX,hl)
      call enthalpy(t,p,Y,phaseY,hv)
      hmix = betaV*hv + betaL*hl
      RHS(1+nc) = (vspec(1) - hmix)/(Rgas*T)
    endif

    call thermo(t,p,X,phaseX,FUGL) ! ln phi liquid
    call thermo(t,p,Y,phaseY,FUGV) ! ln phi vapour
    ! Add to the diagonal
    do i=1,nc
      if (Z(i) > zLimit) then
        RHS(i) = log(Y(i)) - log(X(i)) + FUGV(i) - FUGL(i)
      else
        RHS(i) = 0.0
      endif
    enddo

  end subroutine getStateFunc

  !-----------------------------------------------------------------------------
  !> Calculate two-phase heat capacity at constant pressure.
  !>
  !> \author MHA, 2012-03-20
  !-----------------------------------------------------------------------------
  function dhdt_twoPhase(t,p,Z,betaV,betaL,X,Y,ph) result(dhdt)
    use numconstants, only: small
    use eos, only: enthalpy, thermo
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Liquid phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    integer, dimension(2), optional, intent(in) :: ph !< Phase integers
    real :: dhdt
    ! Locals
    real, dimension(nc+1,nc+1) :: M !< Jacobian matrix
    real, dimension(nc+1) :: RHS !< Right hand side
    real, dimension(2) :: vspec
    integer, dimension(nc+1) :: ipiv
    real, dimension(3*nc) ::  work
    integer :: ifail, phase
    real, dimension(nc) :: FUG
    real :: h

    vspec = 0.0
    call getStateFuncMatrix(t,p,Z,betaV,betaL,X,Y,'PH',vspec,M,RHS,ph)
    RHS(nc+1) = -1.0/(Rgas*t)
    RHS(1:nc) = 0.0

    ! lapack routines
    ! symmetric indefinite matrix
    call dsysv('u', nc+1, 1, M, nc+1, ipiv, RHS, nc+1, work, 3*nc, ifail)

    RHS(nc+1) = t * RHS(nc+1) ! Convert from dlnT -> dT
    if (ifail /= 0 .or. RHS(nc+1) <= small .or. RHS(nc+1) /= RHS(nc+1)) then
      ! Invalid heat capacity (proximity of critical point)
      ! Use single phase value
      call thermo(t,p,z,MINGIBBSPH,FUG,ophase=phase)
      if (phase == SINGLEPH) then
        phase = LIQPH
      endif
      call enthalpy(t,p,z,phase,h,dhdt=dhdt)
    else
      dhdt = 1.0 / RHS(nc+1)
    endif
  end function dhdt_twoPhase

  !-----------------------------------------------------------------------------
  !> Calculate two-phase enthalpy change at constant temperature.
  !>
  !> \author MH, 2014
  !-----------------------------------------------------------------------------
  function dhdp_twoPhase(t,p,Z,betaV,betaL,X,Y,ph) result(dhdp)
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Liquid phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    integer, dimension(2), optional, intent(in) :: ph !< Phase integers
    real :: dhdp
    ! Locals
    real, dimension(nc+1,nc+1) :: M !< Jacobian matrix
    real, dimension(nc+1) :: RHS !< Right hand side
    real, dimension(2) :: vspec
    integer, dimension(nc+1) :: ipiv
    integer :: ifail, phase
    real, dimension(nc) :: FUG
    real :: h
    logical :: lufail

    vspec = 0.0
    call getStateFuncMatrix(t,p,Z,betaV,betaL,X,Y,'TH',vspec,M,RHS,ph)
    RHS(nc+1) = -1.0/(Rgas*t)
    RHS(1:nc) = 0.0

    ! lapack routines
    ! general lu
    lufail = .false.
    call dgetrf(nc+1,nc+1,M,nc+1,ipiv,ifail)
    if (ifail.ne.0) then
      lufail = .true.
    endif
    ! backsubstitute
    if (.not. lufail) then
      call dgetrs('n',nc+1,1,M,nc+1,ipiv,RHS,nc+1,ifail)
      if (ifail.ne.0) then
        lufail = .true.
      endif
    endif

    if (lufail .or. RHS(nc+1) /= RHS(nc+1)) then
      ! Use sinlgle phase values
      call thermo(t,p,z,MINGIBBSPH,FUG,ophase=phase)
      if (phase == SINGLEPH) then
        phase = LIQPH
      endif
      call enthalpy(t,p,z,phase,h,dhdp=dhdp)
    else
      RHS(nc+1) = P * RHS(nc+1) ! Convert from dlnP -> dP
      dhdp = 1.0 / RHS(nc+1)
    endif

  end function dhdp_twoPhase

  !-----------------------------------------------------------------------------
  !> Calculate two-phase temperature differential of specific volume at
  !! constant pressure.
  !!
  !! \author MH, 2012-07-06
  !-----------------------------------------------------------------------------
  function dvdt_twoPhase(t,p,Z,betaV,betaL,X,Y,ph) result(dvdt)
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Liquid phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    integer, dimension(2), optional, intent(in) :: ph !< Phase integers
    real :: dvdt
    ! Locals
    real, dimension(nc+1,nc+1) :: M !< Jacobian matrix
    real, dimension(nc+1) :: RHS !< Right hand side
    real, dimension(2) :: vspec
    integer, dimension(nc+1) :: ipiv
    integer :: ifail, phase
    logical :: lufail
    real, dimension(nc) :: FUG
    real :: v

    vspec = 0.0
    call getStateFuncMatrix(t,p,Z,betaV,betaL,X,Y,'PV',vspec,M,RHS,ph)
    RHS(nc+1) = p/(Rgas*t)
    RHS(1:nc) = 0.0

    ! lapack routines
    ! general lu
    lufail = .false.
    call dgetrf(nc+1,nc+1,M,nc+1,ipiv,ifail)
    if (ifail.ne.0) then
      lufail = .true.
    endif
    ! backsubstitute
    if (.not. lufail) then
      call dgetrs('n',nc+1,1,M,nc+1,ipiv,RHS,nc+1,ifail)
      if (ifail.ne.0) then
        lufail = .true.
      endif
    endif

    if (lufail .or. RHS(nc+1) /= RHS(nc+1) .or. &
         machine_prec > abs(RHS(nc+1))) then
      ! Use sinlgle phase values
      call thermo(t,p,z,MINGIBBSPH,FUG,ophase=phase)
      if (phase == SINGLEPH) then
        phase = LIQPH
      endif
      call specificVolume(t,p,z,phase,v,dvdt=dvdt)
    else
      RHS(nc+1) = t * RHS(nc+1) ! Convert from dlnT -> dT
      dvdt = 1.0 / RHS(nc+1)
    endif
  end function dvdt_twoPhase

  !-----------------------------------------------------------------------------
  !> Calculate two-phase pressure differential of specific volume at
  !! constant temperature.
  !!
  !! \author MH, 2012-07-06
  !-----------------------------------------------------------------------------
  function dvdp_twoPhase(t,p,Z,betaV,betaL,X,Y,ph) result(dvdp)
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Vapour phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    integer, dimension(2), optional, intent(in) :: ph !< Phase integers
    real :: dvdp
    ! Locals
    real, dimension(nc+1,nc+1) :: M !< Jacobian matrix
    real, dimension(nc+1) :: RHS !< Right hand side
    real, dimension(2) :: vspec
    integer, dimension(nc+1) :: ipiv
    real, dimension(3*nc) ::  work
    integer :: ifail, phase
    real, dimension(nc) :: FUG
    real :: v

    vspec = 0.0
    call getStateFuncMatrix(t,p,Z,betaV,betaL,X,Y,'TV',vspec,M,RHS,ph)
    RHS(nc+1) = p/(Rgas*t)
    RHS(1:nc) = 0.0

    ! lapack routines
    ! symmetric indefinite matrix
    call dsysv('u', nc+1, 1, M, nc+1, ipiv, RHS, nc+1, work, 3*nc, ifail)
    if (ifail /= 0 .or. RHS(nc+1) /= RHS(nc+1)) then
      ! Use sinlgle phase values
      call thermo(t,p,z,MINGIBBSPH,FUG,ophase=phase)
      if (phase == SINGLEPH) then
        phase = LIQPH
      endif
      call specificVolume(t,p,z,phase,v,dvdp=dvdp)
    else
      RHS(nc+1) = p * RHS(nc+1) ! Convert from dlnP -> dP
      dvdp = 1.0 / RHS(nc+1)
    endif
  end function dvdp_twoPhase

  !-----------------------------------------------------------------------------
  !> Calculate temperature/pressure differential of mole number (and phase fraction)
  !! at constant pressure
  !! \author MH, 2015-01
  !-----------------------------------------------------------------------------
  subroutine dnvdX(t,p,Z,betaV,betaL,X,Y,dbetadT,dbetadP,dndT,dndP)
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Vapour phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    real, optional, intent(out) :: dbetadT,dbetadP,dndT(nc),dndP(nc) !< Mole based
    ! Locals
    real, dimension(nc,nc) :: M !< Jacobian matrix
    real, dimension(nc,2) :: RHS !< Right hand side
    real, dimension(2) :: vspec
    integer, dimension(nc) :: ipiv
    real, dimension(3*nc) ::  work
    integer :: ifail
    real, dimension(nc) :: FUGL,FUGV,FUGLT,FUGVT,FUGLP,FUGVP

    vspec = 0.0
    call getStateFuncMatrix(t,p,Z,betaV,betaL,X,Y,'TP',vspec,M,RHS(:,1))
    call thermo(t,p,X,LIQPH,FUGL,lnfugt=FUGLT,lnfugp=FUGLP)
    call thermo(t,p,Y,VAPPH,FUGV,lnfugt=FUGVT,lnfugp=FUGVP)
    RHS(:,1) = -(FUGVT - FUGLT)
    RHS(:,2) = -(FUGVP - FUGLP)

    ! lapack routines
    ! symmetric indefinite matrix
    call dsysv('u', nc, 2, M, nc, ipiv, RHS, nc, work, 3*nc, ifail)
    if (ifail /= 0 .or. RHS(1,1) /= RHS(1,1)) then
      RHS = 0.0
    endif
    if (present(dbetadT)) then
      dbetadT = sum(RHS(:,1))
    endif
    if (present(dndT)) then
      dndT = RHS(:,1)
    endif
    if (present(dbetadP)) then
      dbetadP = sum(RHS(:,2))
    endif
    if (present(dndP)) then
      dndP = RHS(:,2)
    endif
  end subroutine dnvdX

  !-----------------------------------------------------------------------------
  !> Calculate single- or two-phase derivatives of entropy and/or molar volume.
  !>
  !> \author MAG, 2015-03-03
  !-----------------------------------------------------------------------------
  subroutine getSVDerivativesTwoPhase(t,p,z,betav,betal,x,y,iphase,dsdp_v,dsdt_v,dvdp_s,dvdt_s)
    implicit none
    real, intent(in) :: betav !< Vapour phase molar fraction [-]
    real, intent(in) :: betal !< Liquid phase molar fraction [-]
    real, dimension(nc), intent(in) :: z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: x !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    integer, intent(in) :: iphase !< Phase flag
    real, optional, intent(out) :: dsdp_v !< Partial derivative of entropy w.r.t. pressure at constant molar volume
    real, optional, intent(out) :: dsdt_v !< Partial derivative of entropy w.r.t. temperature at constant molar volume
    real, optional, intent(out) :: dvdp_s !< Partial derivative of molar volume w.r.t. pressure at constant entropy
    real, optional, intent(out) :: dvdt_s !< Partial derivative of molar volume w.r.t. temperature at constant entropy
    ! Locals
    real, dimension(nc+2,nc+2) :: a,m !< Jacobian matrix
    real, dimension(nc+2) :: rhs !< Right hand side
    real, dimension(2) :: vspec
    integer, dimension(nc+2) :: ipiv
    real :: s,v,dsdp_t,dsdt_p,dvdt_p,dvdp_t
    integer :: ifail,iphase_internal
    logical :: lufail
    !
    iphase_internal = iphase
    !
    if (iphase_internal == TWOPH) then
      ! Two-phase calculation
      !
      vspec = 0.0
      call getStateFuncMatrix(t,p,z,betav,betal,x,y,'SV',vspec,m,rhs)
      !
      if (present(dsdp_v) .or. present(dsdt_v)) then
        rhs(1:nc) = 0.0
        rhs(nc+1) = -1.0/rgas
        rhs(nc+2) = 0.0
        a = m
        ! Solve linear system
        lufail = .false.
        call dgetrf(nc+2,nc+2,a,nc+2,ipiv,ifail)
        if (ifail /= 0) lufail = .true.
        if (.not. lufail) then
          call dgetrs('n',nc+2,1,a,nc+2,ipiv,rhs,nc+2,ifail)
          if (ifail /= 0) lufail = .true.
        endif
        ! Set derivative values
        if (present(dsdt_v)) then
          rhs(nc+1) = t*rhs(nc+1) ! Convert from dlnT -> dT
          dsdt_v = 1.0/rhs(nc+1)
        end if
        if (present(dsdp_v)) then
          rhs(nc+2) = p*rhs(nc+2) ! Convert from dlnP -> dP
          dsdp_v = 1.0/rhs(nc+2)
        end if
      end if
      !
      if (present(dvdp_s) .or. present(dvdt_s)) then
        rhs(1:nc+1) = 0.0
        rhs(nc+2) = p/(rgas*t)
        a = m
        ! Solve linear system
        lufail = .false.
        call dgetrf(nc+2,nc+2,a,nc+2,ipiv,ifail)
        if (ifail /= 0) lufail = .true.
        if (.not. lufail) then
          call dgetrs('n',nc+2,1,a,nc+2,ipiv,rhs,nc+2,ifail)
          if (ifail /= 0) lufail = .true.
        endif
        ! Set derivative values
        if (present(dvdt_s)) then
          rhs(nc+1) = t*rhs(nc+1) ! Convert from dlnT -> dT
          dvdt_s = 1.0/rhs(nc+1)
        end if
        if (present(dvdp_s)) then
          rhs(nc+2) = p*rhs(nc+2) ! Convert from dlnP -> dP
          dvdp_s = 1.0/rhs(nc+2)
        end if
      end if
    else
      ! Single-phase calculation
      !
      if (iphase_internal /= VAPPH) iphase_internal = LIQPH
      call entropy(t,p,z,iphase_internal,s,dsdt=dsdt_p,dsdp=dsdp_t)
      call specificVolume(t,p,z,iphase_internal,v,dvdt=dvdt_p,dvdp=dvdp_t)
      !
      if (present(dsdp_v)) dsdp_v = dsdp_t - dsdt_p*dvdp_t/dvdt_p
      if (present(dsdt_v)) dsdt_v = dsdt_p - dsdp_t*dvdt_p/dvdp_t
      if (present(dvdp_s)) dvdp_s = dvdp_t - dvdt_p*dsdp_t/dsdt_p
      if (present(dvdt_s)) dvdt_s = dvdt_p - dvdp_t*dsdt_p/dsdp_t
    endif
    !
  end subroutine getSVDerivativesTwoPhase

  !-----------------------------------------------------------------------------
  !> Calculate single- or two-phase derivatives of internal energy and/or molar
  !> volume.
  !>
  !> CAUTION: The implementations of dudp_v and dudt_v are covered by the unit
  !>          tests in 3DMF. The implementations of dvdp_u and dvdt_u have not
  !>          been tested yet.
  !>
  !> \author MAG, 2015-11-09
  !-----------------------------------------------------------------------------
  subroutine getUVDerivativesTwoPhase(t,p,z,betav,betal,x,y,iphase, &
       dudp_v,dudt_v,dvdp_u,dvdt_u)
    use eosTV, only: pressure,internal_energy_tv
    implicit none
    real, intent(in) :: betav !< Vapour phase molar fraction [-]
    real, intent(in) :: betal !< Liquid phase molar fraction [-]
    real, dimension(nc), intent(in) :: z !< Overall molar composition [-]
    real, dimension(nc), intent(in) :: x !< Liquid molar composition [-]
    real, dimension(nc), intent(in) :: y !< Vapour molar composition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    integer, intent(in) :: iphase !< Phase flag
    real, optional, intent(out) :: dudp_v !< Partial derivative of internal energy w.r.t. pressure at constant molar volume
    real, optional, intent(out) :: dudt_v !< Partial derivative of internal energy w.r.t. temperature at constant molar volume
    real, optional, intent(out) :: dvdp_u !< Partial derivative of molar volume w.r.t. pressure at constant internal energy
    real, optional, intent(out) :: dvdt_u !< Partial derivative of molar volume w.r.t. temperature at constant internal energy
    ! Local variables:
    real, dimension(nc+2,nc+2) :: a,m !< Jacobian matrix
    real, dimension(nc+2) :: rhs !< Right hand side
    real, dimension(2) :: vspec
    integer, dimension(nc+2) :: ipiv
    real :: u_internal,v,dvdt_p,dvdp_t,dudt_v_internal
    real :: dudv_t,dudv_p,dudp_v_internal
    real :: dpdv_t,dpdv_u,dpdt_v,p_internal
    integer :: ifail,iphase_internal
    logical :: lufail
    !
    iphase_internal = iphase
    !
    if (present(dvdp_u) .or. present(dvdt_u)) then
      call stoperror("The implementations of the derivatives dvdp_u and dvdt_u have not been tested yet.")
    end if
    !
    if (iphase_internal == TWOPH) then
      !
      ! Two-phase calculation
      ! =====================
      !
      vspec = 0.0
      call getStateFuncMatrix(t,p,z,betav,betal,x,y,'UV',vspec,m,rhs)
      !
      if (present(dudp_v) .or. present(dudt_v)) then
        rhs(1:nc) = 0.0
        rhs(nc+1) = -1.0/(rgas*t)
        rhs(nc+2) = 0.0
        a = m
        ! Solve linear system
        lufail = .false.
        call dgetrf(nc+2,nc+2,a,nc+2,ipiv,ifail)
        if (ifail /= 0) lufail = .true.
        if (.not. lufail) then
          call dgetrs('n',nc+2,1,a,nc+2,ipiv,rhs,nc+2,ifail)
          if (ifail /= 0) lufail = .true.
        endif
        ! Set derivative values
        if (present(dudt_v)) then
          rhs(nc+1) = t*rhs(nc+1) ! Convert from dlnT -> dT
          dudt_v = 1.0/rhs(nc+1)
        end if
        if (present(dudp_v)) then
          rhs(nc+2) = p*rhs(nc+2) ! Convert from dlnP -> dP
          dudp_v = 1.0/rhs(nc+2)
        end if
      end if
      !
      if (present(dvdp_u) .or. present(dvdt_u)) then
        rhs(1:nc+1) = 0.0
        rhs(nc+2) = p/(rgas*t)
        a = m
        ! Solve linear system
        lufail = .false.
        call dgetrf(nc+2,nc+2,a,nc+2,ipiv,ifail)
        if (ifail /= 0) lufail = .true.
        if (.not. lufail) then
          call dgetrs('n',nc+2,1,a,nc+2,ipiv,rhs,nc+2,ifail)
          if (ifail /= 0) lufail = .true.
        endif
        ! Set derivative values
        if (present(dvdt_u)) then
          rhs(nc+1) = t*rhs(nc+1) ! Convert from dlnT -> dT
          dvdt_u = 1.0/rhs(nc+1)
        end if
        if (present(dvdp_u)) then
          rhs(nc+2) = p*rhs(nc+2) ! Convert from dlnP -> dP
          dvdp_u = 1.0/rhs(nc+2)
        end if
      end if
    else
      !
      ! Single-phase calculation
      ! ========================
      !
      if (iphase_internal /= VAPPH) iphase_internal = LIQPH
      call specificVolume(t,p,z,iphase_internal,v,dvdt=dvdt_p,dvdp=dvdp_t)
      call internal_energy_tv(t,v,z,u_internal,dudt=dudt_v_internal,dudv=dudv_t)
      p_internal = pressure(t,v,z,dpdv_t,dpdt_v)
      !
      if (present(dudt_v)) dudt_v = dudt_v_internal
      if (present(dvdt_u)) dvdt_u = -dudt_v_internal/dudv_t
      dudv_p = dudv_t - dudt_v_internal*dpdv_t/dpdt_v
      dpdv_u = dpdv_t - dpdt_v*dudv_t/dudt_v_internal
      dudp_v_internal = -dudv_p/dpdv_u
      if (present(dudp_v)) then
        dudp_v = dudp_v_internal
      end if
      if (present(dvdp_u)) dvdp_u = -dudp_v_internal/dudv_p
    endif
    !
  end subroutine getUVDerivativesTwoPhase

  !--------------------------------------------------------------------------
  !> Calculate two-phase Joule-Thompson coefficient
  !! dT/dP at constant enthalpy
  !! \author MH, 2015-09
  !--------------------------------------------------------------------------
  function getJouleThompsonCoeff(t,p,Z,betaV,betaL,X,Y,phase,ph) result(dTdP)
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Vapour phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    integer, intent(in) :: phase !< Phase identifyer
    integer, dimension(2), optional, intent(in) :: ph !< Phase integers
    real :: dTdP
    ! Locals
    real, dimension(nc+1,nc+1) :: M !< Jacobian matrix
    real, dimension(nc+1) :: RHS !< Right hand side
    real, dimension(2) :: vspec !< Dummy
    integer, dimension(nc+1) :: ipiv
    real, dimension(3*nc) ::  work
    integer :: ifail, lphase, phaseX, phaseY
    real, dimension(nc) :: lnfugl,lnfuglp,lnfugv,lnfugvp
    real :: hl,dhldp,hv,dhvdp,h,dhdt,dhdp

    lphase = phase
    if (phase == TWOPH) then
      vspec = 0.0
      call getStateFuncMatrix(t,p,Z,betaV,betaL,X,Y,'PH',vspec,M,RHS,ph)
      if (present(ph)) then
        phaseX = ph(1)
        phaseY = ph(2)
      else
        phaseX = LIQPH
        phaseY = VAPPH
      endif
      call enthalpy(t,p,X,phaseX,hl,dhdp=dhldp)
      call enthalpy(t,p,Y,phaseY,hv,dhdp=dhvdp)
      RHS(1+nc) = (betaV*dhvdp + betaL*dhldp)/(Rgas*T)
      call thermo(t,p,X,phaseX,lnfugl,lnfugp=lnfuglp)
      call thermo(t,p,Y,phaseY,lnfugv,lnfugp=lnfugvp)
      RHS(1:nc) = -(lnfugvp - lnfuglp)
      ! lapack routines
      ! symmetric indefinite matrix
      call dsysv('u', nc+1, 1, M, nc+1, ipiv, RHS, nc+1, work, 3*nc, ifail)
      if (ifail /= 0 .or. RHS(nc+1) /= RHS(nc+1)) then
        ! Use sinlgle phase values
        call thermo(t,p,z,MINGIBBSPH,lnfugl,ophase=lphase)
      else
        dTdP = T * RHS(nc+1) ! Convert from dlnT -> dT
      endif
    endif
    if (lphase /= TWOPH) then
      if (lphase == SINGLEPH) then
        lphase = LIQPH
      endif
      call enthalpy(t,p,Z,lphase,h,dhdt=dhdt,dhdp=dhdp)
      dTdP = -dhdp/dhdt
    endif
  end function getJouleThompsonCoeff

  !--------------------------------------------------------------------------
  !> Calculate temperaure differential of pressure at
  !! constant volume.
  !!
  !! \author MH, 2015-10
  !--------------------------------------------------------------------------
  function dpdt_twoPhase(t,p,Z,betaV,betaL,X,Y,phase,ph) result(dpdt)
    use eos, only: specificVolume, thermo
    implicit none
    real, intent(in) :: betaV !< Vapour phase molar fraction [-]
    real, intent(in) :: betaL !< Vapour phase molar fraction [-]
    real, dimension(nc), intent(in) :: Z !< Overall molar compozition [-]
    real, dimension(nc), intent(in) :: X !< Liquid molar compozition [-]
    real, dimension(nc), intent(in) :: Y !< Vapour molar compozition [-]
    real, intent(in) :: t !< Temperature [K]
    real, intent(in) :: p !< Pressure [Pa]
    integer, intent(in) :: phase !< Phase identifyer
    integer, dimension(2), optional, intent(in) :: ph !< Phase integers
    real :: dpdt
    ! Locals
    real, dimension(nc+1,nc+1) :: M !< Jacobian matrix
    real, dimension(nc+1) :: RHS !< Right hand side
    real, dimension(2) :: vspec
    integer, dimension(nc+1) :: ipiv
    real, dimension(3*nc) ::  work
    integer :: ifail, lphase, phaseX, phaseY
    real, dimension(nc) :: FUGL, FUGV, lnfugvt, lnfuglt
    real :: v, vv, vl, dvvdt, dvldt, dvdt, dvdp

    lphase = phase
    if (phase == TWOPH) then
      vspec = 0.0
      call getStateFuncMatrix(t,p,Z,betaV,betaL,X,Y,'TV',vspec,M,RHS,ph)
      if (present(ph)) then
        phaseX = ph(1)
        phaseY = ph(2)
      else
        phaseX = LIQPH
        phaseY = VAPPH
      endif
      call thermo(t,p,X,phaseX,FUGL,lnfugt=lnfuglt) ! ln phi liquid
      call thermo(t,p,Y,phaseY,FUGV,lnfugt=lnfugvt) ! ln phi vapour
      call specificVolume(t,p,X,phaseX,vl,dvdt=dvldt)
      call specificVolume(t,p,Y,phaseY,vv,dvdt=dvvdt)
      dVdt = betaV*dvvdt + betaL*dvldt
      RHS(nc+1) = dvdt*p/(Rgas*t)
      RHS(1:nc) = (lnfugvt - lnfuglt)
      RHS = - RHS
      ! lapack routines
      ! symmetric indefinite matrix
      call dsysv('u', nc+1, 1, M, nc+1, ipiv, RHS, nc+1, work, 3*nc, ifail)
      if (ifail /= 0 .or. RHS(nc+1) /= RHS(nc+1)) then
        ! Use sinlgle phase values
        call thermo(t,p,z,MINGIBBSPH,FUGL,ophase=lphase)
      else
        dpdt = p * RHS(nc+1) ! Convert from dlnP -> dP
      endif
    endif
    if (lphase /= TWOPH) then
      if (lphase == SINGLEPH) then
        lphase = LIQPH
      endif
      call specificVolume(t,p,z,lphase,v,dvdt=dvdt,dvdp=dvdp)
      dpdt = -dvdt/dvdp
    endif
  end function dpdt_twoPhase

end module state_functions

max_stars_repo_name: SINTEF/Thermopack

max_issues_repo_name: SINTEF/Thermopack

max_forks_repo_name: SINTEF/Thermopack


text:
      subroutine dump_geometries_lg(iunit,iomode,ierror)
C
C
C#######################################################################
C
C     PURPOSE -
C
C    This routine writes the geometry information, surfaces,
C    regions and mregions to a file
C
C
C     INPUT ARGUMENTS -
C
C     iunit  unit number of file
C
C     OUTPUT ARGUMENTS -
C
C
C
C     CHANGE HISTORY -
C
C        $Log: dump_geometries_lg.f,v $
C        Revision 2.00  2007/11/05 19:45:52  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.3   22 Mar 2002 16:41:42   dcg
CPVCS    allow for a mesh object with no geometry
CPVCS
CPVCS       Rev 1.2   07 Feb 2000 16:49:22   dcg

C#######################################################################
C
      implicit none
c
      include 'geom_lg.h'
c
      integer iunit,i,j,ns,nr,nm,ls,lr,lm,len,ierror,length,
     *   icscode
      character *(*) iomode
      character*32 isubname,partname

C
      isubname='dumplggm'
      ierror=0

c
c  write out active info then list of geometry names
c  and general info for each geom
c
      partname='geom_lg'
      call mmfindbk ('geom_names',partname,
     *  ipgeom_names,len,ierror)
      call mmfindbk ('geom_info',partname,
     *  ipgeom_info,len,ierror)
      if(iomode(1:5).eq.'ascii') then
         write(iunit,10) number_of_geometries,nsurf,nregs,
     *   nmregs,maxdef,maxmdef,lastregdef,lastmregdef,
     *   lastsparam
         write(iunit,12)(geom_names(i),i=1,number_of_geometries)
         do i=1,number_of_geometries
           write(iunit,10)(geom_info(j,i),j=1,8)
         enddo
 10      format(9i10)
 12      format(4a32)

      else
         write(iunit) number_of_geometries,nsurf,nregs,
     *   nmregs,maxdef,maxmdef,lastregdef,lastmregdef,
     *   lastsparam,(geom_names(i),i=1,number_of_geometries),
     *   ((geom_info(j,i),j=1,8),
     *     i=1,number_of_geometries)
      endif
c
c  loop through the geometries
c
      do i=1,number_of_geometries
      partname=geom_names(i)
      if(i.ne.number_of_geometries) then
         ns=geom_info(1,i)
         nr=geom_info(2,i)
         nm=geom_info(3,i)
         lr=geom_info(6,i)
         lm=geom_info(7,i)
         ls=geom_info(8,i)
      else
         ns=nsurf
         nr=nregs
         nm=nmregs
         lr=lastregdef
         lm=lastmregdef
         ls=lastsparam
      endif
c
c surface info for this geom - name, type, boundary type,
c definition, offset to definition
c
      call mmfindbk('csall',partname,ipcsall,length,ierror)
      if(ierror.ne.0) ns=0
      call mmfindbk('istype',partname,ipistype,length,ierror)
      call mmfindbk('ibtype',partname,ipibtype,length,ierror)
      call mmfindbk('sheetnm',partname,ipsheetnm,length,ierror)
      call mmfindbk('surfparam',partname,ipsurfparam,length,ierror)
      call mmfindbk('offsparam',partname,ipoffsparam,length,ierror)
      if(iomode(1:5).eq.'ascii') then
         if(ns.gt.0) then
            write(iunit,12) (csall(j),j=1,ns)
            write(iunit,12) (istype(j),j=1,ns)
            write(iunit,12) (ibtype(j),j=1,ns)
            write(iunit,12) (sheetnm(j),j=1,ns)
            write(iunit,20) (surfparam(j),j=1,ls)
 20         format(6e20.13)
            write(iunit,10) (offsparam(j),j=1,ns)
         endif
      else
         if(ns.gt.0) then
            write(iunit)(csall(j),j=1,ns),(istype(j),j=1,ns),
     *      (ibtype(j),j=1,ns), (sheetnm(j),j=1,ns),
     *      (surfparam(j),j=1,ls),(offsparam(j),j=1,ns)
         else
            write(iunit)
         endif
      endif

C
C  write region storage blocks - name , offset to definition,
c  length of definition, definitions
c
      call mmfindbk('cregs',partname,ipcregs,length,ierror)
      call mmfindbk('offregdef',partname,ipoffregdef,length,ierror)
      call mmfindbk('ndefregs',partname,ipndefregs,length,ierror)
      call mmfindbk('regdef',partname,ipregdef,length,ierror)
      if(iomode(1:5).eq.'ascii') then
         if(nr.gt.0) then
            write(iunit,12) (cregs(j),j=1,nr)
            write(iunit,10) (offregdef(j),j=1,nr)
            write(iunit,10) (ndefregs(j),j=1,nr)
            write(iunit,12) (regdef(j),j=1,lr)
         endif
      else
         if(ns.gt.0)
     *     write(iunit)(cregs(j),j=1,nr),(offregdef(j),j=1,nr),
     *    (ndefregs(j),j=1,nr), (regdef(j),j=1,lr)
      endif

C
C  write mregion storage blocks
      call mmfindbk('cmregs',partname,ipcmregs,length,ierror)
      call mmfindbk('offmregdef',partname,ipoffmregdef,length,ierror)
      call mmfindbk('ndefmregs',partname,ipndefmregs,length,ierror)
      call mmfindbk('mregdef',partname,ipmregdef,length,ierror)
      call mmfindbk('matregs',partname,ipmatregs,length,ierror)
      if(iomode(1:5).eq.'ascii') then
         if(nm.gt.0) then
            write(iunit,12) (cmregs(j),j=1,nm)
            write(iunit,10) (offmregdef(j),j=1,nm)
            write(iunit,10) (ndefmregs(j),j=1,nm)
            write(iunit,12) (mregdef(j),j=1,lm)
            write(iunit,10) (matregs(j),j=1,nm)
         endif
      else
         if(nm.gt.0)
     *    write(iunit)(cmregs(j),j=1,nm),(offmregdef(j),j=1,nm),
     *    (ndefmregs(j),j=1,nm), (mregdef(j),j=1,lm),
     *    (matregs(j),j=1,nm)
      endif
      enddo

C
 9999 call mmrelprt(isubname,icscode)
      return
      end

C

max_forks_repo_name: daniellivingston/LaGriT


text:
C
C
*+ make_phbm

       subroutine make_phbm(beam,m,n,ibc,jbc,alpha,status)
C      ---------------------------------------------------
C
C modify the beam for a Prussian Hat clean
C
C Input/Returned
C   BEAM      -      R4(array)       -     beam to modify
C
C Input
C   M/N       -      I4              -     Beam size
C   I/JBC     -      R4              -     max on beam
C   ALPHA     -      R4              -     fractional hat height
*-
       integer    status, m, n, ibc, jbc

       real*4     beam(m,n), alpha

       if (status.ne.0) return

       beam(ibc,jbc) = beam(ibc,jbc) + alpha*beam(ibc,jbc)

       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
	SUBROUTINE GQCVSC  ( dvcvsc, iret )
C************************************************************************
C* GQCVSC								*
C*									*
C* This subroutine returns the curve scaling factor of a device driver.	*
C*									*
C* GQCVSC  ( DVCVSC, IRET )						*
C*									*
C* Output parameters:							*
C*	DVCVSC		REAL		Curve scaling factor		*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* T. Lee/GSC		 7/98						*
C************************************************************************
	INCLUDE		'FUNCCODE.PRM'
	INCLUDE		'ERROR.PRM'
C*
	INTEGER		isend (2)
C------------------------------------------------------------------------
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend (1) = 2
	isend (2) = FQCVSC
	CALL GPUT  ( isend, 2, iret )
	IF  ( iret .ne. NORMAL )  RETURN
C
C*	Get output parameters.
C
	CALL GGET  ( iret, 1, ierr )
	IF  ( ierr .ne. NORMAL )  THEN
	    iret   = ierr
	    RETURN
	END IF
C
	CALL GGETR  ( dvcvsc, 1, ierr )
	IF  ( ierr .ne. NORMAL )  iret = ierr
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine initgr(nowcrp)
c
c     + + + PURPOSE + + +
c     Called to initialize perennial growth after plant has been frozen
c     and is starting to grow.  Estimates the vegetative dry matter at
c     maturity (VDMT).
c
c     Called from PTGRP, INIT1
c     Author(s):  Alberts, Ghiddey, Ferris, Arnold
c     Reference in User Guide:
c
c     Changes:
c
c     Version: This module recoded from WEPP Version 92.25.
c     Date recoded: 08/25/92.
c     Recoded by: Charles R. Meyer.
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
      include 'pntype.inc'
      include 'pmxpln.inc'
      include 'pmxcrp.inc'
      include 'pmxres.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
      include 'pmxhil.inc'
      include 'pmxnsl.inc'
      include 'pmxcut.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
      integer nowcrp
c
c     + + + ARGUMENT DEFINITIONS + + +
c     nowcrp - the current crop
c
c     + + + COMMON BLOCKS + + +
      include 'ccover.inc'
c       read: cancov(mxplan)
c      write: canhgt(mxplan)
c
      include 'ccrpvr2.inc'
c     modify: vdmt(mxplan)
c
      include 'ccrpvr3.inc'
c       read: hmax(ntype), gddmaxi(ntype), bb(ntype), bbb(ntype),
c             xmxlai(ntype)
c     modify: sumgdd(mxplan)
c
      include 'ccrpprm.inc'
c       read: itype(mxcrop,mxplan)
c
      include 'ccrpout.inc'
c     modify: lai(mxplan)
c
      include 'cperen.inc'
c       read: imngmt(mxcrop,mxplan)
c
      include 'cstruc.inc'
c       read: iplane
c
c     + + + END SPECIFICATIONS + + +
c
c
c ---- Compute vegetative dry matter (VDMT), from canopy cover (CANCOV).
      vdmt(iplane) = log(1.0-cancov(iplane)) / (-
     1    bb(itype(nowcrp,iplane)))
      if (vdmt(iplane).lt.0.) vdmt(iplane) = 0.0
c
c     ---- Compute canopy height (CANHGT), from vegetative dry matter (VDMT)
c     and plant height at maturity (HMAX).
      canhgt(iplane) = (1.-exp(-bbb(itype(nowcrp,iplane))*vdmt(iplane)))
     1    * hmax(itype(nowcrp,iplane))
c
c     ---- Compute leaf area index (LAI), from maximum leaf area index
c     (XMXLAI) and vegetative dry matter (VDMT).
c     ** CROPLAND ANNUALS **
      if (imngmt(nowcrp,iplane).eq.1) then
        lai(iplane) = (xmxlai(itype(nowcrp,iplane))*vdmt(iplane)) / (
     1      vdmt(iplane)+0.5512*exp(-6.8*vdmt(iplane)))
c
c     ** PERENNIALS OR RANGELAND **
      else
        lai(iplane) = (xmxlai(itype(nowcrp,iplane))*vdmt(iplane)) / (
     1      vdmt(iplane)+0.2756*exp(-13.6*vdmt(iplane)))
      end if
c
c     ---- Compute cumulative growing degree days (SUMGDD), from growing
c     degree days at maturity (GDDMAX), leaf area index (LAI), and
c     maximum leaf area index (XMXLAI).
      sumgdd(iplane) = gddmax(itype(nowcrp,iplane)) * lai(iplane) /
     1    xmxlai(itype(nowcrp,iplane))
c
c     Handle the case of a winter annual which has relatively low
c     cover -- so low that calculation of sumgdd is below value of
c     the critical value of growing degree days for emergence (CRIT).
c
      if ((jdharv(nowcrp,iplane).lt.jdplt(nowcrp,iplane)).and.(
     1    sumgdd(iplane).lt.crit(itype(nowcrp,iplane))))
     1    sumgdd(iplane) = crit(itype(nowcrp,iplane)) + 0.1
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
! :: ----------------------------------------------------------
! :: summass
! ::             MASS = sum{ vol(i,j)*rho(i,j) }
! ::
! :: INPUTS / OUTPUTS:
! ::  rho        => density field
! ::  rlo,rhi    => index limits of rho aray
! ::  lo,hi      => index limits of grid interior
! ::  dx         => cell size
! ::  mass      <=  total mass
! ::  tmp        => temp column array
! :: ----------------------------------------------------------
! ::

       subroutine summass(rho,r_l1,r_l2,r_h1,r_h2,lo,hi,dx,mass)

       use prob_params_module, only : coord_type

       implicit none
       integer r_l1,r_l2,r_h1,r_h2
       integer lo(2), hi(2)
       double precision mass, dx(2)
       double precision rho(r_l1:r_h1,r_l2:r_h2)
       double precision  tmp(lo(2):hi(2))

       integer i, j
       double precision vol

       vol = dx(1) * dx(2)

       do j = lo(2),hi(2)
          tmp(j) = 0.d0
       enddo

       do i = lo(1), hi(1)
          do j = lo(2), hi(2)
             tmp(j) = tmp(j) + vol * rho(i,j)
          enddo
       enddo

       mass = 0.d0
       do j = lo(2), hi(2)
          mass = mass + tmp(j)
       enddo

       end subroutine summass

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


text:
subroutine printf(t,icompteur,y,ybio)
!*************************************
    use geoclim_write_output_mod, only: geoclim_write_output
    use geographic_write_output_mod, only: geographic_write_output
    use biodiv_write_output_mod,  only: biodiv_write_output
    use dynsoil_write_output_mod, only: dynsoil_write_output
    use dynsoil_offline_computation_mod, only: dynsoil_offline_computation
    use dynsoil_create_restart_mod, only: dynsoil_create_restart
    implicit none
    include 'combine_foam.inc'

!    compteur = compteur + 1
    icompteur = icompteur + 1
    if (float(icompteur).ge.xjump) then   !.or.isolver.eq.2) then
        icompteur=0
!        compteur=0

        !!!!! write geoclim outputs !!!!!
        call geoclim_write_output(t)

        !!!!! write biodiv outputs !!!!!
        if (coupling_ecogeo) then
            print *, 'ECOGEO module not available'
            stop
        end if

        !!!!! check for killing_signal !!!!!
        call read_deathnote(t)
        ! if killing signal is activated, tend <- t-10d6 and ageYprint <- t-10d6  =>  restart will be created, if not already done

    end if

    !!!!! write geographically-distributed outputs !!!!!
    icount_geogprint = icount_geogprint + 1
    if (icount_geogprint >= ijump_geogprint) then
        icount_geogprint = 0
        call geographic_write_output( GEOG_ofile_num, GEOG_ofile_name, GEOG_varout_name,  nlon, nlat, t, &
                                      Tclim, runclim, wth_allsil, wth_litho_wgh, wth_litho, fker, POC_export_rate, fp )
    end if

    !!!!! write dynsoil outputs (and create dynsoil offline variables): !!!!!
    if (coupling_dynsoil) then
        icount_DS_pri = icount_DS_pri + 1
        if ( icount_DS_pri >= ijump_DS_print ) then
            icount_DS_pri=0
            call dynsoil_offline_computation(xlevs, reg_z_prof,reg_tau_prof, reg_tau_surf,reg_x_surf,reg_thick,reg_eros, &
                                             Tclim,runclim,veget_factor, reg_ktop, reg_P_vol,                            &
                                             DS_missingpoints, DS_varout_missval(9),DS_varout_missval(10),               &
                                             reg_x_mean,reg_mean_age                                                     )
            call dynsoil_write_output(DS_ofile_num,DS_ofile_name,DS_varout_name, nlon,nlat, t, Tclim,runclim,  &
                                      reg_thick,reg_x_surf,reg_tau_surf, reg_z_prof,reg_tau_prof,              &
                                      reg_prod,reg_eros,reg_P_diss,reg_P_eros,reg_x_surf_eros,                 &
                                      reg_x_mean,reg_mean_age, reg_Li_Friv,reg_Li_Fsp,reg_Li_driv              )
        end if
    end if

!!!! TEST FOR RESTART FILES: !!!!

    if (t.gt.ageYprint.and.fog.eq.0) then
        fog=1
        ! extensive variables (chemical amounts) => concentration
        do j = 1,11*nbasin
            write(10,*) y(j)/vol(j)
        end do
        ! PCO2 and isopotic variables
        do j = 1+11*nbasin, 19*nbasin
            write(10,*) y(j)
        end do
        ! rest of extensive variables => concentration
        do j = 1+19*nbasin,nvar
            write(10,*) y(j)/vol(j)
        end do
        ! ecological module restart
        if (coupling_ecogeo) then
            print *, 'ECOGEO module not available'
            stop
        end if
        ! DYNSOIL restart files:
        if (coupling_dynsoil)  then
            call dynsoil_create_restart(output_path, run_name, DS_restart_name,                                                 &
                                        DS_varout_name, DS_varout_units, DS_varout_missvalname, DS_varout_missval,              &
                                        xlevs,ref_x_axis,ref_y_axis, reg_thick,reg_x_surf,reg_tau_surf,reg_z_prof,reg_tau_prof, &
                                        t, tsta                                                                                )
        end if
    endif

!   do i=1,nbasin
!       bilan_water(i)=0.
!   enddo

!   do i=1,nbasin
!       do j=1,nbasin
!           bilan_water(i)=bilan_water(i)-F(i,j) &
!                          +F(j,i)
!       end do
!   enddo

!   do i=1,nbasin
!       write(*,*)float(i),bilan_water(i)/(1.d+6*31.536d+6)
!   end do
!   stop

end subroutine


text:
subroutine c_filenames()
#include "cpp_options.h"

#ifdef one_file_per_step

!Need to hard code the file names below for one file per step scenario.
!Pay attention to the calculation of time increment

    use global, only: filenames,Nrecs,filename_increment,FnPartiInit
    implicit none
    character(len=5) :: fnn
    character(len=5) :: variables(6)
    integer*8 :: i,j

    print*, trim(FnPartiInit)
    !Use space to make sure " " has the same length as declared above, len=5 in this case. (There should be a clever way doing
    !this.)

    variables=(/"U    ","V    ","W    ","T    ","S    ","G    "/)

    !The variable 'filenames' is a two-dimensional array contaning the filenames for all variables at all steps.
    !The first dimension in filenames(i,j) corresponds to the number of files/steps for each varaible. Dimension size = Nrecs.
    !The second dimension corresponds to varialbe defined in 'variables' above. We have 6 variables here corresponding to
    !U/V/W/T/S/G where 'G' is neutral density. It was used in a paper related to DIMES where we release particles on a neutral
    !density surface. You can leave fn_GAMMA ='', G will not be used.
    DO i = 1, Nrecs
        write(fnn,"(I5.5)") i*filename_increment
        DO j = 1, 6
            !filenames(i,j) = trim(variables(j))//'/'//trim(variables(j))//"_90x300x640_"//trim(fnn)//'.data'
            filenames(i,j) = "rot_"//trim(variables(j))//"_"//trim(fnn)//'.bin'
        ENDDO
    ENDDO
    print*, "filename(1,1) c_filename",trim(filenames(1,1))

#else
    use global, only: fn_UVEL,fn_VVEL,fn_THETA,fn_WVEL,fn_SALT,fn_GAMMA, filenames,Nrecs,filename_increment,FnPartiInit

    filenames(1,1)=fn_UVEL
    filenames(1,2)=fn_VVEL
    filenames(1,3)=fn_WVEL
    filenames(1,4)=fn_THETA
    filenames(1,5)=fn_SALT
    filenames(1,6)=fn_GAMMA

#endif

end subroutine c_filenames

max_stars_repo_name: jinbow/Octupus

max_issues_repo_name: jinbow/Octupus

max_forks_repo_name: jinbow/Octupus


text:
! Generated by make_vegetable_driver. DO NOT EDIT
program main
    implicit none

    call run()
contains
    subroutine run()
        use compiler_test, only: &
                compiler_ref_reference => &
                    test_ref_reference
        use usage_test, only: &
                usage_usage => &
                    test_usage
        use vegetables, only: test_item_t, test_that, run_tests

        type(test_item_t) :: tests
        type(test_item_t) :: individual_tests(2)

        individual_tests(1) = compiler_ref_reference()
        individual_tests(2) = usage_usage()
        tests = test_that(individual_tests)

        call run_tests(tests)

    end subroutine
end program

max_stars_repo_name: wyphan/reference-counter

max_issues_repo_name: wyphan/reference-counter

max_forks_repo_name: wyphan/reference-counter


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine resultsnoddir(nk,v,nactdof,b,ipompc,nodempc,coefmpc,
     &  nmpc,mi)
!
!     copying the dof-values into (idir,node) format
!     (for sensitivity results; SPC's correspond to zero)
!
      implicit none
!
      integer mi(*),nactdof(0:mi(2),*),ipompc(*),nodempc(3,*),nk,i,j,
     &  nmpc,ist,ndir,node,index
!
      real*8 v(0:mi(2),*),b(*),coefmpc(*),bnac,fixed_disp
!
!     copying the dof-values into (idir,node) format
!
      do i=1,nk
         do j=1,mi(2)
            if(nactdof(j,i).gt.0) then
               bnac=b(nactdof(j,i))
            else
               cycle
            endif
            v(j,i)=bnac
         enddo
      enddo
!
!     treating the MPC's
!
      do i=1,nmpc
         ist=ipompc(i)
         node=nodempc(1,ist)
         ndir=nodempc(2,ist)
         index=nodempc(3,ist)
         fixed_disp=0.d0
         if(index.ne.0) then
            do
               fixed_disp=fixed_disp-coefmpc(index)*
     &              v(nodempc(2,index),nodempc(1,index))
               index=nodempc(3,index)
               if(index.eq.0) exit
            enddo
         endif
         fixed_disp=fixed_disp/coefmpc(ist)
         v(ndir,node)=fixed_disp
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
c { dg-do compile }
      Subroutine FndSph(Alpha,Rad)
      Dimension Rad(100),RadInp(100)
      Do I = 1, NSphInp
        Rad(I) = RadInp(I)
        Alpha = 1.2
      End Do
      End


! { dg-do run }
! Tests the fix for PR33888, in which the character length of
! the elemental function myfunc was not being calculated before
! the temporary for the array result was allocated.
!
! Contributed by Walter Spector <w6ws@earthlink.net>
!
program ftn95bug
  implicit none

  character(8) :: indata(4) =  &
              (/ '12344321', '98766789', 'abcdefgh', 'ABCDEFGH' /)

  call process (myfunc (indata))  ! <- This caused a gfortran ICE !

contains

  elemental function myfunc (s)
    character(*), intent(in) :: s
    character(len (s)) :: myfunc

    myfunc = s

  end function

  subroutine process (strings)
    character(*), intent(in) :: strings(:)

    if (any (strings .ne. indata)) STOP 1

  end subroutine

end program


text:
      PROGRAM TOP RANK	!Just do it.
      CHARACTER*28 HEADING(4)	!The first line is a header.
      TYPE GRIST		!But this is the stuff.
       CHARACTER*28 NAME		!Arbitrary sizes.
       CHARACTER*6 ID			!Possibly imperfect.
       REAL*8 SALARY			!Single precision is a bit thin.
       CHARACTER*6 DEPARTMENT		!Not a number.
      END TYPE GRIST		!So much for the aggregate.
      INTEGER HORDE		!Some parameterisation.
      PARAMETER (HORDE = 66)	!This should suffice.
      TYPE(GRIST) EMPLOYEE(HORDE),HIC	!An extra for the sort.
      LOGICAL CURSE		!Possible early completion.
      INTEGER I,N,H		!Steppers.
      INTEGER II,R		!Needed for the results.
      INTEGER KBD,MSG,IN	!I/O unit numbers.

      KBD = 5	!Standard input.
      MSG = 6	!Standard output.
      IN = 10	!Suitable for an input file.
      WRITE (MSG,1)
    1 FORMAT ("Reads a set of employee information from TopRank.csv"/
     1"Then for each department code, shows the N highest paid.")
      OPEN (IN,NAME = "TopRank.csv",FORM = "FORMATTED")
      READ (IN,*) HEADING	!Column headings: the "Salary" heading is not numeric.
Chug through the input.
      N = 0		!None so far.
   10 READ (IN,*,END = 20) EMPLOYEE(N + 1)	!Get the next record.
      N = N + 1					!We did. Count it in.
      IF (N.GT.HORDE) STOP "Too many employee records!"	!Ah, suspicion.
      GO TO 10					!Perhaps there will be another.
Collate the collection.
   20 CLOSE(IN)			!Finished with the input.
Crank up a comb sort, which requires only one comparison statement. Especially good for compound fields.
      H = N - 1			!Last - first, and not + 1.
   21 H = MAX(1,H*10/13)	!The special feature.
      CURSE = .FALSE.		!So far, so good.
      DO I = N - H,1,-1		!If H = 1, this is a Bubblesort.
        IF (EMPLOYEE(I).DEPARTMENT.LT.EMPLOYEE(I + H).DEPARTMENT) CYCLE	!In order by department.
        IF (EMPLOYEE(I).DEPARTMENT.EQ.EMPLOYEE(I + H).DEPARTMENT	!Or, Equal department,
     *  .AND. EMPLOYEE(I).SALARY.GT.EMPLOYEE(I + H).SALARY) CYCLE	!And in decreasing order by salary.
        CURSE = .TRUE.			!No escape. the elements are in the wrong order.
        HIC = EMPLOYEE(I)		!So a SWAP statement would be good. But alas.
        EMPLOYEE(I) = EMPLOYEE(I + H)	!For large data aggregates, an indexed sort would be good.
        EMPLOYEE(I + H) = HIC		!But, just slog away.
      END DO				!Consider the next pairing.
      IF (CURSE .OR. H.GT.1) GO TO 21	!Work remains?

Cast forth results.
   30 WRITE (6,31) N	!Announce, and solicit a parameter.
   31 FORMAT (I0," employees. How many per dept? ",$)	!The $ sez "don't start a new line."
      READ (KBD,*) R	!The parameter.
      IF (R.LE.0) STOP	!bah.
      WRITE (MSG,32) "Rank",HEADING
   32 FORMAT (/,A6,1X,A28,2X,A12,1X,A10,A)	!Compare to FORMAT 33.
      HIC.DEPARTMENT = "...Not this"	!Different from all departmental codes.
      DO I = 1,N	!Scan the sorted data.
        IF (HIC.DEPARTMENT.EQ.EMPLOYEE(I).DEPARTMENT) THEN	!Another?
          II = II + 1		!Same department, so count another adherent.
         ELSE			!But with a change of department code,
          II = 1		!Start a fresh count.
          HIC.DEPARTMENT = EMPLOYEE(I).DEPARTMENT	!And remember the new code.
          WRITE (MSG,*) "For ",HIC.DEPARTMENT	!Announce the new department's code.
        END IF			!So much for grouping.
        IF (II.LE.R) WRITE (MSG,33) II,EMPLOYEE(I)	!Still within the desired rank?
   33   FORMAT (I6,1X,A28,1X,A12,F11.2,1X,A)	!Some layout. Includes the repeated departmental code name.
      END DO			!On to the next.

      END	!That was straightforward.

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
      SUBROUTINE SSG1
C
      INTEGER         PG,SLT,BGPDT,CSTM,SIL,ECPT,MPT,GPTT,EDT,CASECC,
     1                CORE(166),SYSBUF,IWORD(4),MCB(7),SUBNAM(2)
      DIMENSION       PG(7),ILIST(360),ARY(1),DEFML(2),IDEFML(2),
     1                IARY(1),GVECT(1080)
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /LOADX / LC,SLT,BGPDT,OLD,CSTM,SIL,ISIL,ECPT,MPT,GPTT,EDT,
     1                N(3),LODC,MASS,NOBLD,IDIT,DUM6(6)
      COMMON /BLANK / NROWSP,LOADNN
      COMMON /SYSTEM/ SYSBUF,IOTPE,DUM53(53),ITHERM
      COMMON /LOADS / NLOAD,IPTR
      COMMON /ZZZZZZ/ ICORE(1)
      EQUIVALENCE     (CORE(1),ICORE(1),IARY(1),ARY(1)),
     1                (DEFML(1),IDEFML(1))
      DATA    IWORD / 4,6,7,162/
      DATA    SUBNAM/ 4HSSG1,4H    /
C
C     MODIFY OPEN CORE POINTER IPTR FOR MAGNETICS PROBLEM
C
      IPTR  = MAX0(NROWSP,166)
      MCB(1)= 105
      CALL RDTRL (MCB(1))
      IF (MCB(1) .GT. 0) IPTR = MAX0(3*NROWSP,3*MCB(2),166)
C
C     INITIALIZE.
C
      LC    = KORSZ(ICORE(1))
      NLLST = LC - 2*SYSBUF
      SLT   = 101
      BGPDT = 102
      CSTM  = 103
      SIL   = 104
      ECPT  = 105
      MPT   = 106
      GPTT  = 107
      EDT   = 108
      MASS  = 109
      CASECC= 110
      IDIT  = 111
      LODC  = 201
C             205 = NEWSLT (THERMAL)
      PG(1) = 301
      ICR2  = 302
      ICR3  = 303
      DO 10 I = 2,7
   10 PG(I) = 0
      PG(3) = NROWSP
      PG(4) = 2
      PG(5) = 1
C
C     AVOID CALCULATING UNUSED LOADS
C
C     NEDT  = NUMBER OF ELEMENT  DEFORMATIONS
C     NTEMP = NUMBER OF THERMAL LOADS
C     NCENT = NUMBER OF CENTRIFUGAL LOADS
C
      CALL SSG1A (N1,ILIST(1),NEDT,NTEMP,NCENT,CASECC,IHARM)
      N1A = N1 + 1
      LC  = LC - SYSBUF
      CALL OPEN (*310,PG(1),ICORE(LC+1),1)
      CALL WRITE (PG(1),PG(1),2,1)
      NGRAV = 0
      NEX   = N1 + NTEMP + NEDT + NCENT
      IF (N1 .EQ. 0) GO TO 21
C
C     MODIFY SLT -QVOL-, -QBDY1-, -QBDY2-, AND -QVECT- CARDS.
C
      NEWSLT = ICR3
      IF (ITHERM .NE. 0) NEWSLT = 205
      ISLT = SLT
      CALL SSGSLT (SLT,NEWSLT,ECPT)
      SLT = NEWSLT
      CALL EXTERN (NEX,NGRAV,GVECT(1),ILIST(1),PG(1),N1,IHARM)
C
C     RESET -SLT- TO ORIGINAL SLT DATA BLOCK
C
      SLT = ISLT
      N1  = N1 - NGRAV
   21 IF (NTEMP) 30,40,30
   30 CALL TEMPL (NTEMP,ILIST(N1+1),PG(1))
      N1 = N1 + NTEMP
   40 IF (NEDT) 50,60,50
   50 CALL EDTL (NEDT,ILIST(N1+1),PG(1))
      N1 = N1 + NEDT
   60 CALL CLOSE (PG,1)
      CALL WRTTRL (PG(1))
      IF (NGRAV) 90,100,90
   90 CONTINUE
C
C     CHECK TO SEE IF THE MASS MATRIX IS PURGED
C
      MCB(1) = MASS
      CALL RDTRL (MCB(1))
      IF (MCB(1) .LE. 0) CALL MESAGE (-56,0,IWORD)
      CALL GRAVL1 (NGRAV,GVECT(1),ICR2,IHARM)
C
C     USE LOAD FILE AS SCRATCH NOTHING ON IT NOW
C
      CALL SSG2B (MASS,ICR2,0,ICR3,0,1,1,LODC)
      CALL GRAVL2 (NGRAV,ICR3,PG(1))
      N1 = N1 + NGRAV
  100 IPONT1 = IPTR + 2
      IPONT  = IPTR + 1
      NLOAD  = 0
      DO 110 I = 1,NLLST
  110 IARY(I) = 0
      CALL OPEN (*320,CASECC,ICORE(LC+1),0)
      LC1  = LC - SYSBUF
      ISLT = 0
      CALL OPEN (*130,SLT,ICORE(LC1+1),0)
      ISLT = 1
      DO 120 I = 1,N1A
      CALL FWDREC (*270,SLT)
  120 CONTINUE
  130 DO 140 I = 1,LOADNN
      CALL FWDREC (*320,CASECC)
  140 CONTINUE
      IFRST = 0
  150 CALL READ (*250,*250,CASECC,CORE(1),166,1,FLAG)
      IF (IFRST .NE. 0) GO TO 151
      IFRST = 1
      ISPCN = CORE(3)
      MPCN  = CORE(2)
  151 CONTINUE
C
C     TEST FOR SYMMETRY, BUCKLING OR DIFFERENTIAL STIFFNESS.
C
      IF (CORE(16).NE.0 .OR. CORE(5).NE.0 .OR. CORE(138).NE.0) GO TO 150
      IF (CORE(3).NE.ISPCN .OR. CORE(2).NE.MPCN) GO TO 250
      INULL = 0
      DO 230 K = 1,4
      I = IWORD(K)
      IF (ITHERM.NE.0 .AND. I.EQ.7) GO TO 230
      IF (CORE(I) .EQ. 0) GO TO 230
      DO 160 J = 1,N1
      IF (CORE(I) .EQ. ILIST(J)) GO TO 220
  160 CONTINUE
C
C     COMBINATION CARD
C
      INULL = 1
  170 CALL READ (*270,*330,SLT,IDEFML(1),2,0,IFLAG)
      IF (CORE(I) .EQ. IDEFML(1)) GO TO 190
  180 CALL READ (*270,*330,SLT,IDEFML(1),2,0,IFLAG)
      IF (IDEFML(2) .EQ. -1) GO TO 170
      GO TO 180
  190 A = DEFML(2)
  200 CALL READ (*270,*330,SLT,IDEFML(1),2,0,IFLAG)
      IF (IDEFML(2) .EQ.  -1) GO TO 210
      IF (IPONT+1 .GT. NLLST) GO TO 340
      IARY(IPONT  ) = IARY(IPONT) + 1
      IARY(IPONT1 ) = IDEFML(2)
      ARY(IPONT1+1) = A*DEFML(1)
      IPONT1 = IPONT1 + 2
      GO TO 200
  210 CALL BCKREC (SLT)
      GO TO 230
  220 IARY(IPONT) = IARY(IPONT) + 1
      IF (IPONT+1 .GT. NLLST) GO TO 340
      IARY(IPONT1 ) = CORE(I)
      ARY(IPONT1+1) = 1.0
      IPONT1 = IPONT1 + 2
      INULL  = 1
  230 CONTINUE
      IF (INULL .EQ. 0) GO TO 260
  240 IPONT = IPONT + IARY(IPONT)*2 + 1
      NLOAD = NLOAD + 1
      IPONT1= IPONT1+ 1
      GO TO 150
  250 CALL CLOSE (CASECC,1)
      IF (ISLT .EQ. 1) CALL CLOSE (SLT,1)
      CALL COMBIN (PG(1),ILIST(1),N1)
      RETURN
C
  260 IARY(IPONT) = 1
      IF (IPONT+1 .GT. NLLST) GO TO 340
      IARY(IPONT1 ) =-1
      ARY(IPONT1+1) = 1.0
      IPONT1 = IPONT1 + 2
      GO TO 240
C
  270 IP1 = SLT
  280 IP2 =-1
  290 CALL MESAGE (IP2,IP1,SUBNAM)
      IP1 = CASECC
      GO TO 280
  310 IP1 = PG(1)
      GO TO 280
  320 IP1 = CASECC
      GO TO 280
  330 IP2 =-2
      IP1 = SLT
      GO TO 290
C
  340 I = ICORE(I)
      NWDS = 0
  350 CALL READ (*330,*360,SLT,CORE(1),LC,0,IFLAG)
      NWDS = NWDS + LC
      GO TO 350
  360 NWDS = NWDS + IFLAG
      WRITE  (IOTPE,370) UFM,I,NLLST,NWDS
  370 FORMAT (A23,' 3176, INSUFFICIENT OPEN CORE AVAILABLE TO PROCESS ',
     1       'ALL LOAD CARD COMBINATIONS IN MODULE SSG1.',
     2       /32X,'CURRENT LOAD ID BEING PROCESSED IS',I9,1H.,
     3       /32X,'OPEN CORE AVAILABLE IS',I9,' WORDS.',
     4       /32X,'ADDITIONAL OPEN CORE REQUIRED IS',I9,' WORDS.')
      IP1 = 0
      IP2 =-61
      GO TO 290
      END


text:
Nicholas Nick Park is a candidate for ASUCD ASUCD Senate Senate in the Fall 2013 ASUCD Election running on the NOW slate.

 Candidate Statement

Hey there Aggies!

My name is Nick Park and I’m excited to be running for Senate this quarter with the NOW Slate. I am a secondyear Managerial Economics major from Pleasanton, California and I am determined to make UC Davis the best campus possible for both current and future Aggies. I currently serve as an intern to ASUCD Senator Felicia Ong and ASUCD President Carly Sandstrom, and have worked on initiatives such as the 5K Stride for Aggie Pride and the ASUCD Scholarship.

I want to be your next ASUCD Senator because I am passionate about improving campus and bettering the Aggie experience!

IMPROVE ORIENTATION FOR INCOMING STUDENTS
Orientation is an important time for incoming students to learn about our campus, get excited to attend UC Davis, and make new friends as they transition to college. Currently, orientation falls short in fulfilling these goals. I want to enhance firstyear and transfer orientation by planning interactive tours of campus, utilizing the quad and downtown for social activities, and getting students the information they need through online webinars instead of fourhour long lectures in SciLec. I plan to work in partnership with Student Housing to accomplish this and ensure that our future Aggies leave orientation excited to return to campus in the fall.

INCREASE CALFRESH AID TO STUDENTS
The CalFresh program currently offers funds for food assistance to students approved for a work study award or to those who meet certain income and work requirements. However, most students who qualify for CalFresh are unaware of their eligibility. I want to increase awareness of the CalFresh Program to eligible UC Davis students by working with the ASUCD Pantry and the campus Financial Aid office so that students are aware of the program and its benefits. Ensuring that every student is afforded the opportunity to an education without sacrificing their personal wellbeing is a cause that I am deeply committed to.

I look forward to serving you and enhancing the college experience for both future and current Aggies! Please remember to VOTE NICK PARK #1 FOR ASUCD SENATE AND NOW #26!


text:
      block data bkdeqp
c
c     This is the block data subroutine for EQPT.
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
      end


text:
! sample compile command to generate *.o object file:
!      gfortran -fopenmp -c Fig_11.6_hist.f90

program main
   use omp_lib
   implicit none

   integer, parameter :: num_trials = 1000000	 ! number of x values
   integer, parameter :: num_bins = 100		   ! number of bins in histogram
   real*8, save :: xlow = 0.0;		   ! low end of x range
   real*8, save :: xhi = 10.0;		   ! high end of x range

   real*8 :: x
   integer*8 :: hist(num_bins)   ! the histogram
   integer*8 :: ival, i
   real*8 :: bin_width    ! the width of each bin in the hi, stogram
   real*8 :: sumh, sumhsq, ave, std_dev
   integer(kind=omp_lock_kind) :: hist_lcks(num_bins)   ! array of locks, one per bucket

   interface
      function drandom() result(val)
         real*8 :: val
      end function
      subroutine seed(low_in, hi_in)
         real*8, intent(in) :: low_in, hi_in
      end subroutine
   end interface

   call seed(xlow, xhi)     ! seed random generator over range of x
   bin_width = (xhi - xlow) / dble(num_bins)

   ! initialize the histogram and the array of locks
   !$omp parallel do schedule(static)
      do i = 1, num_bins
         hist(i) = 0
         call omp_init_lock(hist_lcks(i))
      enddo
   !$omp end parallel do
   ! test uniform pseudorandom sequence by assigning values
   ! to the right histogram bin
   !$omp parallel do schedule(static) private(x,ival)
      do i = 1, num_trials
         x = drandom()
         ival = int8((x - xlow)/bin_width)
         ! protect histogram bins. Low overhead due to uncontended locks
         call omp_set_lock(hist_lcks(ival))
         hist(ival) = hist(ival) + 1
         call omp_unset_lock(hist_lcks(ival))
      enddo
   !$omp end parallel do

   sumh = 0.0
   sumhsq = 0.0
   ! compute statistics (ave, std_dev) and destroy locks
   !$omp parallel do schedule(static) reduction(+:sumh,sumhsq)
      do i = 1, num_bins
         sumh  = sumh + hist(i)
         sumhsq = sumhsq + hist(i)*hist(i)
         call omp_destroy_lock(hist_lcks(i))
      enddo
   !$omp end parallel do

   ave = sumh / dble(num_bins)
   std_dev = sqrt(sumhsq /(dble(num_bins) - ave * ave))
end program main

max_stars_repo_name: zafar-hussain/OmpCommonCore

max_issues_repo_name: zafar-hussain/OmpCommonCore

max_forks_repo_name: zafar-hussain/OmpCommonCore