text:
!> Functionality for locating points on the saturation curve based on
!> auxiliary properties such as specific volume or entropy.
!> The isentropeEnvelopeCross routines should eventually be moved from
!> saturation to saturation_point_locators
!>
!> \author Ailo 2015-12
!> \author MH 2018-06
module saturation_point_locators ! Give a more descriptive name..
  use saturation
  use saturation_curve
  use solid_saturation
  use eos, only: thermo, entropy, specificVolume
  use thermopack_constants
  use thermopack_var, only: nc, nph, get_active_thermo_model, thermo_model
  use nonlinear_solvers
  use numconstants, only: machine_prec
  use puresaturation, only: puresat
  implicit none
  save

  integer, parameter :: locate_from_entropy = 1
  integer, parameter :: locate_from_lnvol = 2
  integer, parameter :: locate_from_enthalpy = 3
  integer, parameter :: locate_from_temperature = 4
  integer, parameter :: locate_from_pressure = 5

contains

  !-----------------------------------------------------------------------------
  !> Bracket solver for finding the exact point on the phase envelope
  !>
  !> What property is is determined by propflag.
  !>
  !> \author MH 2013-10, Ailo 2015-12
  !-----------------------------------------------------------------------------
  subroutine bracketSolveForProperty(n,Z,beta,propflag,propspec,X,Xold,dXdS,s,ierr,mode)
    implicit none
    integer, intent(in) :: n
    real, intent(in) :: Z(nc)
    real, intent(in) :: beta
    integer, intent(in) :: propflag
    real, intent(in) :: propspec
    real, intent(out) :: X(n)
    real, intent(in) :: Xold(n)
    real, intent(in) :: dXdS(n)
    integer, intent(in) :: s
    integer, optional, intent(out) :: ierr
    integer, optional, intent(in) :: mode
    ! Locals
    real, dimension(nc+3*n+5) :: param
    type(nonlinear_solver) :: solver
    real :: Xs,Xsmax,Xsmin
    !
    if (verbose) then
      print *,'In bracketSolveForProperty....'
    endif
    param(1) = propspec
    param(2) = beta
    param(3) = real(s)
    param(4:nc+3) = Z
    param(nc+4:nc+n+3) = Xold
    param(nc+n+4:nc+2*n+3) = dXdS
    param(nc+2*n+4:nc+3*n+3) = X
    param(nc+3*n+4) = real(propflag)
    if (present(mode)) then
      param(nc+3*n+5) = real(mode)
    endif
    solver%abs_tol = 1e-8
    solver%max_it = 30
    if (X(s) > Xold(s)) then
      Xsmax = X(s)
      Xsmin = Xold(s)
    else
      Xsmax = Xold(s)
      Xsmin = X(s)
    endif
    Xs = X(s)
    if (n == nc + 2) then
      call pegasus(Xsmin,Xsmax,propertyFunctionWrapper,Xs,solver,param)
    else if (n == 2) then
      call pegasus(Xsmin,Xsmax,propertyFunctionWrapperSubl,Xs,solver,param)
    else if (n == nc + 3) then
      call pegasus(Xsmin,Xsmax,propertyFunctionWrapperTriple,Xs,solver,param)
    endif
    if (present(ierr)) then
      ierr = solver%exitflag
    else
      if (solver%exitflag /= 0) then
        write(*,*) "Exitflag: ",solver%exitflag
        call stoperror("bracketSolveForProperty failed.")
      endif
    endif
    X = param(nc+2*n+4:nc+3*n+3)
  end subroutine bracketSolveForProperty

  !-----------------------------------------------------------------------------
  !> Bracket solver for finding the exact point on saturation line
  !>
  !> What property is is determined by propflag.
  !>
  !> \author MH 2019-05
  !-----------------------------------------------------------------------------
  subroutine bracketSolveForPropertySingle(ic,propflag,propspec,phase,T1,P1,dPdT1,Ts,Ps,ierr)
    implicit none
    integer, intent(in) :: ic, phase
    real, intent(in) :: T1,P1,dPdT1
    real, intent(inout) :: Ts,Ps
    integer, intent(in) :: propflag
    real, intent(in) :: propspec
    integer, optional, intent(out) :: ierr
    ! Locals
    real, dimension(9) :: param
    type(nonlinear_solver) :: solver
    real :: Xs,Xsmax,Xsmin
    !
    if (verbose) then
      print *,'In bracketSolveForPropertySingle....'
    endif
    param(1) = propspec
    param(2) = real(phase)
    param(3) = real(ic)
    param(4) = dPdT1
    param(5) = T1
    param(6) = P1
    param(7) = Ts
    param(8) = real(propflag)
    param(9) = 0.0
    solver%abs_tol = 1e-8
    Xsmax = max(Ps,P1)
    Xsmin = min(Ps,P1)
    Xs = P1
    call pegasus(Xsmin,Xsmax,propertyFunctionWrapperSingle,Xs,solver,param)
    if (present(ierr)) then
      ierr = solver%exitflag
    else
      if (solver%exitflag /= 0) then
        write(*,*) "Exitflag: ",solver%exitflag
        call stoperror("bracketSolveForProperty failed.")
      endif
    endif
    Ps = Xs
    Ts = param(9)
  end subroutine bracketSolveForPropertySingle

  !-----------------------------------------------------------------------------
  !> Function used to solve for point on saturation line having property value
  !> propspec.
  !>
  !> \author MH 2019-05
  !-----------------------------------------------------------------------------
  function propertyFunctionWrapperSingle(Xs,param) result(fun)
    implicit none
    real, intent(inout) :: Xs
    real, dimension(9), intent(inout) :: param
    real :: fun
    ! Locals:
    integer :: phase
    real, dimension(nc) :: Z
    real, dimension(1) :: X, Xmin, Xmax
    real :: t, prop, propspec, p, T0, P0, dPdT, Tmax
    integer :: propflag, ic
    type(nonlinear_solver) :: solver
    real :: paramSat(4)

    ! Unpack the param vector
    propspec = param(1)
    phase = nint(param(2))
    ic = nint(param(3))
    dPdT = param(4)
    T0 = param(5)
    P0 = param(6)
    Tmax = param(7)
    propflag = nint(param(8))

    ! Composition
    Z = 0
    Z(ic) = 1

    ! Extrapolate for better initial values
    p = Xs
    T = T0 + (P-P0)/dPdT
    X = log(T)

    ! Find point on sat line corresponding to P equal to Xs
    paramSat(1) = specP
    paramSat(3) = Xs
    paramSat(2) = real(ic)
    paramSat(4) = 0.0
    Xmin = log(min(T0,Tmax))
    Xmax = log(max(T0,Tmax))
    X = min(max(X,Xmin),Xmax)
    solver%rel_tol = 1.0e-20
    solver%abs_tol = 1.0e-10
    solver%max_it = 100
    call nonlinear_solve(solver,sat_fun_single_if,sat_diff_single,&
         sat_diff_single,limit_dx,premReturn,setXv,X,xmin,xmax,paramSat)
    t = exp(X(1))

    ! Update t-slot in param vector with new t value
    param(9) = t

    ! Calculate the property value at the new saturation point
    call genericProperty(t,p,Z,phase,propflag,prop)

    ! Compute new value of objective function
    fun = (prop - propspec)/max(abs(propspec), 1.0)
  end function propertyFunctionWrapperSingle

  !-----------------------------------------------------------------------------
  !> Function used to solve for point on phase envelope having property value
  !> propspec.
  !>
  !> \author MH 2013-10, Ailo 2015-12
  !-----------------------------------------------------------------------------
  function propertyFunctionWrapper(Xs,param) result(fun)
    implicit none
    real, intent(inout) :: Xs
    real, dimension(4*nc+10), intent(inout) :: param
    real :: fun
    ! Locals:
    integer :: iter, s, phase
    real, dimension(nc) :: Z,K
    real, dimension(nc+2) :: X,Xold,dXdS
    real :: t, prop, propspec, beta, ln_spec, p
    integer :: propflag

    ! Unpack the param vector
    propspec = param(1)
    beta = param(2)
    s = int(param(3))
    Z = param(4:nc+3)
    Xold = param(nc+4:2*nc+5)
    dXdS = param(2*nc+6:3*nc+7) ! is constant each time function is called
    X = param(3*nc+8:4*nc+9)
    propflag = int(param(4*nc+10))

    ! Extrapolate for better initial values in sat_newton
    X = Xold + dXdS*(Xs-Xold(s))
    K = exp(X(1:nc))
    t = exp(X(nc+1))
    p = exp(X(nc+2))

    ! Find point on envelope corresponding to ln(property value) equal to Xs
    ln_spec = Xs
    iter = sat_newton(Z,K,t,p,beta,s,ln_spec)
    X(1:nc) = log(K)
    X(nc+1) = log(T)
    X(nc+2) = log(p)

    ! Update Xold-slot in param vector with new X value
    param(3*nc+8:4*nc+9) = X

    ! Get correct phase
    if (beta > 0.5) then
      phase = VAPPH
    else
      phase = LIQPH
    endif

    ! Calculate the property value at the new saturation point
    call genericProperty(t,p,Z,phase,propflag,prop)

    ! Compute new value of objective function
    fun = (prop - propspec)/max(abs(propspec), 1.0)
  end function propertyFunctionWrapper

  subroutine genericProperty(t,p,Z,phase,propflag,prop)
    use eos
    use trend_solver, only: trend_density
    !$ use omp_lib
    implicit none
    integer, intent(in) :: phase !< Phase identifier
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: p !< Pa - Pressure
    real, dimension(1:nc), intent(in) :: Z !< Compozition
    integer, intent(in) :: propflag !< Flag determining what property to return
    real, intent(out) :: prop !< Property

    if (propflag == locate_from_entropy) then
      call entropy(t,p,Z,phase,prop)
    else if (propflag == locate_from_lnvol) then
      call specificVolume(t,p,Z,phase,prop)
      prop = log(prop)
    else if (propflag == locate_from_enthalpy) then
      call enthalpy(t,p,Z,phase,prop)
    else if (propflag == locate_from_temperature) then
      prop = T
    else if (propflag == locate_from_pressure) then
      prop = P
    end if

  end subroutine genericProperty

  !-----------------------------------------------------------------------------
  !> Locate the points on the phase envelope corresponding to the property values
  !> in prop_grid. The grid must be sorted according to whether one
  !> starts on the dew curve or the bubble curve.
  !>
  !> The routine starts at the dew point at given initial pressure, and then
  !> traverses the phase envelope until it has bracketed a grid value; it then
  !> locates the exact saturation point corresponding to this value, before
  !> continuing.
  !>
  !> This routine can perhaps be generalized to subsume the routines
  !> envelopePlot and isentropeEnvelopeCross.
  !>
  !> To do: test the routine when starting from bubble curve.
  !>
  !> \author Ailo 2015-12
  !> \author Morten 2018-06
  !-----------------------------------------------------------------------------
  subroutine sat_points_based_on_prop(Z,T0,P0,x0,y0,n_grid,propflag,prop_grid,&
       T_grid,P_grid,phase_grid,wi_grid,n_grid_found,&
       dS_override,phase_in,ierr_out,sgn_in,spec_in,&
       tMin,tMax,pMin,pMax,normal_grid,sequential_mode)
    use thermo_utils, only: isSingleComp
    implicit none
    ! Input:
    real, intent(in) :: Z(nc)                  !< total composition
    real, intent(in) :: t0, p0                 !< init. point has pressure p0
    real, dimension(nc), intent(in) :: x0,y0   !< init. liq./vap. compo. guess
    integer, intent(in) :: n_grid              !< number of grid points
    integer, intent(in) :: propflag            !< specified property 1:s, 2:lnv
    real, intent(inout) :: prop_grid(n_grid)   !< property grid (non sequential mode may reorder array)
    real, optional, intent(in) :: dS_override  !< step length along envelope
    integer, optional, intent(in) :: phase_in  !< start search from this phase
    integer, optional, intent(out) :: ierr_out !< error flag; nonzero if error
    real, optional, intent(in) :: sgn_in       !< Override search direction
    integer, optional, intent(in) :: spec_in   !< Override initial search direction
    real, optional, intent(in) :: tMin, tMax   !< Temperature limits (K)
    real, optional, intent(in) :: pMin, pMax   !< Pressure limits (Pa)
    logical, optional, intent(in) :: sequential_mode !< Default sequential
    ! Output:
    real, intent(out) :: T_grid(n_grid)        !< t at the grid points
    real, intent(out) :: P_grid(n_grid)        !< p at the grid points
    integer, intent(out) :: phase_grid(n_grid) !< incumbent phase at grid points
    real, intent(out) :: wi_grid(nc,n_grid)    !< Incipient phase at boundary
    integer, intent(out) :: n_grid_found       !< Number of grid points found
    real, optional, intent(out) :: normal_grid(2,n_grid) !< Normal vector (from tangent) pointing into two-phase region
    ! Locals:
    real :: t, p
    integer :: specification
    real, dimension(nc) :: K, lnfug, lnfugG, lnfugL
    real, dimension(nc+2) :: dXdS, dXdSold
    real, dimension(nc+2) :: Xvar, Xold, Xtemp
    real :: dS, tuning, sgn, Pstart
    integer :: iter,s,n,ophase,ierr,j
    integer, dimension(1) :: smax,zmax
    real :: curr_prop_val, old_prop_val, ln_spec, beta
    logical :: have_switched_formulation, seq_mode
    logical :: should_switch_formulation, reset_dS
    integer :: grid_idx, phase, start_phase, n_grid_not_found, n_found
    real :: sspec ! property value to bracket
    real :: dSmin, dSmax
    real :: prop_grid_local(n_grid), prop_grid_list(n_grid)      !< property grids
    logical :: exit_after_testing, recalculate
    real :: t_c,v_c,p_c,s_c,lnv_c,h_c,crit_prop

    ! Start envelope traversal from dew or bubble curve?
    if ( present(phase_in) ) then
      start_phase = phase_in
    else
      start_phase = VAPPH ! dew curve is default
    end if

    ! Locate critical point and set property at critical point
    call get_crit_prop(Z,t_c,v_c,p_c,s_c,h_c,lnv_c,propflag,crit_prop)

    ! Set mode of operation
    if (present(sequential_mode)) then
      seq_mode = sequential_mode
    else
      seq_mode = .true.
    endif
    if (.not. seq_mode) then
      ! Copy to local prop_grid
      prop_grid_local = prop_grid
      n_grid_not_found = n_grid
    endif

    ! Special consideration for pure fluids
    zmax = maxloc(Z)
    if (isSingleComp(Z)) then
      call sat_points_based_on_prop_single(Z,P0,n_grid,propflag,prop_grid,&
           T_grid,P_grid,phase_grid,n_grid_found,ierr_out)
      do j=1,n_grid_found
        wi_grid(:,j) = Z
      enddo
      return
    endif

    ! Initialize error flag. If something went wrong, set ierr_out /= 0.
    if (present(ierr_out)) ierr_out = 0

    ! Initialize variables needed to solve for initial point on envelope
    if (start_phase == VAPPH) then
      beta = 1.0
    else
      beta = 0.0
    end if
    t = t0
    p = p0
    call thermo(t0,p0,y0,VAPPH,lnfugG)
    call thermo(t0,p0,x0,LIQPH,lnfugL)
    K = exp(lnfugL-lnfugG)

    ! Solve for initial saturation point temperature and equilibrium factors K
    tuning = 1.2
    if (present(sgn_in)) then
      sgn = sgn_in
    else
      sgn = 1.0
    endif
    if (present(spec_in)) then
      specification = spec_in
    else
      specification = specP
    endif
    Xvar(1:nc) = log(K)
    Xvar(nc+1) = log(t)
    Xvar(nc+2) = log(p)
    s = ispec(specification) ! index of fixed variable in Xvar below (here nc+2)
    ln_spec = Xvar(s)
    iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
    if (ierr /= 0) then
      if (present(ierr_out)) then
        ierr_out = ierr
        return
      else
        call stoperror('sat grid: Initial point not found.')
      endif
    endif

    ! Retrieve solved variables
    Xvar(1:nc) = log(K)
    Xvar(nc+1) = log(t)
    Xvar(nc+2) = log(p)
    Pstart = p

    ! Compute property value of initial point on phase envelope
    ! lnv usually decreases
    call genericProperty(t,p,Z,start_phase,propflag,old_prop_val)
    if (.not. present(sgn_in) .and. seq_mode) then
      if ((start_phase == VAPPH .and. old_prop_val < prop_grid(1))) then
        print *, "Initial saturation point T, P0, phase:", t, p, start_phase
        print *, "Initial property value:               ", old_prop_val
        print *, "First property value in grid:         ", prop_grid(1)
        call stoperror("sat grid: Adjust P0 so grid_prop(1) is on dew curve.")
      else if ((start_phase == LIQPH .and. old_prop_val > prop_grid(1))) then
        print *, "Initial saturation point T and P0:", t, p
        print *, "Initial property value: \t", old_prop_val
        print *, "First property value in grid: \t", prop_grid(1)
        call stoperror("sat grid: Adjust P0 so grid_prop(1) is on bubble curve.")
      endif
    endif

    ! When traversing the phase envelope, we may eventually go from the dew
    ! curve to the bubble curve, forcing us to invert what is incumbent and
    ! incipient phase. These variables control this switch.
    have_switched_formulation = .false.
    should_switch_formulation = .false.

    ! Set steplength along phase envelope
    if (present(dS_override)) then
      dSmax = dS_override
      dSmin = min(0.25*dSmax,0.025)
    else
      dSmax = 0.15
      dSmin = 0.025
    endif
    dS = dSmin

    ! Traverse the phase envelope and bracket the property grid values.
    n = 1
    grid_idx = 1
    n_grid_found = 0
    sspec = prop_grid(grid_idx)
    dXdS = 0.0
    exit_after_testing = .false.
    recalculate = .false.
    do while (.true.)
      dXdSold = dXdS
      call newton_extrapolate(Z,Xvar,dXdS,beta,s,ln_spec)
      smax = maxloc(abs(dXdS))
      ! Rescaling the sensitivities
      if ((.not. smax(1) == s) .and. n > 1) then
        s = smax(1)
        ! Rescaling the sensitivities
        sgn = sign(1.0,Xvar(s) - Xold(s))
        dXdS = dXdS / dXdS(s)
        dXdSold = dXdSold / dXdSold(s)
      endif
      Xold = Xvar
      Xvar = Xold + dXdS*dS*sgn
      K = exp(Xvar(1:nc))
      t = exp(Xvar(nc+1))
      p = exp(Xvar(nc+2))
      ln_spec = Xvar(s)
      iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
      reset_dS = .false.
      if (ierr /= 0) then
        ! Something went wrong.
        ! Attempt to make the step shorter.
        Xvar = Xold + dXdS*dS*sgn*0.5
        K = exp(Xvar(1:nc))
        t = exp(Xvar(nc+1))
        p = exp(Xvar(nc+2))
        ln_spec = Xvar(s)
        iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
        reset_dS = .true.
      endif
      if (ierr /= 0) then
        ! Something went wrong.
        ! Attempt to make the step longer.
        Xvar = Xold + dXdS*dS*sgn*2.0
        K = exp(Xvar(1:nc))
        t = exp(Xvar(nc+1))
        p = exp(Xvar(nc+2))
        ln_spec = Xvar(s)
        iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
        reset_dS = .true.
      endif
      if (ierr /= 0 .and. s > nc) then
        ! Set K-value to try passing critical point
        smax = maxloc(abs(dXdS(1:nc)))
        sgn = sign(1.0,sgn*dXdS(s)*dXdS(smax(1)))
        s = smax(1)
        dXdS = dXdS / dXdS(s)
        if (abs(Xold(s) + dS*sgn) < 1.0e-4) then
          ! Adapt dS to overshoot critical point
          dS = 1.0e-4 + abs(Xold(s))
        endif
        ! Keep pressure and temperature constant
        dXdS(nc+1:nc+2) = 0.0
        Xvar = Xold + dXdS*dS*sgn
        K = exp(Xvar(1:nc))
        t = exp(Xvar(nc+1))
        p = exp(Xvar(nc+2))
        ln_spec = Xvar(s)
        iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
        reset_dS = .true.
      endif
      if (ierr /= 0) then
        if (present(ierr_out)) then
          ierr_out = ierr
          return
        else
          call stoperror("grid error: Neither decreasing nor &
               &increasing the step helped")
        endif
      endif

      if (reset_dS) then
        ! Reset dS
        dS = dSmin
      endif
      Xvar(1:nc) = log(K)
      Xvar(nc+1) = log(t)
      Xvar(nc+2) = log(p)
      n = n + 1

      !Exit at specified minimum temperature
      if (present(tMin)) then
        ! Is temperature decreasing? - And below Tmin?
        if (Xvar(nc+1) - Xold(nc+1) < 0.0 .and. T < Tmin) then
          s = nc+1
          ln_spec = log(Tmin)
          recalculate = .true.
          if (present(ierr_out)) then
            ierr_out = -1
          endif
        endif
      endif
      if (recalculate) then
        ! Extrapolate from previous point
        call newton_extrapolate(Z,Xold,dXdS,beta,s,ln_spec)
        dS = ln_spec - Xold(s) ! Sign included
        Xvar = Xold + dXdS*dS
        K = exp(Xvar(1:nc))
        t = exp(Xvar(nc+1))
        p = exp(Xvar(nc+2))
        iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
        Xvar(1:nc) = log(K)
        Xvar(nc+1) = log(t)
        Xvar(nc+2) = log(p)
        exit_after_testing = .true.
      endif

      if (beta > 0.5) then
        phase = VAPPH
      else
        phase = LIQPH
      endif

      call genericProperty(t,p,Z,phase,propflag,curr_prop_val)
      if (seq_mode) then
        ! Check to see if sspec is bracketed.
        if ( sspec >= curr_prop_val .and. sspec <= old_prop_val .OR.&
             sspec <= curr_prop_val .and. sspec >= old_prop_val) then
          call locate_sat_prop(propflag,Z,s,sspec,Xvar,Xold,dXds,beta,ierr)
          if (ierr_demands_return(ierr,ierr_out)) return
          call store_point(Z,Xvar,phase,beta,grid_idx,n_grid,&
               T_grid,P_grid,phase_grid,wi_grid,n_grid_found,&
               normal_grid)
          if ( grid_idx > n_grid ) then
            return
          else
            curr_prop_val = sspec
            sspec = prop_grid(grid_idx)
          endif
        endif
      else ! .not. seq_mode
        ! Check to see if some of the specifications are bracketed.
        if ( isAnyPropBracketed(curr_prop_val,old_prop_val,n_grid,prop_grid_local,&
             n_grid_not_found,prop_grid_list,n_found)) then
          do j=1,n_found
            sspec = prop_grid_list(j)
            call locate_saturation_prop_point(Z,phase,s,sspec,&
                 propflag,crit_prop,curr_prop_val,old_prop_val,&
                 Xvar,Xold,dXds,beta,t_c,p_c,Xtemp,ierr)
            if (ierr_demands_return(ierr,ierr_out)) return
            prop_grid(grid_idx) = sspec ! Reorder array
            call store_point(Z,Xtemp,phase,beta,grid_idx,n_grid,&
                 T_grid,P_grid,phase_grid,wi_grid,n_grid_found,&
                 normal_grid)
          enddo
          if ( n_grid_found == n_grid ) then
            return
          endif
        endif
      endif
      ! update old_prop_val and continue
      old_prop_val = curr_prop_val

      if (exit_after_testing) then
        return
      endif
      ! Abort if we are outside boundaries
      if (present(tMax)) then
        if (t > tMax) then
          if (present(ierr_out)) then
            ierr_out = -1
          endif
          return
        endif
      endif
      if (present(pMin)) then
        if (p < pMin) then
          if (present(ierr_out)) then
            ierr_out = -1
          endif
          return
        endif
      endif
      if (present(pMax)) then
        if (p > pMax) then
          if (present(ierr_out)) then
            ierr_out = -1
          endif
          return
        endif
      endif

      ! While the possibility to switch formulation may have arisen,
      ! it is not necessarily yet possible to do it error-free.
      ! If the flag should_switch_formulation has been raised, attempt
      ! to switch until it can be done error-free.
      if (should_switch_formulation) then
        call changeFormulation(beta,Xvar,Xold,K,s,ln_spec,sgn)
        iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
        ! If error-free, change phase now. If not, try again next time.
        if (ierr == 0) then
          Xvar(1:nc) = log(K)
          Xvar(nc+1) = log(t)
          Xvar(nc+2) = log(p)
          have_switched_formulation = .true.
          should_switch_formulation = .false.
        else
          call changeFormulation(beta,Xvar,Xold,K,s,ln_spec,sgn)
        endif
      endif

      if (.not. have_switched_formulation) then
        ! Check if we have entered the SINGLEPH-area.
        ! If so, we can change formulation (beta)
        call thermo(t,p,z,MINGIBBSPH,lnfug,ophase=ophase)
        if (ophase == SINGLEPH) then
          ! Wait until next iteration to actually switch formulation.
          should_switch_formulation = .true.
        endif
      endif

      ! Tune dS up or down based on how fast sat_newton converged
      if (iter < 3) then
        dS = dS * tuning
      else if (iter > 5) then
        dS = dS * (1.0 - tuning)
      endif
      dS = max(min(dS,dSmax),dSmin)
    enddo
  contains
    function ierr_demands_return(ierr,ierr_out) result(do_return)
      integer, intent(in) :: ierr
      integer, optional, intent(out) :: ierr_out
      !
      Logical :: do_return
      do_return = .false.
      if (ierr /= 0) then
        if (present(ierr_out)) then
          if (ierr == 100) then
            ierr_out = ierr
            do_return = .true.
          else
            ! Accept tolerance error:
            ierr_out = 0
            if (ierr /= 100) then
              print *,'sat_points_based_on_prop:: Tolerance not met for bracket solver. Continuing'
            endif
          endif
        else
          if (ierr == 100) then
            print *,'sat_points_based_on_prop:: No solution for bracket solver'
            call stoperror("grid error: Bracket solver failed!")
          else
            print *,"grid error: Bracket solver failed! ierr=",ierr
          endif
        endif
      endif
    end function ierr_demands_return
  end subroutine sat_points_based_on_prop

  !-----------------------------------------------------------------------------
  !> Locate property on saturation line
  !!
  !! \author Morten Hammer 2020-01
  !-----------------------------------------------------------------------------
  subroutine locate_saturation_prop_point(Z,phase,s,sspec,&
       propflag,crit_prop,curr_prop_val,old_prop_val,&
       Xvar,Xold,dXds_in,beta,T_c,p_c,Xtemp,ierr)
    ! Input
    real, intent(in) :: Z(nc)                  !< total composition
    integer, intent(in) :: propflag            !< specified property 1:s, 2:lnv
    integer, intent(in) :: phase   !<
    real, intent(in) :: sspec  !<
    integer, intent(in) :: s   !<
    real, intent(in) :: crit_prop,curr_prop_val,old_prop_val !<
    real, intent(in) :: beta  !<
    real, intent(in) :: Xvar(nc+2), Xold(nc+2), dXdS_in(nc+2)  !<
    real, intent(in) :: t_c  !<
    real, intent(in) :: p_c  !<
    ! Output
    real, intent(out) :: Xtemp(nc+2)
    integer, intent(out) :: ierr   !<
    ! Locals
    integer :: l, iter
    real :: dsc, Xc, iter_prop_val, temp_prop_val, ln_spec
    logical :: critIsCrossed, brackets_solution
    real, dimension(nc+2) :: XtempOld, dXds
    real :: t, p
    real, dimension(nc) :: K
    critIsCrossed = ((crit_prop-curr_prop_val)*(crit_prop-old_prop_val) < 0)
    !critIsCrossed .false.
    if (critIsCrossed) then
      if (s <= nc) then
        Xc = 0.0
      elseif (s == nc + 1) then
        Xc = log(T_c)
      else!if (s == nc + 2) then
        Xc = log(P_c)
      endif
      if ((crit_prop-curr_prop_val)*(crit_prop-sspec) > 0) then
        ! Located on curr_prop_val side
        Xtemp = Xvar
        temp_prop_val = curr_prop_val
      else
        ! Located on old_prop_val side
        Xtemp = Xold
        temp_prop_val = old_prop_val
      endif
      XtempOld = Xtemp
      ! Approach critical point
      brackets_solution = .false.
      do l=1,5
        call newton_extrapolate(Z,Xtemp,dXdS,beta,s,ln_spec)
        dsc = (Xc - Xtemp(s))*0.5**l
        Xtemp = Xtemp + dXdS*dSc
        K = exp(Xtemp(1:nc))
        t = exp(Xtemp(nc+1))
        p = exp(Xtemp(nc+2))
        ln_spec = Xtemp(s)
        iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
        if (ierr /= 0) then
          ! Reset and include critical point in search
          Xtemp = XtempOld
          exit
        endif
        Xtemp(1:nc) = log(K)
        Xtemp(nc+1) = log(t)
        Xtemp(nc+2) = log(p)
        call genericProperty(t,p,Z,phase,propflag,iter_prop_val)
        if ((temp_prop_val-sspec)*(iter_prop_val-sspec) < 0) then
          brackets_solution = .true.
          exit
        endif
      enddo
      if (.not. brackets_solution) then
        ! Fallback to search between Xold and Xvar
        XtempOld = Xold
        Xtemp = Xvar
      endif
    else
      Xtemp = Xvar
      XtempOld = Xold
      dXds = dXds_in
    endif
    call locate_sat_prop(propflag,Z,s,sspec,Xtemp,XtempOld,dXds,beta,ierr)
  end subroutine locate_saturation_prop_point

  !-----------------------------------------------------------------------------
  !> Locate critical point and set property at critical point
  !!
  !! \author Morten Hammer 2020-01
  !-----------------------------------------------------------------------------
  subroutine get_crit_prop(Z,t_c,v_c,p_c,s_c,h_c,lnv_c,propflag,crit_prop)
    use critical, only: calcCriticalTV, calcCritical
    use eosTV, only: enthalpyTV, entropyTV
    use eos, only: specificvolume
    ! Input
    real, intent(in) :: Z(nc)
    integer, intent(in) :: propflag
    ! Output
    real, intent(out) :: t_c,v_c,p_c,s_c,h_c,lnv_c,crit_prop
    ! Locals
    integer :: nenv, ierr
    integer, parameter :: n_max = 500
    real :: Ta(n_max), pa(n_max), Ki(n_max,nc)
    real :: betai(n_max)
    real :: crit(2)
    type(thermo_model), pointer :: act_mod_ptr

    act_mod_ptr => get_active_thermo_model()
    ! Locate critical point and set property at critical point
    if (act_mod_ptr%EoSlib /= TREND) then
      t_c = 0.0
      v_c = 0.0
      call calcCriticalTV(t_c,v_c,Z,ierr,tol=1.0e-7,p=p_c)
    else
      ierr = 0
      ! Get phase envelope, including critical point in T,p
      call envelopePlot(z,300.0,1.0e5,spec=1,beta_in=1.0,Pmax=2.0e7,nmax=n_max,&
           Ta=Ta,pa=pa,Ki=Ki,betai=betai,n=nenv,crit=crit)
      T_c = crit(1)
      p_c = crit(2)
      call specificvolume(t_c,p_c,Z,LIQPH,v_c)
    endif
    if (ierr == 0) then
      lnv_c = log(v_c)
      call entropyTV(t_c,v_c,z,s_c)
      call enthalpyTV(t_c,v_c,z,h_c)
    else
      ! Set un-physical values
      lnv_c = 0.0
      s_c = 1.0e20
      h_c = 1.0e20
      p_c = 0.0
      t_c = 0.0
    endif
    select case(propflag)
    case(locate_from_entropy)
      crit_prop = s_c
    case(locate_from_lnvol)
      crit_prop = lnv_c
    case(locate_from_enthalpy)
      crit_prop = h_c
    case(locate_from_temperature)
      crit_prop = t_c
    case(locate_from_pressure)
      crit_prop = p_c
    end select
  end subroutine get_crit_prop

  !-----------------------------------------------------------------------------
  !> Locate the points on pure fluid saturation line corresponding to
  !! the property values in prop_grid.
  !!
  !! \author Morten 2019-05
  !-----------------------------------------------------------------------------
  subroutine sat_points_based_on_prop_single(Z,P0,n_grid,propflag,prop_grid,&
       T_grid,P_grid,phase_grid,n_grid_found,ierr_out)
    use eos, only: getCriticalParam
    implicit none
    ! Input:
    real, intent(in) :: Z(nc)                  !< total composition
    real, intent(in) :: p0                     !< init. point has pressure p0
    integer, intent(in) :: n_grid              !< number of grid points
    integer, intent(in) :: propflag            !< specified property 1:s, 2:lnv
    real, intent(inout) :: prop_grid(n_grid)   !< property grid (non sequential mode may reorder array)
    integer, optional, intent(out) :: ierr_out !< error flag; nonzero if error
    ! Output:
    real, intent(out) :: T_grid(n_grid)        !< t at the grid points
    real, intent(out) :: P_grid(n_grid)        !< p at the grid points
    integer, intent(out) :: phase_grid(n_grid) !< incumbent phase at grid points
    integer, intent(out) :: n_grid_found       !< Number of grid points found
    ! Locals:
    real :: t, p, Ps, Ts
    integer, dimension(1) :: zmax
    real, dimension(nc) :: X,Y
    real, dimension(4) :: param
    real, dimension(1) :: XX,Xmax,Xmin
    type(nonlinear_solver) :: solver
    integer :: ic, iter, ierr, j
    real :: f,dfdt,dfdp,tci,pci,oi,dpdt,dP
    integer , parameter :: nmax = 200
    real :: Ta(nmax*2), Pa(nmax*2), propa(nmax*2)
    integer :: phasea(nmax*2)
    real :: curr_prop_val, old_prop_val, crit_prop_val
    integer :: grid_idx, n_grid_not_found, n_found
    real :: sspec ! property value to bracket
    real :: prop_grid_local(n_grid), prop_grid_list(n_grid)      !< property grids

    ! Copy to local prop_grid
    prop_grid_local = prop_grid
    n_grid_not_found = n_grid

    ! Initialize error flag. If something went wrong, set ierr_out /= 0.
    if (present(ierr_out)) ierr_out = 0

    ! Initialize composition
    X = 0.0
    zmax = maxloc(Z)
    ic = zmax(1)
    X(ic) = 1.0
    Y=X

    ! Initialize variables needed to solve for initial point on envelope
    t = safe_bubT(p0,X,Y,ierr_out)

    call getCriticalParam(ic,tci,pci,oi)
    ! Find initlal point
    p = p0
    dP = (pci-P)/nmax
    xmax = log(tci)
    xmin = log(t)

    ! Initial point
    Ta(1) = t
    Pa(1) = p
    phasea(1) = VAPPH
    call genericProperty(t,p,Z,VAPPH,propflag,old_prop_val)
    propa(1) = old_prop_val
    Ta(nmax*2) = t
    Pa(nmax*2) = p
    phasea(nmax*2) = LIQPH
    call genericProperty(t,p,Z,LIQPH,propflag,curr_prop_val)
    propa(nmax*2) = curr_prop_val

    ! Critical point
    call genericProperty(tci,pci,X,LIQPH,propflag,crit_prop_val)
    Ta(nmax) = tci
    Pa(nmax) = pci
    phasea(nmax) = VAPPH
    propa(nmax) = crit_prop_val
    Ta(nmax+1) = tci
    Pa(nmax+1) = pci
    phasea(nmax+1) = LIQPH
    propa(nmax+1) = crit_prop_val

    ! Solver paramaters
    solver%rel_tol = 1.0e-20
    solver%abs_tol = 1.0e-10
    solver%max_it = 10000
    param(1) = specP
    param(2) = ic
    param(4) = 0.0
    do iter = 2,nmax-1
      call sat_fun_single(ic,t,p,f,dfdt,dfdp,.false.)
      dpdt = -dfdt/dfdp
      param(3) = p0 + dP*(iter-1)
      XX(1) = log(t+dP/dpdt)
      call nonlinear_solve(solver,sat_fun_single_if,sat_diff_single,&
           sat_diff_single,limit_dx,premReturn,setXv,XX,xmin,xmax,param)
      if (solver%exitflag == 0) then
        T = exp(XX(1))
        P = param(3)
      else
        if (present(ierr_out)) then
          ierr_out = solver%exitflag
          return
        else
          call stoperror("sat_points_based_on_prop_single: Not able to locate point on saturation line")
        endif
      endif
      Ta(iter) = t
      Pa(iter) = p
      phasea(iter) = VAPPH
      call genericProperty(t,p,X,VAPPH,propflag,curr_prop_val)
      propa(iter) = curr_prop_val
      Ta(2*nmax-iter+1) = t
      Pa(2*nmax-iter+1) = p
      phasea(2*nmax-iter+1) = LIQPH
      call genericProperty(t,p,X,LIQPH,propflag,curr_prop_val)
      propa(2*nmax-iter+1) = curr_prop_val
    enddo

    ! Traverse the saturation line and bracket the property grid values.
    grid_idx = 1
    n_grid_found = 0
    old_prop_val = propa(1)
    do iter = 2,2*nmax
      curr_prop_val = propa(iter)
      ! Check to see if some of the specifications are bracketed.
      if ( isAnyPropBracketed(curr_prop_val,old_prop_val,n_grid,prop_grid_local,&
           n_grid_not_found,prop_grid_list,n_found)) then
        do j=1,n_found
          sspec = prop_grid_list(j)
          Ts=Ta(iter)
          Ps=Pa(iter)
          call locate_sat_prop_single(propflag,ic,sspec,&
               Ta(iter-1),Pa(iter-1),phasea(iter-1),Ts,Ps,ierr)
          !print *,iter,ierr,Ts,Ps,phasea(iter-1)
          if (ierr /= 0) then
            if (present(ierr_out)) then
              ierr_out = ierr
              return
            else
              if (ierr == 100) then
                print *,'No solution for bracket solver'
              else
                print *,'ierr',ierr
              endif
              call stoperror("grid error: Bracket solver failed!")
            endif
          endif
          prop_grid(grid_idx) = sspec ! Reorder array
          call store_point_single(Z,Ts,Ps,phasea(iter),grid_idx,n_grid,&
               T_grid,P_grid,phase_grid,n_grid_found)
          if ( n_grid_found == n_grid ) then
            return
          end if
        enddo
      endif
      old_prop_val = curr_prop_val
    enddo
  end subroutine sat_points_based_on_prop_single

  !-----------------------------------------------------------------------------
  !> Locate property - envelope intersect
  !>
  !> \author MH, 2018-06
  !-----------------------------------------------------------------------------
  subroutine locate_sat_prop(propflag,Z,s,sspec,Xvar,Xold,dXds,beta,ierr)
    implicit none
    integer, intent(in) :: propflag
    integer, intent(in) :: s
    real, intent(in) :: Z(nc)
    real, intent(in) :: sspec
    real, intent(in) :: beta
    real, intent(inout) :: Xvar(nc+2),Xold(nc+2),dXds(nc+2)
    integer, intent(out) :: ierr
    ! Locals
    integer :: sl, iter
    real :: K(nc), ln_spec, p, t
    if ( propflag == locate_from_temperature .OR. &
         propflag == locate_from_pressure) then
      if (propflag == locate_from_temperature) then
        sl = nc+1
      else
        sl = nc+2
      endif
      ln_spec = log(sspec)
      Xvar = Xold + dXdS*(ln_spec - Xold(sl))
      K = exp(Xvar(1:nc))
      t = exp(Xvar(nc+1))
      p = exp(Xvar(nc+2))
      iter = sat_newton(Z,K,t,p,beta,sl,ln_spec,ierr)
    else
      ! sspec is bracketed, so solve for exact property value and return
      call bracketSolveForProperty(n=nc+2,Z=Z,beta=beta,propflag=propflag,&
           propspec=sspec,X=Xvar,Xold=Xold,dXdS=dXdS,s=s,ierr=ierr)
      ! Restart at located point
      !dS = dSmin
    endif
  end subroutine locate_sat_prop

  !-----------------------------------------------------------------------------
  !> Locate property - saturation line intersect
  !>
  !> \author MH, 2019-05
  !-----------------------------------------------------------------------------
  subroutine locate_sat_prop_single(propflag,ic,sspec,T1,P1,phase,Ts,Ps,ierr)
    implicit none
    integer, intent(in) :: propflag
    integer, intent(in) :: ic, phase
    real, intent(in) :: T1,P1
    real, intent(in) :: sspec
    real, intent(inout) :: Ts,Ps
    integer, intent(out) :: ierr
    ! Locals
    real :: f,dfdt,dfdp, p, t, dpdt, dT, dP
    real :: XX(1), param(4), xmax(1), xmin(1)
    type(nonlinear_solver) :: solver
    t = T1
    p = p1
    call sat_fun_single(ic,t,p,f,dfdt,dfdp,.false.)
    dpdt = -dfdt/dfdp
    if ( propflag == locate_from_temperature .OR. &
         propflag == locate_from_pressure) then
      if (propflag == locate_from_temperature) then
        param(1) = specT
        param(3) = sspec
        dT = sspec - t
        xmin(1) = log(min(P1,Ps))
        xmax(1) = log(max(P1,Ps))
        XX(1) = log(p+dT*dpdt)
      else
        param(1) = specP
        param(3) = sspec
        dP = sspec - p
        xmin(1) = log(min(T1,Ts))
        xmax(1) = log(max(T1,Ts))
        XX(1) = log(t+dP/dpdt)
      endif
      XX(1) = min(max(XX(1),xmin(1)),xmax(1))
      param(2) = ic
      param(4) = 0.0
      ! Solver paramaters
      solver%rel_tol = 1.0e-20
      solver%abs_tol = 1.0e-10
      solver%max_it = 100
      call nonlinear_solve(solver,sat_fun_single_if,sat_diff_single,&
           sat_diff_single,limit_dx,premReturn,setXv,XX,xmin,xmax,param)
      ierr = solver%exitflag
      if (ierr == 0) then
        if (propflag == locate_from_temperature) then
          Ts = sspec
          Ps = exp(XX(1))
        else
          Ps = sspec
          Ts = exp(XX(1))
        endif
      endif
    else
      ! sspec is bracketed, so solve for exact property value and return
      call bracketSolveForPropertySingle(ic,propflag,sspec,phase,&
           T1,P1,dPdT,Ts,Ps,ierr)
    endif
  end subroutine locate_sat_prop_single

  !-----------------------------------------------------------------------------
  !> Store point in output vector
  !>
  !> \author MH, 2018-06
  !-----------------------------------------------------------------------------
  subroutine store_point(Z,Xvar,phase,beta,grid_idx,n_grid,&
       T_grid,P_grid,phase_grid,wi_grid,n_grid_found,&
       normal_grid)
    use solid_saturation, only: sat_diff_newton_threePh
    implicit none
    ! Input:
    real, dimension(nc), intent(in) :: Z       !< Overall composition
    real, dimension(nc+2) :: Xvar              !< Variable vector
    real, intent(in) :: beta                   !< Two-phase gas phase fraction (mol/mol)
    integer, intent(in) :: phase               !< Single phase at envelope crossing
    integer, intent(inout) :: grid_idx         !< Current index, increased by one before exit
    integer, intent(in) :: n_grid              !< number of grid points
    ! Output:
    real, intent(out) :: T_grid(n_grid)        !< t at the grid points
    real, intent(out) :: P_grid(n_grid)        !< p at the grid points
    integer, intent(out) :: phase_grid(n_grid) !< incumbent phase at grid points
    real, intent(out) :: wi_grid(nc,n_grid)    !< Incipient phase at boundary
    integer, intent(out) :: n_grid_found       !< Number of grid points found
    real, optional, intent(out) :: normal_grid(2,n_grid) !< Normal vector (from tangent) pointing into two-phase region
    ! Locals:
    real :: K(nc), tangent(2), normal(2), dlnTdX(2), dlnPdX(2), dbeta(2)
    real, dimension(nc+3) :: Xvar3ph !< Variable vector
    real, dimension(nc+3,nc+3) :: Jac3ph !< Function differentials
    real, dimension(nc+5) :: param3ph !< Parameter vector
    real, dimension(nc+2,2) :: rhs !< Variable vector
    real, dimension(nc+2,nc+2) :: Jac2ph !< Function differentials
    !real, dimension(nc+5) :: param3ph !< Parameter vector
    integer, dimension(nc+2) :: INDX
    integer :: INFO

    K = exp(Xvar(1:nc))
    T_grid(grid_idx) = exp(Xvar(nc+1))
    P_grid(grid_idx) = exp(Xvar(nc+2))
    phase_grid(grid_idx) = phase
    if (beta > 0.5) then
      wi_grid(:,grid_idx) = Z/(1-beta+beta*K)
    else
      wi_grid(:,grid_idx) = K*Z/(1-beta+beta*K)
    endif
    if (present(normal_grid)) then
      Xvar3ph(1:nc+2) = Xvar
      Xvar3ph(nc+3) = beta
      param3ph(1:nc) = z
      param3ph(nc+1) = 0.0  ! Solid fraction
      param3ph(nc+2) = nc+2 ! Const P
      param3ph(nc+3) = 0.0  ! Dummy spec.
      param3ph(nc+4) = real(VARFLUID_DUMMY)
      param3ph(nc+5) = 1.0 ! Dummy solid component
      call sat_diff_newton_threePh(Jac3ph,Xvar3ph,param3ph)
      Jac2ph = Jac3ph(1:nc+2,1:nc+2)
      rhs(:,1) = -Jac3ph(1:nc+2,nc+3)
      rhs(:,2) = 0.0
      rhs(nc+2,2) = 1.0 ! Perturbation in pressure
      ! Solve equation system
      call DGESV( nc+2, 2, Jac2ph, nc+2, INDX, rhs, nc+2, INFO )
      dlnTdX(1) = rhs(nc+1,2) ! dlnT_dlnP
      dlnTdX(2) = rhs(nc+1,1) ! dlnT_dbeta
      !
      Jac2ph = -Jac3ph(1:nc+2,1:nc+2)
      Jac2ph(nc+2,nc+2) = 0.0
      Jac2ph(nc+2,nc+1) = 1.0 ! Fixate temperature
      rhs(:,1) = Jac3ph(1:nc+2,nc+3)
      rhs(:,2) = 0.0
      rhs(nc+2,2) = 1.0 ! Perturbation in temperature
      ! Solve equation system
      call DGESV( nc+2, 2, Jac2ph, nc+2, INDX, rhs, nc+2, INFO )
      dlnPdX(1) = rhs(nc+2,2) ! dlnP_dlnT
      dlnPdX(2) = rhs(nc+2,1) ! dlnP_dbeta
      !
      if (dlnPdX(1) > dlnTdX(1)) then
        tangent = (/dlnTdX(1), 1.0/)
      else
        tangent = (/1.0, dlnPdX(1)/)
      endif
      normal(1) = - tangent(2)
      normal(2) = tangent(1)
      dbeta(1) = 1.0/dlnTdX(2)
      dbeta(2) = 1.0/dlnPdX(2)
      if (beta > 0.5 .and. sum(normal*dbeta) > 0.0) then
        normal = - normal
      else if (beta < 0.5 .and. sum(normal*dbeta) < 0.0) then
        normal = - normal
      endif
      normal = normal/sqrt(dot_product(normal,normal))
      normal_grid(:,grid_idx) = normal
    endif
    n_grid_found = grid_idx
    grid_idx = grid_idx + 1
  end subroutine store_point

  !-----------------------------------------------------------------------------
  !> Store point in output vector for pure fluid
  !>
  !> \author MH, 2019-05
  !-----------------------------------------------------------------------------
  subroutine store_point_single(Z,T,P,phase,grid_idx,n_grid,&
       T_grid,P_grid,phase_grid,n_grid_found)
    implicit none
    ! Input:
    real, dimension(nc), intent(in) :: Z       !< Overall composition
    real, intent(in)                :: T,P     !<
    integer, intent(in) :: phase               !< Single phase at envelope crossing
    integer, intent(inout) :: grid_idx         !< Current index, increased by one before exit
    integer, intent(in) :: n_grid              !< number of grid points
    ! Output:
    real, intent(out) :: T_grid(n_grid)        !< t at the grid points
    real, intent(out) :: P_grid(n_grid)        !< p at the grid points
    integer, intent(out) :: phase_grid(n_grid) !< incumbent phase at grid points
    integer, intent(out) :: n_grid_found       !< Number of grid points found
    ! Locals:

    T_grid(grid_idx) = T
    P_grid(grid_idx) = P
    phase_grid(grid_idx) = phase
    n_grid_found = grid_idx
    grid_idx = grid_idx + 1
  end subroutine store_point_single

  !-----------------------------------------------------------------------------
  !> Test numdiff beta
  !>
  !> \author MH, 2018-06
  !-----------------------------------------------------------------------------
  subroutine num_diff_beta(Z,Xvar,beta,dbeta)
    implicit none
    ! Input:
    real, dimension(nc), intent(in) :: Z       !< Overall composition
    real, dimension(nc+2) :: Xvar              !< Variable vector
    real, intent(in) :: beta                   !< Two-phase gas phase fraction (mol/mol)
    ! Output:
    real, intent(out) :: dbeta(2)              !< dbeta
    ! Locals:
    real :: K(nc), t0, p0, t, p, beta_eps, db, ln_spec
    real, parameter :: eps = 1.0e-5
    integer :: ierr, iter, s
    K = exp(Xvar(1:nc))
    t0 = exp(Xvar(nc+1))
    p0 = exp(Xvar(nc+2))
    ln_spec = Xvar(nc+2)
    s = nc+2
    iter = sat_newton(Z,K,t0,p0,beta,s,ln_spec,ierr)
    t = t0
    p = p0
    if (beta < 0.5) then
      db = eps
    else
      db = - eps
    endif
    beta_eps = beta + db
    iter = sat_newton(Z,K,t,p,beta_eps,s,ln_spec,ierr)
    dbeta(1) = db/(log(t) - log(t0))
    ln_spec = Xvar(nc+1)
    s = nc+1
    iter = sat_newton(Z,K,t,p,beta_eps,s,ln_spec,ierr)
    dbeta(2) = db/(log(p) - log(p0))
  end subroutine num_diff_beta

  !-----------------------------------------------------------------------------
  !> TNumerical dlntdlnp
  !>
  !> \author MH, 2018-06
  !-----------------------------------------------------------------------------
  subroutine num_diff_dlntdlnp(Z,Xvar,beta,dlntdlnp)
    implicit none
    ! Input:
    real, dimension(nc), intent(in) :: Z       !< Overall composition
    real, dimension(nc+2) :: Xvar              !< Variable vector
    real, intent(in) :: beta                   !< Two-phase gas phase fraction (mol/mol)
    ! Output:
    real, intent(out) :: dlntdlnp              !< dlntdlnp
    ! Locals:
    real :: K(nc), t0, p0, t, p, ln_spec
    real, parameter :: eps = 1.0e-5
    integer :: ierr, iter, s
    K = exp(Xvar(1:nc))
    t0 = exp(Xvar(nc+1))
    p0 = exp(Xvar(nc+2))
    ln_spec = Xvar(nc+2)
    s = nc+2
    iter = sat_newton(Z,K,t0,p0,beta,s,ln_spec,ierr)
    p = p0 + p0*eps
    ln_spec = log(p)
    t = t0
    iter = sat_newton(Z,K,t,p,beta,s,ln_spec,ierr)
    dlntdlnp = (log(t) - log(t0))/(log(p) - log(p0))
  end subroutine num_diff_dlntdlnp

  !-----------------------------------------------------------------------------
  !> Iso-line crossing of saturation line
  !>
  !> \author MH, 2017-06
  !-----------------------------------------------------------------------------
  subroutine iso_cross_saturation_line(ts,ps,t,p,X,Y,Z,beta,&
       prop_spec,prop_specID,phase,ierr)
    use stringmod, only: str_eq
    use eos, only: enthalpy, entropy
    use saturation_curve, only: extrapolate_to_saturation_line
    implicit none
    real, dimension(nc), intent(in) :: Z
    real, intent(in) :: ts !< Single phase temperature (K)
    real, intent(in) :: ps !< Single phase pressure (Pa)
    real, intent(inout) :: t !< Two-phase temperature (K)
    real, intent(inout) :: p !< Two-phase pressure (Pa)
    real, intent(inout) :: beta !< Two-phase gas phase fraction (mol/mol)
    real, dimension(nc), intent(in) :: X !< Liquid composition
    real, dimension(nc), intent(in) :: Y !< Gas composition
    real, intent(in) :: prop_spec ! Enthalpy, entropy or temperature
    character, intent(in) :: prop_specID !< "h", "s", "t", "p"
    integer, intent(out) :: phase !< Single phase at envelope crossing
    integer, intent(out) :: ierr !< Error indicator, zero on success
    ! Locals
    real :: dpds_sat, dtds_sat, h, dhdt, dhdp, s, dsdt, dsdp
    real :: dtdp, dpdt, dh, ds, dT, dP, sgn
    integer :: extrap ! 1 = isobaric, 2 = isothermal
    integer :: spec, propflag
    integer, parameter :: n_grid = 1
    real :: T_grid(n_grid), P_grid(n_grid), prop_grid(n_grid)
    real :: wi_grid(nc,n_grid)
    integer :: phase_grid(n_grid), n_grid_found
    real, dimension(nc) :: X_cpy, Y_cpy
    real :: tMax, tMin, pMax, pMin
    ! How to extrapolate?
    if (str_eq(prop_specID,"T")) then
      extrap = ISO_T ! isothermal
    else if (str_eq(prop_specID,"P")) then
      extrap = ISO_P ! isobaric
    else if (str_eq(prop_specID,"S")) then
      extrap = ISO_S ! isentropic
    else if (str_eq(prop_specID,"H")) then
      extrap = ISO_H ! isenthalpic
    else
      print *,"Wrong extrapolation variable specified"
      return
    endif
    X_cpy = X
    Y_cpy = Y
    call extrapolate_to_saturation_line(t,p,X_cpy,Y_cpy,Z,beta,&
         extrap,ierr,dtds_sat,dpds_sat)
    if (ierr == 0) then
      if (abs(dtds_sat)/t > abs(dpds_sat)/p) then
        spec = specT
        dpdt = dpds_sat/dtds_sat
      else
        spec = specP
        dtdp = dtds_sat/dpds_sat
      endif
      if (beta > 0.5) then
        phase = VAPPH
      else
        phase = LIQPH
      endif
      if (prop_specID == "h" .or. prop_specID == "H") then
        propflag = locate_from_enthalpy
        call enthalpy(T,P,Z,phase,h,dhdt,dhdp)
        if (spec == specT) then
          dh = dhdt + dhdp*dpdt
          dT = (prop_spec-h)/dh
          sgn = sign(1.0,dT)
        else ! spec == specP
          dh = dhdt*dtdp + dhdp
          dP = (prop_spec-h)/dh
          sgn = sign(1.0,dP)
        endif
      else if (prop_specID == "s" .or. prop_specID == "S") then
        propflag = locate_from_entropy
        call entropy(T,P,Z,phase,s,dsdt,dsdp)
        if (spec == specT) then
          ds = dsdt + dsdp*dpdt
          dT = (prop_spec-s)/ds
          sgn = sign(1.0,dT)
        else ! spec == specP
          ds = dsdt*dtdp + dsdp
          dP = (prop_spec-s)/ds
          sgn = sign(1.0,dP)
        endif
      else
        ! Isotherm and isobar
        if ( (extrap == 2 .and. abs(prop_spec - T)/prop_spec < machine_prec*1.0e7) .or.&
             (extrap == 1 .and. abs(prop_spec - P)/prop_spec < machine_prec*1.0e7)) then
          ! Converged
          return
        else
          ! Need to search for T or P
          if (extrap == 2) then
            propflag = locate_from_temperature
            if (spec == specT) then
              sgn = sign(1.0, prop_spec - T)
            else !spec == specP
              sgn = sign(1.0, (prop_spec - T)*dpdt)
            endif
          else !extrap == 1
            propflag = locate_from_pressure
            if (spec == specP) then
              sgn = sign(1.0, prop_spec - P)
            else !spec == specT
              sgn = sign(1.0, (prop_spec - P)*dtdp)
            endif
          endif
        endif
      endif
      ! Search for property on envelope
      prop_grid = prop_spec
      tMax = ts + 5.0
      tMin = ts - 5.0
      pMax = ps + 5.0e5
      pMin = max(ps - 5.0e5, 1.0e5)
      call sat_points_based_on_prop(Z,T,P,X_cpy,Y_cpy,n_grid,&
           propflag,prop_grid,T_grid,P_grid,phase_grid,&
           wi_grid,n_grid_found,&
           0.05,phase,ierr,sgn,spec,tMin=tMin,tMax=tMax,&
           pMin=pMin,pMax=Pmax)
      T = T_grid(1)
      P = P_grid(1)
      phase = phase_grid(1)
    endif
  end subroutine iso_cross_saturation_line

  !-----------------------------------------------------------------------------
  !> Locate the points on the sublimation line corresponding to the property values
  !! in prop_grid.
  !!
  !! The routine starts at the "triple" point, and then
  !! traverses the sublimation line until it has bracketed one or more grid value; it then
  !! locates the exact sublimation point corresponding to this value, before
  !! continuing.
  !!
  !! \author Morten 2018-06
  !-----------------------------------------------------------------------------
  subroutine subl_points_based_on_prop(Z,Ttr,Ptr,Pmin,n_grid,propflag,prop_grid,&
       T_grid,P_grid,n_grid_found,ierr_out)
    use solideos, only: solidComp, nSolid
    implicit none
    ! Input:
    real, intent(in) :: Z(nc)                  !< total composition
    real, intent(in) :: ttr, ptr               !< triple point
    integer, intent(in) :: n_grid              !< number of grid points
    integer, intent(in) :: propflag            !< specified property 1:s, 2:lnv
    real, intent(inout) :: prop_grid(n_grid)   !< property grid
    real, intent(in) :: pMin                   !< Pressure limit (Pa)
    ! Output:
    real, intent(out) :: T_grid(n_grid)        !< t at the grid points
    real, intent(out) :: P_grid(n_grid)        !< p at the grid points
    integer, intent(out) :: n_grid_found       !< gird points found
    integer, optional, intent(out) :: ierr_out !< error flag; nonzero if error
    ! Locals:
    real :: t, p
    integer :: i,iter,is,ierr,s,phase,j
    real :: curr_prop_val, old_prop_val, ln_spec, betaSol, dln_s, pn, sspec
    real, dimension(2) :: Xvar, Xold, dXds, Xtemp
    integer :: n_found, nsol, n_grid_not_found
    real, parameter :: dPmax = 0.2e5
    real :: prop_grid_local(n_grid), prop_grid_list(n_grid)      !< property grids
    ! Mapping solid-gas line starting at triple point
    if (nSolid /= 1) then
      call stoperror('saturation::solidEnvelopePlot: No solid model initialized!')
    else
      is = solidComp(1)
    endif
    ! Initial state
    t = Ttr
    p = Ptr
    phase = VAPPH

    ! Compute property value of initial point on sublimation line
    call genericProperty(t,p,Z,phase,propflag,old_prop_val)
    prop_grid_local = prop_grid
    n_grid_found = 0
    n_grid_not_found = n_grid

    betaSol = 0.0
    s = 2 ! Specify pressure
    nsol = max(int((Ptr-Pmin)/dPmax)+1,50)
    do i=1,nsol
      Xold = (/ log(t), log(p)/)
      pn = Ptr - (Ptr-Pmin)*real(i-1)/real(nsol-1)
      ln_spec = log(pn)
      dln_s = ln_spec - log(p)
      call newton_extrapolate_solid(Z,t,p,betaSol,VAPPH,s,is,ln_spec,dln_s,dXds)
      iter = solidPointOnEnvelope(Z,t,p,phase,betaSol,is,s,ln_spec,ierr)
      if (ierr /= 0) then
        if (present(ierr_out)) then
          ierr_out = ierr
          return
        else
          call stoperror("grid error: Not able to solve for sublimation line point")
        endif
      endif
      Xvar = (/ log(t), log(p)/)
      call genericProperty(t,p,Z,phase,propflag,curr_prop_val)

      ! Check to see if sspec is bracketed.
      if ( isAnyPropBracketed(curr_prop_val,old_prop_val,n_grid,prop_grid_local,&
           n_grid_not_found,prop_grid_list,n_found)) then
        do j=1,n_found
          sspec = prop_grid_list(j)
          if ( propflag == locate_from_temperature .OR. &
               propflag == locate_from_pressure) then
            if (propflag == locate_from_temperature) then
              s = 1
            else
              s = 2
            endif
            ln_spec = log(sspec)
            Xtemp = Xold + dXdS*(ln_spec - Xold(s))
            t = exp(Xtemp(1))
            p = exp(Xtemp(2))
            iter = solidPointOnEnvelope(Z,t,p,phase,betaSol,is,s,ln_spec,ierr)
          else
            Xtemp = Xvar
            ! sspec is bracketed, so solve for exact property value and return
            call bracketSolveForProperty(n=2,Z=Z,beta=1.0,propflag=propflag,&
                 propspec=sspec,X=Xtemp,Xold=Xold,dXdS=dXdS,s=s,ierr=ierr)
            t = exp(Xtemp(1))
            p = exp(Xtemp(2))
          endif
          if (ierr /= 0) then
            if (present(ierr_out)) then
              ierr_out = ierr
              return
            else
              if (ierr == 100) then
                print *,'No solution for bracket solver'
              else
                print *,'ierr',ierr
              endif
              call stoperror("grid error: Bracket solver failed!")
            endif
          endif

          call store_subl_point(Z,t,p,sspec,n_grid,T_grid,P_grid,n_grid_found,prop_grid)
          if ( n_grid_found == n_grid ) then
            return
          end if
        enddo
      else
        ! sspec not bracketed, so update old_prop_val and continue
        old_prop_val = curr_prop_val
      endif

    enddo

  end subroutine subl_points_based_on_prop

  !-----------------------------------------------------------------------------
  !> Store point in output vector
  !>
  !> \author MH, 2018-06
  !-----------------------------------------------------------------------------
  subroutine store_subl_point(Z,t,p,prop,n_grid,T_grid,P_grid,n_grid_found,prop_grid)
    implicit none
    ! Input:
    real, dimension(nc), intent(in) :: Z       !< Overall composition
    real, intent(in) :: t,p                    !< Two-phase gas phase fraction (mol/mol)
    integer, intent(in) :: n_grid              !< number of grid points
    real, intent(in) :: prop                   !< Proporty at point
    ! Output:
    real, intent(inout) :: T_grid(n_grid)        !< t at the grid points
    real, intent(inout) :: P_grid(n_grid)        !< p at the grid points
    integer, intent(inout) :: n_grid_found       !< Number of grid points found
    real, intent(inout) :: prop_grid(n_grid)     !< property grid
    ! Locals:
    n_grid_found = n_grid_found + 1
    T_grid(n_grid_found) = t
    P_grid(n_grid_found) = p
    prop_grid(n_grid_found) = prop
  end subroutine store_subl_point

  !-----------------------------------------------------------------------------
  !> Function used to solve for point on sublimation line having property value
  !> propspec.
  !>
  !> \author MH 2018-6
  !-----------------------------------------------------------------------------
  function propertyFunctionWrapperSubl(Xs,param) result(fun)
    use solideos, only: solidComp
    implicit none
    real, intent(inout) :: Xs
    real, dimension(nc+10), intent(inout) :: param
    real :: fun
    ! Locals:
    integer :: iter, s, phase, is, ierr
    real, dimension(nc) :: Z
    real, dimension(2) :: X,Xold,dXdS
    real :: t, prop, propspec, ln_spec, p, beta, betaSol
    integer :: propflag

    ! Unpack the param vector
    propspec = param(1)
    beta = param(2)
    s = int(param(3))
    Z = param(4:nc+3)
    Xold = param(nc+4:nc+5)
    dXdS = param(nc+6:nc+7) ! is constant each time function is called
    X = param(nc+8:nc+9)
    propflag = int(param(nc+10))

    ! Extrapolate for better initial values in sat_newton
    X = Xold + dXdS*(Xs-Xold(s))
    t = exp(X(1))
    p = exp(X(2))

    ! Get correct phase
    if (beta > 0.5) then
      phase = VAPPH
    else
      phase = LIQPH
    endif
    betaSol = 0.0
    is = solidComp(1)

    ! Find point on envelope corresponding to ln(property value) equal to Xs
    ln_spec = Xs
    iter = solidPointOnEnvelope(Z,t,p,phase,betaSol,is,s,ln_spec,ierr)
    X(1) = log(T)
    X(2) = log(p)

    ! Update Xold-slot in param vector with new X value
    param(nc+8:nc+9) = X

    ! Calculate the property value at the new saturation point
    call genericProperty(t,p,Z,phase,propflag,prop)

    ! Compute new value of objective function
    fun = (prop - propspec)/max(abs(propspec), 1.0)
  end function propertyFunctionWrapperSubl

  !-----------------------------------------------------------------------------
  !> Function used to solve for point on sublimation line having property value
  !> propspec.
  !>
  !> \author MH 2018-6
  !-----------------------------------------------------------------------------
  function isAnyPropBracketed(curr_prop_val,old_prop_val,n_grid,prop_grid,&
       n_grid_not_found,prop_list,n_found) result(isBracketed)
    implicit none
    real, intent(in) :: curr_prop_val,old_prop_val
    integer, intent(in) :: n_grid
    real, dimension(n_grid), intent(inout) :: prop_grid
    integer, intent(inout) :: n_grid_not_found
    real, dimension(n_grid), intent(out) :: prop_list
    integer, intent(out) :: n_found
    logical :: isBracketed
    ! Locals:
    integer :: i,j,k
    isBracketed = .false.
    n_found = 0
    do i=1,n_grid_not_found
      ! Check to see if specification is bracketed.
      if ( prop_grid(i) >= curr_prop_val .and. prop_grid(i) <= old_prop_val .OR. &
           prop_grid(i) <= curr_prop_val .and. prop_grid(i) >= old_prop_val) then
        isBracketed = .true.
        n_found = n_found + 1
        prop_list(n_found) = prop_grid(i)
      endif
    enddo
    ! Remove bracketed property specifications
    if (isBracketed) then
      do i=1,n_found
        do j=1,n_grid_not_found
          if (prop_list(i) == prop_grid(j)) then
            do k=j,n_grid_not_found-1
              prop_grid(k) = prop_grid(k+1)
            enddo
            exit
          endif
        enddo
        n_grid_not_found = n_grid_not_found - 1
      enddo
    endif
  end function isAnyPropBracketed

  !-----------------------------------------------------------------------------
  !> Locate the points on the triple line corresponding to the property values
  !! in prop_grid.
  !!
  !! The routine starts at the "triple" point, and then
  !! traverses the triple area line until it has bracketed one or more grid value; it then
  !! locates the exact sublimation point corresponding to this value, before
  !! continuing.
  !!
  !! Solid appearance line is specified, setting betaSol=0.0
  !! Zero liquid line is specified, setting betaSol<0.0
  !!
  !! \author Morten 2018-06
  !-----------------------------------------------------------------------------
  subroutine triple_area_points_based_on_prop(Z,Ttr,Ptr,Pmax,betaSol,n_grid,&
       propflag,prop_grid,T_grid,P_grid,beta_grid,X_grid,Y_grid,n_grid_found,ierr_out)
    use solideos, only: solidComp, nSolid
    implicit none
    ! Input:
    real, intent(in) :: Z(nc)                  !< total composition
    real, intent(in) :: ttr, ptr               !< triple point
    integer, intent(in) :: n_grid              !< number of grid points
    integer, intent(in) :: propflag            !< specified property 1:s, 2:lnv
    real, intent(inout) :: prop_grid(n_grid)   !< property grid
    real, intent(in) :: pMax                   !< Pressure limit (Pa)
    real, intent(in) :: betaSol                !< -1.0 interpreted as maximum solid amount
    ! Output:
    real, intent(out) :: T_grid(n_grid)        !< t at the grid points
    real, intent(out) :: P_grid(n_grid)        !< p at the grid points
    real, intent(out) :: beta_grid(2,n_grid)   !< gas and solid fraction
    real, intent(out) :: X_grid(nc,n_grid)     !< Liquid composition
    real, intent(out) :: Y_grid(nc,n_grid)     !< Gas composition
    integer, intent(out) :: n_grid_found       !< gird points found
    integer, optional, intent(out) :: ierr_out !< error flag; nonzero if error
    ! Locals:
    real :: t, p
    integer :: i,iter,is,ierr,s,j,sl
    real :: curr_prop_val, old_prop_val, ln_spec, dln_s, pn, sspec
    real, dimension(nc+3) :: Xvar, Xold, dXds, Xtemp
    integer :: n_found, nsol, n_grid_not_found, mode
    real, parameter :: dPmax = 0.2e5
    real :: prop_grid_local(n_grid), prop_grid_list(n_grid)      !< property grids
    real :: Pdew, beta
    integer, parameter :: nd = 3
    real :: XX(nph,nc), betaVec(nph), K(nc), bval, bs
    integer :: phaseVec(nph)
    if (present(ierr_out)) ierr_out = 0
    ! Mapping solid-gas line starting at triple point
    if (nSolid /= 1) then
      call stoperror('saturation::solidEnvelopePlot: No solid model initialized!')
    else
      is = solidComp(1)
    endif
    ! Initial state
    t = Ttr
    p = Ptr
    phaseVec(1:nd) = (/VAPPH, LIQPH, SOLIDPH/)
    XX(3,:) = 0.0
    XX(3,is) = 1.0
    betaVec(1:nd) = (/1.0,0.0,0.0/)
    XX(1,:) = Z
    Pdew = safe_dewP(T,XX(2,:),XX(1,:),ierr)
    K = XX(1,:)/XX(2,:)

    ! Compute property value of initial point on sublimation line
    call genericPropertyMP(nd,t,p,Z,XX,betaVec,phaseVec,propflag,old_prop_val)
    prop_grid_local = prop_grid
    n_grid_found = 0
    n_grid_not_found = n_grid
    if (betaSol < 0.0) then
      mode = VARSOLID
      s = nc+2 ! Specify pressure
      nsol = max(int((Pmax-Ptr)/dPmax)+1,50)
    else
      mode = VARFLUID
      s = nc+3 ! Specify beta
      nsol = 100
    endif
    ! Start from lower triple area edge
    beta = 1.0
    bs = 0.0
    do i=1,nsol
      call fillXvar(mode,bs,beta,t,p,K,Xold)
      if (mode == VARSOLID) then
        pn = Ptr + (Pmax-Ptr)*real(i-1)/real(nsol-1)
        ln_spec = log(pn)
        dln_s = ln_spec - log(p)
        !call newton_extrapolate_threePh(Z,K,t,p,beta,bs,s,is,ln_spec,dln_s,mode)
      else ! mode == VARFLUID
        beta = 1.0 - real(i-1)/real(nsol-1)
        ln_spec = beta
      endif
      iter = solidPointOnEnvelopeThreePh(Z,K,t,p,beta,bs,is,s,ln_spec,mode,ierr)
      if (ierr /= 0) then
        if (present(ierr_out)) then
          ierr_out = ierr
          return
        else
          call stoperror("grid error: Not able to solve for triple area-line point")
        endif
      endif
      call fillXvar(mode,bs,beta,t,p,K,Xvar)
      call propFromXvar(Xvar,Z,is,mode,bs,beta,t,p,K,betaVec,XX)
      call genericPropertyMP(nd,t,p,Z,XX,betaVec,phaseVec,propflag,curr_prop_val)
      ! Check to see if sspec is bracketed.
      if ( isAnyPropBracketed(curr_prop_val,old_prop_val,n_grid,prop_grid_local,&
           n_grid_not_found,prop_grid_list,n_found)) then
        do j=1,n_found
          sspec = prop_grid_list(j)
          if ( propflag == locate_from_temperature .OR. &
               propflag == locate_from_pressure) then
            if (propflag == locate_from_temperature) then
              sl = nc + 1
            else
              sl = nc + 2
            endif
            ln_spec = log(sspec)
            Xtemp = Xold + dXdS*(ln_spec - Xold(sl))
            call propFromXvar(Xtemp,Z,is,mode,bs,beta,t,p,K,betaVec,XX)
            iter = solidPointOnEnvelopeThreePh(Z,K,t,p,beta,bs,is,sl,ln_spec,mode,ierr)
          else
            Xtemp = Xvar
            ! sspec is bracketed, so solve for exact property value and return
            if (mode == VARSOLID) then
              bval = beta
            else ! mode == VARFLUID
              bval = betaSol
            endif
            call bracketSolveForProperty(n=nc+3,Z=Z,beta=bval,propflag=propflag,&
                 propspec=sspec,X=Xtemp,Xold=Xold,dXdS=dXdS,s=s,ierr=ierr,mode=mode)
            call propFromXvar(Xtemp,Z,is,mode,bs,beta,t,p,K,betaVec,XX)
          endif
          if (ierr /= 0) then
            if (present(ierr_out)) then
              ierr_out = ierr
              return
            else
              if (ierr == 100) then
                print *,'No solution for bracket solver'
              else
                print *,'ierr',ierr
              endif
              call stoperror("grid error: Bracket solver failed!")
            endif
          endif

          call store_point_triple(nd,t,p,XX,beta,bs,sspec,n_grid,T_grid,P_grid,beta_grid,&
               X_grid,Y_grid,n_grid_found,prop_grid)
          if ( n_grid_found == n_grid ) then
            return
          end if
        enddo
      endif
      ! Update old_prop_val and continue
      old_prop_val = curr_prop_val

    enddo

  end subroutine triple_area_points_based_on_prop

  !-----------------------------------------------------------------------------
  !> Store point in output vector for triple area line
  !>
  !> \author MH, 2018-06
  !-----------------------------------------------------------------------------
  subroutine store_point_triple(nd,t,p,XX,beta,betaSol,prop,n_grid,T_grid,P_grid,&
       beta_grid,X_grid,Y_grid,n_grid_found,prop_grid)
    implicit none
    ! Input:
    integer, intent(in) :: nd                  !< number of phases
    real, intent(in) :: t,p                    !< Two-phase gas phase fraction (mol/mol)
    real, intent(in) :: beta,betaSol           !< Phase fractions
    real, intent(in) :: XX(nph,nc)             !< Phase compositions
    integer, intent(in) :: n_grid              !< number of grid points
    real, intent(in) :: prop                   !< Property located
    ! Output:
    real, intent(inout) :: T_grid(n_grid)        !< t at the grid points
    real, intent(inout) :: P_grid(n_grid)        !< p at the grid points
    real, intent(inout) :: beta_grid(2,n_grid)   !< gas and solid fraction
    real, intent(inout) :: X_grid(nc,n_grid)     !< Liquid composition
    real, intent(inout) :: Y_grid(nc,n_grid)     !< Gas composition
    integer, intent(inout) :: n_grid_found       !< Number of grid points found
    real, intent(inout) :: prop_grid(n_grid)     !< Property grid
    ! Locals:
    n_grid_found = n_grid_found + 1
    T_grid(n_grid_found) = t
    P_grid(n_grid_found) = p
    beta_grid(:,n_grid_found) = (/ (1.0-betaSol)*beta, betaSol/)
    Y_grid(:,n_grid_found) = XX(1,:)
    X_grid(:,n_grid_found) = XX(2,:)
    prop_grid(n_grid_found) = prop
  end subroutine store_point_triple

  !-----------------------------------------------------------------------------
  !> Get properties from Xvar
  !>
  !> \author MH, 2018-06
  !-----------------------------------------------------------------------------
  subroutine propFromXvar(Xvar,Z,is,mode,betaSol,beta,t,p,K,betaVec,XX)
    implicit none
    ! Input:
    real, dimension(nc+3), intent(in) :: Xvar !< Variable vector
    real, dimension(nc), intent(in) :: Z      !< Overall composition
    real, intent(inout) :: betaSol,beta       !< Solid and gas phase fraction (mol/mol)
    integer, intent(in) :: mode               !< VARSOLID or VALFLUID
    integer, intent(in) :: is                 !< solid index
    ! Output:
    real, intent(out) :: t                    !< temperature
    real, intent(out) :: p                    !< pressure
    real, intent(out) :: K(nc)                !< equilibrium constants
    real, intent(out) :: betaVec(nph)         !< Phase mol fractions
    real, intent(inout) :: XX(nph,nc)         !< Phase compositions
    !
    real :: Zstar(nc)
    K = exp(Xvar(1:nc))
    t = exp(Xvar(nc+1))
    p = exp(Xvar(nc+2))
    if (mode == VARSOLID) then
      betaSol = Xvar(nc+3)
      betaVec(1) = (1.0-betaSol)*beta
      betaVec(3) = betaSol
    else ! mode == VARFLUID
      beta = Xvar(nc+3)
      betaVec(1) = (1.0-betaSol)*beta
      betaVec(3) = betaSol
    endif
    betaVec(2) = max(0.0,1.0 - betaVec(1) - betaVec(3))
    Zstar = Z
    Zstar(is) = Z(is) - betaSol
    Zstar = Zstar/(1.0-betaSol)
    XX(2,:) = Zstar/(1-beta+beta*K)
    XX(1,:) = K*XX(2,:)
  end subroutine propFromXvar

  !-----------------------------------------------------------------------------
  !> Fill Xvar from properties
  !>
  !> \author MH, 2018-06
  !-----------------------------------------------------------------------------
  subroutine fillXvar(mode,betaSol,beta,t,p,K,Xvar)
    implicit none
    ! Input:
    real, intent(in) :: t                      !< temperature
    real, intent(in) :: p                      !< pressure
    real, intent(in) :: K(nc)                  !< equilibrium constants
    real, intent(in) :: betaSol,beta           !< Solid and gas phase fraction (mol/mol)
    integer, intent(in) :: mode                !< VARSOLID or VALFLUID
    ! Output:
    real, dimension(nc+3), intent(out) :: Xvar !< Variable vector
    !
    if (mode == VARSOLID) then
      Xvar = (/ log(K), log(t), log(p), betaSol/)
    else ! mode == VARFLUID
      Xvar = (/ log(K), log(t), log(p), beta/)
    endif
  end subroutine fillXvar

  subroutine genericPropertyMP(nd,t,p,Z,XX,betaVec,phaseVec,propflag,prop)
    use vls, only: mpSpecificVolume, mpEnthalpy, mpEntropy
    implicit none
    integer, intent(in) :: nd !< Phase identifier
    integer, intent(in) :: phaseVec(nph) !< Phase identifier
    real, intent(in) :: betaVec(nph) !< Phase fractions
    real, intent(in) :: XX(nph,nc) !< Phase compositions
    real, intent(in) :: t !< K - Temperature
    real, intent(in) :: p !< Pa - Pressure
    real, dimension(1:nc), intent(in) :: Z !< Compozition
    integer, intent(in) :: propflag !< Flag determining what property to return
    real, intent(out) :: prop !< Property

    if (propflag == locate_from_entropy) then
      prop = mpEntropy(nd,t,p,betaVec,XX,phaseVec)
    else if (propflag == locate_from_lnvol) then
      prop = mpSpecificVolume(nd,t,p,betaVec,XX,phaseVec)
      prop = log(prop)
    else if (propflag == locate_from_enthalpy) then
      prop = mpEnthalpy(nd,t,p,betaVec,XX,phaseVec)
    else if (propflag == locate_from_temperature) then
      prop = T
    else if (propflag == locate_from_pressure) then
      prop = P
    end if

  end subroutine genericPropertyMP

  !-----------------------------------------------------------------------------
  !> Function used to solve for point on sublimation line having property value
  !> propspec.
  !>
  !> \author MH 2018-6
  !-----------------------------------------------------------------------------
  function propertyFunctionWrapperTriple(Xs,param) result(fun)
    use solideos, only: solidComp
    implicit none
    real, intent(inout) :: Xs
    real, dimension(4*nc+14), intent(inout) :: param
    real :: fun
    ! Locals:
    integer :: iter, s, is, ierr
    real, dimension(nc) :: Z, K
    real, dimension(nc+3) :: X,Xold,dXdS
    real :: t, prop, propspec, ln_spec, p, beta, betaSol
    integer :: propflag, mode
    real :: XX(nph,nc), betaVec(nph)
    integer :: phaseVec(nph)
    integer, parameter :: nd = 3
    ! Unpack the param vector
    propspec = param(1)
    mode = int(param(4*nc+14))
    s = int(param(3))
    Z = param(4:nc+3)
    Xold = param(nc+4:2*nc+6)
    dXdS = param(2*nc+7:3*nc+9) ! is constant each time function is called
    X = param(3*nc+10:4*nc+12)
    propflag = int(param(4*nc+13))
    !
    is = solidComp(1)

    ! Extrapolate for better initial values in sat_newton
    X = Xold !+ dXdS*(Xs-Xold(s))
    if (mode == VARSOLID) then
      beta = param(2)
      betaSol = X(nc+3)
    else
      betaSol = param(2)
      beta = X(nc+3)
    endif
    call propFromXvar(X,Z,is,mode,betaSol,beta,t,p,K,betaVec,XX)

    ! Find point on envelope corresponding to ln(property value) equal to Xs
    ln_spec = Xs
    iter = solidPointOnEnvelopeThreePh(Z,K,t,p,beta,betaSol,is,s,ln_spec,mode,ierr)
    call fillXvar(mode,betaSol,beta,t,p,K,X)
    call propFromXvar(X,Z,is,mode,betaSol,beta,t,p,K,betaVec,XX)
    phaseVec(1:nd) = (/VAPPH, LIQPH, SOLIDPH/)
    XX(3,:) = 0.0
    XX(3,is) = 1.0

    ! Update Xold-slot in param vector with new X value
    param(3*nc+10:4*nc+12) = X

    ! Calculate the property value at the new saturation point
    call genericPropertyMP(nd,t,p,Z,XX,betaVec,phaseVec,propflag,prop)

    ! Compute new value of objective function
    fun = (prop - propspec)/max(abs(propspec), 1.0)
  end function propertyFunctionWrapperTriple

end module saturation_point_locators

max_stars_repo_name: jabirali/thermopack

max_issues_repo_name: jabirali/thermopack

max_forks_repo_name: jabirali/thermopack


text:
      SUBROUTINE WF_TI_MTRX
     &          (isotropic)

c     Computes t - the coefficient matrix of the linear system of
c                         dp/dz=-ik*t*p.

      implicit real*8 (a-h,o-z)

c     Wave Fields
      real     *  8 dtheta,lub,thtinc,alpha,h,prec
      complex  * 16 c1,s1,c2,s2,p1,p2,dp1dh,dp2dh,
     &              t11,t31,t41,t12,t32,t42,t14,t34,t44,
     &              tm11,tm31,tm41,tm12,tm32,tm42,tm14,tm34,tm44

      common/wf_iter/
     &              dtheta(2),lub(2),thtinc,alpha,h,prec,
     &              iter_flag,mxiter,nriter
     &      /wf_pmtx/
     &              c1,c2,s1,s2,p1(8),p2(8),dp1dh(8),dp2dh(8)
     &      /wf_tmtx/
     &              t11,t31,t41,t12,t32,t42,t14,t34,t44,
     &              tm11,tm31,tm41,tm12,tm32,tm42,tm14,tm34,tm44

      t41=tm41
      t32=tm32+c2**2
      t14=1.d0-s2**2*tm14
      t11=-tm11*s2
      t44=-tm44*s2
      if (isotropic .eq. 0) then
         t31=tm31
         t12=tm12*s2
         t42=tm42
         t34=tm34*s2
      end if
      RETURN
      END      ! WF_TI_MTRX

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


text:
*     TB04CD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, MMAX, PMAX, NPZMAX
      PARAMETER        ( NMAX = 20, MMAX = 20, PMAX = 20,
     $                   NPZMAX = NMAX )
      INTEGER          PMNMAX
      PARAMETER        ( PMNMAX = PMAX*MMAX*NPZMAX )
      INTEGER          LDA, LDB, LDC, LDD, LDGAIN, LDNP, LDNZ
      PARAMETER        ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
     $                   LDD = PMAX, LDGAIN = PMAX, LDNP = PMAX,
     $                   LDNZ = PMAX )
      INTEGER          LIWORK
      PARAMETER        ( LIWORK = NMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = NMAX*( NMAX + PMAX ) +
     $                            MAX( NMAX + MAX( NMAX, PMAX ),
     $                                 NMAX*( 2*NMAX + 3 ) ) )
*     .. Local Scalars ..
      DOUBLE PRECISION TOL
      INTEGER          I, IJ, INFO, J, K, M, N, NPZ, P
      CHARACTER*1      JOBD, EQUIL
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
     $                 D(LDD,MMAX), DWORK(LDWORK), GAINS(LDGAIN,MMAX),
     $                 POLESI(PMNMAX), POLESR(PMNMAX), ZEROSI(PMNMAX),
     $                 ZEROSR(PMNMAX)
      INTEGER          IWORK(LIWORK), NP(LDNP,MMAX), NZ(LDNZ,MMAX)
*     .. External Subroutines ..
      EXTERNAL         TB04CD
*     .. Intrinsic Functions ..
      INTRINSIC        MAX
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, M, P, TOL, JOBD, EQUIL
      NPZ = N
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99992 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
            WRITE ( NOUT, FMT = 99991 ) M
         ELSE
            READ ( NIN, FMT = * ) ( ( B(I,J), I = 1,N ), J = 1,M )
            IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
               WRITE ( NOUT, FMT = 99990 ) P
            ELSE
               READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
               READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,P )
*              Find the transfer matrix T(s) of (A,B,C,D) in the
*              pole-zero-gain form.
               CALL TB04CD( JOBD, EQUIL, N, M, P, NPZ, A, LDA, B, LDB,
     $                      C, LDC, D, LDD, NZ, LDNZ, NP, LDNP, ZEROSR,
     $                      ZEROSI, POLESR, POLESI, GAINS, LDGAIN, TOL,
     $                      IWORK, DWORK, LDWORK, INFO )
*
               IF ( INFO.NE.0 ) THEN
                  WRITE ( NOUT, FMT = 99998 ) INFO
               ELSE
                  WRITE ( NOUT, FMT = 99997 )
                  DO 60 J = 1, M
                     DO 40 I = 1, P
                        IJ = ( (J-1)*P + I-1 )*NPZ + 1
                        IF ( NZ(I,J).EQ.0 ) THEN
                           WRITE ( NOUT, FMT = 99996 ) I, J
                        ELSE
                           WRITE ( NOUT, FMT = 99995 ) I, J,
     $                        ( ZEROSR(K), ZEROSI(K),
     $                                 K = IJ,IJ+NZ(I,J)-1 )
                        END IF
                        WRITE ( NOUT, FMT = 99994 ) I, J,
     $                     ( POLESR(K), POLESI(K), K = IJ,IJ+NP(I,J)-1 )
                        WRITE ( NOUT, FMT = 99993 ) I, J, ( GAINS(I,J) )
   40                CONTINUE
   60             CONTINUE
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' TB04CD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TB04CD = ',I2)
99997 FORMAT (/' The poles, zeros and gains of the transfer matrix',
     $         ' elements: ')
99996 FORMAT (/' no zeros for element (',I2,',',I2,')')
99995 FORMAT (/' zeros of element (',I2,',',I2,') are ',//
     $         '   real part     imag part '// (2X,F9.4,5X,F9.4))
99994 FORMAT (/' poles of element (',I2,',',I2,') are ',//
     $         '   real part     imag part '// (2X,F9.4,5X,F9.4))
99993 FORMAT (/' gain of element (',I2,',',I2,') is ', F9.4)
99992 FORMAT (/' N is out of range.',/' N = ',I5)
99991 FORMAT (/' M is out of range.',/' M = ',I5)
99990 FORMAT (/' P is out of range.',/' P = ',I5)
      END

ext: f

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference


text:
C
C	$Id: gqfaci.f,v 1.5 2008-07-27 00:20:59 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE GQFACI(ERRIND,COLI)
C
C  INQUIRE FILL AREA COLOUR INDEX
C
      include 'gkscom.h'
C
      INTEGER ERRIND,COLI
C
C  Check if GKS is in the proper state.
C
      CALL GZCKST(8,-1,ERRIND)
      IF (ERRIND.EQ.0) THEN
       COLI = CFACI
      ELSE
       COLI = -1
      ENDIF
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
About Me

Name: Jeff

Details
Phone Number: Long distance at the moment

Office: Beijing

Personality: normal

Location: address(1600 Pennsylvania Avenue, Washington, DC 20500)

signed up in April of 2007.

Say something

20070412 22:52:29 nbsp Beijing is prettyUsers/jeff

20070412 22:52:48 nbsp Asians are good cooks Users/jeff

20070413 16:07:25 nbsp Nice HMDS Users/StevenDaubert

20070413 18:07:26 nbsp hey jeff, that image is likely copyright of someone other than you seeing as it said collegehumor.com on it. Almost everything on the Davis Wiki is under the Creative Commons License (see the bottom of the page), so in order to post content belonging to other people you need to obtain their permission. Users/ArlenAbraham

20070415 21:35:08 nbsp Thanks for the help i am still new to wikis Users/jeff

20111017 06:47:20 nbsp Hi Jeff! Users/AshleyHamidi

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


      module sibalb_coeff
      implicit none

      contains

      FUNCTION COEFFSIB(TABLE, NTABL, LAI ,DX, DY)

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NTABL, LAI

      REAL, INTENT(IN) :: DX, DY
      REAL, INTENT(IN), DIMENSION(NTABL,2) :: TABLE
      REAL COEFFSIB

      COEFFSIB = (TABLE(LAI,  1)                                                  &
             + (TABLE(LAI  ,2) - TABLE(LAI  ,1)) * DY ) * (1.0-DX)             &
             + (TABLE(LAI+1,1)                                                 &
             + (TABLE(LAI+1,2) - TABLE(LAI+1,1)) * DY ) * DX

      RETURN
      END FUNCTION COEFFSIB

    end module sibalb_coeff


text:
!
!     function dex2: formate of basis function array in setsto3/4
!     for EHT calculations
!
!     part of QMDFF
!
double precision function dex2(m)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IF(M .LT. 2) THEN
   DEX2=1
ELSE
   DEX2=1
   do I=1,M,2
      DEX2=DEX2*I
   end do
ENDIF
RETURN
end function dex2

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
C ==========================================================
      real Function NLnFnc(
C ==========================================================
C group (F)
     &       nU, U,
C group (ANI)
     &       XYP, IPE, IEE, HesP, hStar, status,
     &       lE, iEs, XYPs, IPEs, HesPs, detGs, qEs,
     &       nPw, nEw, XYPw, HesPw, IPEw,
     &       MetricFunction, flagAnalytic,
     &       iSE, rSE, iP1, iFNCs, calCrv,
     &       L1Et, L2Et, tE,
     &       nL2t, nStept, nEt, nCrvFnc, LFnc, ILt)
C ==========================================================
      include 'status.fd'
C ==========================================================
C group (F)
      real U(*)

C group (ANI)
      Integer IPE(3, *), IEE(3, *)
      Integer iEs(*), IPEs(3, *), IPEw(3, *), iSE(*), status
      real  XYP(2, *), HesP(3, *), hStar

      real  XYPs(2),   HesPs(*),  detGs,  qEs(*)
      real  XYPw(2, *), HesPw(3, *), rSE(*)
      real  tE(*)

      EXTERNAL calCrv

      Integer L1Et(2, *), L2Et(*)
      Integer nL2t(*), nStept(4, *), nEt(*)
      Integer LFnc(*), ILt(*)

      Logical flagAnalytic

      Integer  MetricFunction
      EXTERNAL MetricFunction

C group (Local variables)
      Integer ip(4)
      real  prjXYPs(2), XYPo(2), tc, qEt
      Logical check1j, ifXnode

C ==========================================================
      ip(1) = 1
      ip(2) = 2
      ip(3) = 3
      ip(4) = 1

      If(nU.EQ.2) Then
         Do i = 1, 2
            XYPs(i) = U(i)
         End do
      Else if(nU.EQ.1) Then
         tc = U(1)

         Call aniCrv(tc, XYPs, iFNCs, calCrv)
      End if

      Call findSE(nCrvFnc, LFnc, iFNCs, k)
      If(k.GT.0) Then
         ir = ILt(k)
         Call prjCrv(XYPs, prjXYPs, iFNCs, tc, calCrv,
     &               L1Et(1, ir), L2Et(ir), nL2t(k), nStept(1, k),
     &               nEt(k), tE(ir))

         If(.NOT.flagAnalytic) Then
         Call LINTRP(nEw, IPEw, nPw, XYPw, 3, HesPw, 1,
     &               prjXYPs, HesPs, iSE, rSE, .FALSE.)
      Else
            Call scaleBack(prjXYPs, XYPo)
            Call iniQ_analytic(1, XYPo, MetricFunction, HesPs)
         End if
      Else
         If(.NOT.flagAnalytic) Then
         Call LINTRP(nEw, IPEw, nPw, XYPw, 3, HesPw, 1,
     &               XYPs, HesPs, iSE, rSE, .FALSE.)
         Else
            Call scaleBack(XYPs, XYPo)
            Call iniQ_analytic(1, XYPo, MetricFunction, HesPs)
         End if
      End if

      Call calDet(HesPs, detGs)

      Do 50 n = 1, lE
         iE = iEs(n)
         If(iE.LE.0) goto 50

         Do i1 = 1, 3
            If(IPEs(i1, n).EQ.iP1) Then
               i2 = ip(i1 + 1)
               i3 = ip(i2 + 1)

               iPa = IPEs(i2, n)
               iPb = IPEs(i3, n)

               Call calQE(
     &              HesPs,        XYPs,
     &              HesP(1, iPa), XYP(1, iPa),
     &              HesP(1, iPb), XYP(1, iPb),
     &              hStar, qEs(n))

c  ...  updating the quality
               If(ifXnode(status, ANISmoothMesh)) Then
                  Do 10 j1 = 1, 3
                     jE = IEE(j1, iE)
                     If(jE.LE.0) goto 10

                     jP1 = IPE(1, jE)
                     jP2 = IPE(2, jE)
                     jP3 = IPE(3, jE)

                     If(check1j(iP1, IPE(1, jE))) Then
                        if(jP2.EQ.iP1) jP2 = jP1
                        if(jP3.EQ.iP1) jP3 = jP1

                        Call calQE(
     &                       HesPs,        XYPs,
     &                       HesP(1, jP2), XYP(1, iPa),
     &                       HesP(1, jP3), XYP(1, iPb),
     &                       hStar, qEt)

                        qEs(n) = min(qEs(n), qEt)
                     Else
                        Call calQE(
     &                       HesP(1, jP1), XYPs,
     &                       HesP(1, jP2), XYP(1, iPa),
     &                       HesP(1, jP3), XYP(1, iPb),
     &                       hStar, qEt)

                        qEs(n) = min(qEs(n), qEt)

                        Call calQE(
     &                       HesPs,        XYP(1, jP1),
     &                       HesP(1, iPa), XYP(1, jP2),
     &                       HesP(1, iPb), XYP(1, jP3),
     &                       hStar, qEt)

                        Call findSE(lE, iEs, jE, m)
                        If(m.LE.0) Call errMes(6002,
     &                                  'NLnFnc', 'System error')
                        qEs(m) = min(qEs(m), qEt)
                     End if
 10               Continue
               End if

               goto 50
            End if
         End do
 50   Continue

      NLnFnc = 1D0
      Do n = 1, lE
         NLnFnc = min(NLnFnc, qEs(n))
      End do
      NLnFnc = 1D0 - NLnFnc
      Return
      End


text:
! ======================================================================
! NAME
!
!     NONLINEAR
!
! DESCRIPTION
!
!     NONLINEAR is a supporting program of GPECN1DTDPBC that generates
!     a uniform, time-independent, and real-valued coupling constant
!     that parameterizes the strength of a (simple, single-component)
!     Bose-Einstein condensate’s nonlinear mean-field interaction
!     arising from two-body contact-like interactions between atoms
!     within the condensate.
!
! AUTHOR
!
!     Marty Kandes, Ph.D.
!     Computational & Data Science Research Specialist
!     User Services Group
!     San Diego Supercomputer Center
!     University of California, San Diego
!
! COPYRIGHT
!
!     Copyright (c) 2010 - 2021 Martin Charles Kandes
!
! LICENSE
!
!     The MIT License (MIT)
!
! LAST UPDATED
!
!     Tuesday, April 4th, 2017
!
! ----------------------------------------------------------------------
      PROGRAM NONLINEAR
      IMPLICIT NONE

!     Parameter Declarations:
!     -----------------------
      REAL, PARAMETER :: PI = 3.1415926535897932384626433832795028841971E0

!     Input Variable and Array Declarations:
!     --------------------------------------
      INTEGER :: NUMBER_OF_GRID_POINTS
      INTEGER :: NONLINEAR_COUPLING_SWITCH

      REAL :: RADIUS
      REAL :: NONLINEAR_COUPLING_AMPLITUDE

!     Internal Variable and Array Declarations:
!     -----------------------------------------
      CHARACTER(80) :: BUFFER

      INTEGER :: CURRENT_GRID_POINT
      INTEGER :: LENGTH
      INTEGER :: INFO

      REAL :: ANGULAR_GRID_SPACING_SIZE

      COMPLEX, ALLOCATABLE, DIMENSION(:) :: NONLINEAR_COUPLING

!     Input Variable Assignment Read from Command-line Arguments:
!     -----------------------------------------------------------
      CALL GET_COMMAND_ARGUMENT(1,BUFFER,LENGTH,INFO)
      READ(BUFFER,*) NUMBER_OF_GRID_POINTS

      CALL GET_COMMAND_ARGUMENT(2,BUFFER,LENGTH,INFO)
      READ(BUFFER,*) NONLINEAR_COUPLING_SWITCH

      CALL GET_COMMAND_ARGUMENT(3,BUFFER,LENGTH,INFO)
      READ(BUFFER,*) RADIUS

      CALL GET_COMMAND_ARGUMENT(4,BUFFER,LENGTH,INFO)
      READ(BUFFER,*) NONLINEAR_COUPLING_AMPLITUDE

!     Allocate Arrays:
!     ----------------
      ALLOCATE(NONLINEAR_COUPLING(NUMBER_OF_GRID_POINTS))

!     Internal Variable Assignment:
!     -----------------------------
      ANGULAR_GRID_SPACING_SIZE = 2.0E0*PI/FLOAT(NUMBER_OF_GRID_POINTS)

!     Main Program:
!     -------------
      IF (NONLINEAR_COUPLING_SWITCH.EQ.0) THEN
         NONLINEAR_COUPLING = CMPLX(NONLINEAR_COUPLING_AMPLITUDE,0.0E0)
      ENDIF

      OPEN(UNIT=998,ACTION='WRITE',FORM='UNFORMATTED')
      WRITE(UNIT=998) NONLINEAR_COUPLING
      CLOSE(UNIT=998)

!     Deallocate Arrays:
!     ------------------
      DEALLOCATE(NONLINEAR_COUPLING)

      STOP
      ENDPROGRAM NONLINEAR
! ======================================================================


text:
program print_asterisks
    implicit none
    integer :: i = 123456
    character(len=6) :: c = 'abcdef'

    print '(I8)', i
    print '(I4)', i

    print '(A8)', c
    print '(A4)', c

end program print_asterisks

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
menjacket2 is a spammer who likes yoga pants and whose account has been BannedGroup shuttered.
 In menjacket2s defense, who doesnt like yoga pants?
  The poor, defenseless yogans, whose hide is ruthlessly harvested in increasing numbers every year by the faceless yogaindustrial complex?


text:
Terminology
Mixed Race is a term used by people who wish to acknowledge their mixed race heritage.  The term mixed race is used as opposed to biracial because people may be more than two races.  Biracial tends to imply that one parent is one race and the other parent is a different race and can be leading, but both parents could be of the same mixed heritage (such as both parents being White and Black, themselves biracial and thus their children are biracial).

Multiethnic is a term to acknowledge the mixed culture and history they have grown up with.  Anyone who has to identify as (African, Asian, Mexican, Brazilian, Indian, etc)American will have a multiethnic background because they are negotiating a mixed of cultures between their home/native/ancestral culture and the U.S. culture.  Multiethnic can also be used by people who grow up as one race but identifies with another races culture (such as a White person who identifies with Hip Hop or Mexican culture).

History

First and foremost, race and ethnicity are social constructs.  There is no race gene that determines ones race.  Race and ethnicity will vary in different locations and social arenas.  The definitions provided above...

Many white people may actually be mixed race if their heritage comes from different European nations, but because white has been constructed to include most of Western Europe, people will just claim White instead of each heritage.  The same may be applied to Black people whose mixed race heritage may be lost due to the social construct of Black.

There is a long history of Mixed Race that needs to be researched and written about.

When Africans were forced to come to the U.S. as slaves, there was racial mixing between Whites and Blacks that lead to offspring.  Because of the One Drop Rule, any person who had even one drop of Black blood would be racialized as Black.  When the Chinese, Japanese, Pilipinos and Koreans settled in the U.S., they were largely bachelor societies with very few women because of exclusionary laws passed by the U.S. government.  Asian then began to marry White, Black, Mexican and Native Americans to start a family.

JapaneseAmericans have the highest rate of outmarriages amongst the Asian American community.  Outmarriage means to marry someone different from ones own race.

The 2000 Census was the first time that people were allowed to mark off more than one box for racial identification.  Prior to this, people were forced to select only one box.

Lamar Heystek, UCD alumni who is mixed race with Dutch and Indonesian ancestry, currently sits on the Davis City Council.

University

Currently at UCD, students are allowed to check all of the races they identify with.  BUT when the university inputs this info into their system, the students race with least representation at the school is marked as their only race.  So unfortunately, a statistic of students identifying as mixed race is not available.

Example:  Users/JoAnnaRich is mixed and marked off the Pilipino, White/Caucasian, Native American boxes but on University records will be marked as Native American.
Other Examples:
  Black and White > Black
  Mexican and Black > Black
  Asian and Mexican > Mexican

The Cross Cultural Center has a multiethnic intern.

CAPS has a Multicultural Immersion Program (MIP) featuring workshops such as Check One Box: Mixed Race Identity and Mixing It Up: Interracial Dating and Relationships.

Classes

ASA 120  Multiracial Asian Pacific American Issues
NAS 134  Race and Sex: Race Mixture and Mixed Peoples
AMS 154  Race, Culture and Society in the United States
AAS 100  Survey of Ethnicity in the U.S.
AAS 156  Language and Identity in Africa and the African Diaspora
AAS 180  Race and Ethnicity in Latin America
HIS 7A  History of Latin America to 1700
HIS 7B  History of Latin America, 17001900
HIS 7C  History of Latin America, 1900present

Student Organizations and Events

Mixed Student Union Mixed Student Union (MSU) is the only student organization on campus devoted to the mixed race, multiethnic, and transracial adoptee community.
One Drop is an AS Papers publication that focuses on mixed heritage.
Mixed Heritage Week
Motley Movie Medley  Annual showcase of films for, by and about mixed heritage people.

Library

Shields Library offers a few books about mixed heritage in the E184 section.  Currently all off the MAVIN magazines, which focus on mixed heritage issues, is being bound and will be available again in Spring quarter.

Professors
Several UCD professors are mixed heritage, research, teach classes, and/or have written books about being mixed:
Kevin Johnson King Law School  How Did You Get To Be A Mexican?; A Reader on Race, Civil Rights, and the Law: A Multiracial Approach; Mixed Race America and the Law: A Reader
Jack Forbes  Native American Studies
Caroline Valverde  Asian American Studies
Clarence Walker  History
David G. Smith  Anthropology
Yvette FloresOrtiz  Chicana/o Studies

Graduate Students
Tamu Nolfo is biracial PhD candidate and is currently doing research on mixed race youth.
Isis L. Golden is a PhD candidate in Dramatic Art/Performance Studies.  http://theatredance.ucdavis.edu/grad_prog/gradbio_igolden.html The majority of her research in contemporary Native American Theater, explores issues of biracial/multiracial ethnicity and sovereignty and autonomy.


MODULE module_sf_sstskin

CONTAINS

   SUBROUTINE sst_skin_update(xland,glw,gsw,hfx,qfx,tsk,ust,emiss, &
                dtw1,sstsk,dt,stbolt, &
                ids, ide, jds, jde, kds, kde, &
                ims, ime, jms, jme, kms, kme, &
                its, ite, jts, jte, kts, kte )

   USE module_wrf_error
   IMPLICIT NONE

   INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
                                     ims, ime, jms, jme, kms, kme, &
                                     its, ite, jts, jte, kts, kte

   REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: xland, glw, gsw
   REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: hfx, qfx, tsk
   REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: ust, emiss
   REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT ) :: dtw1
   REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT ) :: sstsk
   REAL, INTENT(IN ) :: DT
   REAL, INTENT(IN ) :: STBOLT

   REAL :: lw, sw, q, qn, zeta, dep, dtw3, skinmax, skinmin
   REAL :: fs, con1, con2, con3, con4, con5, zlan, q2, ts, phi, qn1
   REAL :: usw, qo, swo, us, tb, dtc, dtw, alw, dtwo, delt, f1
   INTEGER :: i, j, k

   INTEGER , PARAMETER :: n=1152
   REAL , PARAMETER :: z1=3.,an=.3,zk=.4,rho=1.2,rhow=1025.,cw=4190.
   REAL , PARAMETER :: g=9.8,znuw=1.e-6,zkw=1.4e-7,sdate=1201.6667
      skinmax=-9999.
      skinmin=9999.
      do i=its,ite
      do j=jts,jte
      if(xland(i,j).ge.1.5) then
      qo=glw(i,j)-emiss(i,j)*stbolt*(sstsk(i,j)**4)-2.5e6*qfx(i,j)-hfx(i,j)
      swo=gsw(i,j)
      us=MAX(ust(i,j), 0.01)
      tb=tsk(i,j)-273.15
      dtwo=dtw1(i,j)
      delt=dt
      q=qo/(rhow*cw)
      sw=swo/(rhow*cw)
      f1=1. -0.27*exp(-2.8*z1)-0.45*exp(-0.07*z1)
      dtc=0.0
      alw=1.e-5*max(tb,1.)
      con4=16.*g*alw*znuw**3/zkw**2
      usw=sqrt(rho/rhow)*us
      con5=con4/usw**4
      q2=max(1./(rhow*cw),-q)
      zlan=6./(1.+(con5*q2)**0.75)**0.333
      dep=zlan*znuw/usw
      fs=0.065+11.*dep-(6.6e-5/dep)*(1.-exp(-dep/8.e-4))
      fs=max(fs,0.01)
      dtc=dep*(q+sw*fs)/zkw
      dtc=min(dtc,0.)
      dtw=0.0
      alw=1.e-5*max(tb,1.)
      con1=sqrt(5.*z1*g*alw/an)
      con2=zk*g*alw
      qn=q+sw*f1
      usw=sqrt(rho/rhow)*us
      if(dtwo.gt.0..and.qn.lt.0.) then
         qn1=sqrt(dtwo)*usw**2/con1
         qn=max(qn,qn1)
      endif
      zeta=z1*con2*qn/usw**3
      if(zeta.gt.0.) then
         phi=1.+5.*zeta
      else
         phi=1./sqrt(1.-16.*zeta)
      endif
      con3=zk*usw/(z1*phi)
      dtw=(dtwo+(an+1.)/an*(q+sw*f1)* &
                          delt/z1)/(1.+(an+1.)*con3*delt)
      dtw=max(0.,dtw)
      dtwo=dtw
      ts = tb + dtw + dtc
      skinmax=amax1(skinmax,ts-tb)
      skinmin=amin1(skinmin,ts-tb)
      sstsk(i,j)=ts+273.15
      dtw1(i,j)=dtw
      endif
      end do
      end do
      return
   END SUBROUTINE sst_skin_update
END MODULE module_sf_sstskin