text:
module single_component
  !> Interface to single-component equations of state
  use eos_parameters, only: single_eos
  implicit none
  save

  public :: Zfac_single, enthalpy_single, entropy_single
  public :: lnfugCoeff_single, Gres_single, pressure_single
  public :: Fid_single, Fres_single

contains

  !---------------------------------------------------------------------- >
  !> This function calculates the Compressibility factor with derivatives.
  !!
  !! \param T The temperature [K]
  !! \param P The pressure [Pa]
  !! \param Z The mole fraction [-]
  !! \param phase The phase, 1=liquid, 2=vapour
  !! \param Zfac The Z-factor [-]
  !! \param dZdt Temperature derivative [1/K]
  !! \param dZdp Pressure derivative [1/Pa]
  !! \param dZdz Composition derivative [-]
  !!
  !!
  !! \author MH
  subroutine Zfac_single(nc,seos,T,P,Z,phase,Zfac,dZdt,dZdp,dZdz)
    use eosdata
    use mbwr_additional, only: mbwr_zfac, mbwr_volume
    implicit none
    integer, intent(in) :: nc
    class(single_eos), intent(inout) :: seos
    real, dimension(nc), intent(in) :: Z
    real, intent(in) :: T, P
    integer, intent(in) :: phase
    real, intent(out) :: Zfac
    real, optional, intent(out) :: dZdt, dZdp
    real, optional, dimension(nc), intent(out) :: dZdz
    ! Locals
    real :: v
    !
    !-------- Specific for each equation of state ------------------------
    Choice_EoS: select case (seos%subeosidx) ! choose the Equation of State
    case (meosMbwr19, meosMbwr32)
      v = mbwr_volume (t,p,Z(1),phase,seos%mbwr_meos(1))
      call MBWR_zfac(seos%mbwr_meos(1), t, p, v, Z(1), zfac, dzdt,dzdp,dzdz)
    case (meosNist)
      call seos%nist(1)%meos%calc_zfac(t,p,Z,phase, zfac, dzdt, dzdp, dzdz)
    end select Choice_EoS

  end subroutine Zfac_single

  !---------------------------------------------------------------------- >
  !> This function calculates the Enthalpy and the derivatives.
  !!
  !! \param T The temperature [K]
  !! \param P The pressure [Pa]
  !! \param Z The mole fraction [-]
  !! \param phase The phase, 1=liquid, 2=vapour
  !! \param residual Return only the residual part if .true.
  !! \param enthalpy Enthalpy [J/mol]
  !! \param dhdt Temperature derivative [J/moleK]
  !! \param dhdp Pressure derivative [J/molePa]
  !! \param dhdz Composition derivative [J/mole^2]
  !!
  !! \author MH
  subroutine enthalpy_single(nc,comp,seos,T,P,Z,phase,residual,&
       enthalpy,dhdt,dhdp,dhdz)
    use eosdata
    use compdata, only: gendata_pointer
    use mbwr_additional, only: mbwr_hres, mbwr_volume
    use ideal, only: Hideal_mix
    implicit none
    integer, intent(in) :: nc
    type (gendata_pointer), intent(in), dimension(nc) :: comp
    class(single_eos), intent(inout) :: seos
    real, dimension(nc), intent(in) :: Z
    real, intent(in) :: T, P
    real, intent(out) :: enthalpy
    real, optional, intent(out) :: dhdt, dhdp
    real, optional, dimension(nc), intent(out) :: dhdz
    integer, intent(in) :: phase
    logical, intent(in) :: residual
    ! Locals
    real :: h_id, dhdt_id, dhdp_id
    real, dimension(nc) :: dhdz_id
    real :: v
    !
    !-------- Specific for each equation of state ------------------------
    Choice_EoS: select case (seos%subeosidx) ! choose the Equation of State
    case (meosMbwr19, meosMbwr32)
      if (.not. residual) then
        call Hideal_mix(nc, comp, T, Z, h_id, dhdt, dhdp, dhdz)
        if(present(dhdt)) dhdt_id = dhdt
        if(present(dhdp)) dhdp_id = dhdp
        if(present(dhdz)) dhdz_id = dhdz
      endif
      v = mbwr_volume (t,p,Z(1),phase,seos%mbwr_meos(1))
      call MBWR_Hres(seos%mbwr_meos(1), t, p, v, Z(1), enthalpy, dhdt,dhdp,dhdz)
      if (.not. residual) then
        enthalpy = enthalpy + h_id
        if(present(dhdt)) dhdt = dhdt + dhdt_id
        if(present(dhdp)) dhdp_id = dhdp + dhdp_id
        if(present(dhdz)) dhdz_id = dhdz + dhdz_id
      endif
     case (meosNist)
      call seos%nist(1)%meos%calc_enthalpy(t, p, Z, phase, enthalpy, &
           dhdt, dhdp, dhdz, residual=residual)
    end select Choice_EoS
  end subroutine enthalpy_single

  !---------------------------------------------------------------------- >
  !> This function calculates the Entropy and the derivatives.
  !!
  !! \param T The temperature [K]
  !! \param P The pressure [Pa]
  !! \param Z The overall mole fraction [-]
  !! \param phase The phase, 1=liquid, 2=vapour
  !! \param residual Return only the residual part if .true.
  !! \param entropy The entropy [J/mol K]
  !! \param dsdt Temperature derivative [J/mol K]
  !! \param dsdp Pressure derivative [J/mol Pa]
  !! \param dsdz Composition derivative [J/mol^2]
  !!
  !! \author MH
  subroutine entropy_single(nc,comp,seos,T,P,Z,phase,residual,&
       entropy,dsdt,dsdp,dsdz)
    use eosdata
    use compdata, only: gendata_pointer
    use mbwr_additional, only: mbwr_sres, mbwr_volume
    use ideal, only: TP_Sideal_mix
    implicit none
    integer, intent(in) :: nc
    type (gendata_pointer), intent(in), dimension(nc) :: comp
    class(single_eos), intent(inout) :: seos
    real, dimension(nc), intent(in) :: Z
    real, intent(in) :: T, P
    logical, intent(in) :: residual
    real, intent(out) :: entropy
    real, optional, intent(out) :: dsdt, dsdp
    real, optional, dimension(nc), intent(out) :: dsdz
    integer, intent(in) :: phase
    ! Locals
    real :: s_id, dsdt_id, dsdp_id
    real, dimension(nc) :: dsdz_id
    real :: v
    !
    !-------- Specific for each equation of state ------------------------
    Choice_EoS: select case (seos%subeosidx) ! choose the Equation of State
    case (meosMbwr19, meosMbwr32)
      if (.not. residual) then
        call TP_Sideal_mix(nc, comp, T, P, Z, S_id, dsdt, dsdp, dsdz)
        if(present(dsdt)) dsdt_id = dsdt
        if(present(dsdp)) dsdp_id = dsdp
        if(present(dsdz)) dsdz_id = dsdz
      endif
      v = mbwr_volume (t,p,Z(1),phase,seos%mbwr_meos(1))
      call MBWR_Sres(seos%mbwr_meos(1), t, p, v, Z(1), entropy, dsdt, dsdp, dsdz)
      if (.not. residual) then
        entropy = entropy + s_id
        if(present(dsdt)) dsdt = dsdt + dsdt_id
        if(present(dsdp)) dsdp_id = dsdp + dsdp_id
        if(present(dsdz)) dsdz_id = dsdz + dsdz_id
      endif
    case (meosNist)
      call seos%nist(1)%meos%calc_entropy(t, p, Z, phase, entropy, &
           dsdt, dsdp, dsdz, residual=residual)
    end select Choice_EoS

  end subroutine entropy_single

  !---------------------------------------------------------------------- >
  !> This function calculates the Fugacity coefficient and the derivatives.
  !!
  !! \param T The temperature [K]
  !! \param P The pressure [Pa]
  !! \param Z The overall mole fraction [-]
  !! \param phase The phase, 1=liquid, 2=vapour
  !! \param dlnfdt Temperature derivative [1/K] (dlog(f)/dT)
  !! \param dlnfdp Pressure derivative [1/Pa]   (dlog(f)/dP)
  !! \param dlnfdz Composition derivative [1/mol] (dlog(f)/dNi)
  !! \param numder Analytical derivatives if true (default false)
  !! \param lnfug The fugacity coefficients [-]
  !!
  !! \author MH
  subroutine lnfugCoeff_single(nc,seos,T,P,Z,phase,lnfug,dlnfdt,dlnfdp,dlnfdz)
    use eosdata
    use mbwr_additional, only: mbwr_lnphi, mbwr_volume
    implicit none
    integer, intent(in) :: nc
    class(single_eos), intent(inout) :: seos
    real, dimension(nc), intent(in) :: Z
    real, intent(in) :: T, P
    integer, intent(in) :: phase
    real, dimension(nc), intent(out) :: lnfug
    real, optional, intent(out) :: dlnfdt(nc), dlnfdp(nc), dlnfdz(nc,nc)
    ! Locals
    real :: v

    !-------- Specific for each equation of state ------------------------
    Choice_EoS: select case (seos%subeosidx) ! choose the Equation of State
    case (meosMbwr19, meosMbwr32)
      v = mbwr_volume (t,p,Z(1),phase,seos%mbwr_meos(1))
      call MBWR_lnphi(seos%mbwr_meos(1), t, p, v, Z(1), lnfug, &
           dlnfdt, dlnfdp, dlnfdz)
    case (meosNist)
      call seos%nist(1)%meos%calc_lnphi(t, p, Z, phase, lnfug, &
           dlnfdt, dlnfdp, dlnfdz)
    end select Choice_EoS

  end subroutine lnfugCoeff_single

  !---------------------------------------------------------------------- >
  !> This subroutine calculates the residual gibbs energy and the derivatives.
  !!
  !! \param T The temperature [K]
  !! \param P The pressure [Pa]
  !! \param Z The overall mole fraction [-]
  !! \param phase The phase, 1=liquid, 2=vapour
  !! \param gr Residual gibbs energy [J/mol]
  !! \param dgrdt Temperature derivative [J/mol/K]
  !! \param dgrdp Pressure derivative [J/mol/Pa]
  !!
  !! \author Morten H.
  subroutine Gres_single(nc,seos,T,P,Z,phase,gr,dgrdt,dgrdp)
    use eosdata
    use mbwr_additional, only: MBWR_volume, MBWR_Gres
    implicit none
    integer, intent(in) :: nc
    class(single_eos), intent(inout) :: seos
    real, dimension(nc), intent(in) :: Z
    real, intent(in) :: T, P
    integer, intent(in) :: phase
    real, intent(out) :: gr
    real, optional, intent(out) :: dgrdt, dgrdp
    ! Locals
    real :: v

    !-------- Specific for each equation of state ------------------------
    Choice_EoS: select case (seos%subeosidx) ! choose the Equation of State
    case (meosMbwr19, meosMbwr32)
      v = mbwr_volume (t,p,Z(1),phase,seos%mbwr_meos(1))
      call MBWR_Gres(seos%mbwr_meos(1), t, p, v, Z(1), gr)
      if (present(dgrdt)) call stoperror("dgrdt not implemented in MBWR_Gres")
      if (present(dgrdp)) call stoperror("dgrdp not implemented in MBWR_Gres")
    case (meosNist)
      call seos%nist(1)%meos%calc_resgibbs(t, p, Z, phase, gr, dgrdt, dgrdp)
    end select Choice_EoS

  end subroutine Gres_Single

  !> Calculate pressure given composition, temperature and density
  !!
  !! \param T - Temprature [K]
  !! \param v - Specific volume [m3/mol]
  !!
  subroutine pressure_single(nc,seos,T,v,n,p,dpdv,dpdt,dpdz)
    use eosdata
    use mbwr_additional, only: MBWR_press
    implicit none
    integer, intent(in) :: nc
    class(single_eos), intent(inout) :: seos
    real, intent(in) :: t, v, n(nc)
    real, intent(inout) :: p
    real, optional, intent(inout) :: dpdv, dpdt
    real, dimension(nc), optional, intent(out) :: dpdz
    ! Locals
    real :: rho
    !---------------------------------------------------------------------
    rho = 1.0/v
    !-------- Specific for each equation of state ------------------------
    Choice_EoS: select case (seos%subeosidx) ! choose the Equation of State
    case (meosMbwr19, meosMbwr32)
      call MBWR_press(seos%mbwr_meos(1), T, v, sum(n), p, dpdv, dpdt)
    case (meosNist)
      call seos%nist(1)%meos%mp_pressure(rho=rho,t=T,p=p,p_rho=dpdv,p_T=dpdt)
      if (present(dpdv)) then
        dpdv = -dpdv/v**2
      end if
    end select Choice_EoS
    if (present(dpdz)) then
      call stoperror("single_component::pressure_single: dpdz not yet implemented")
    endif
  end subroutine pressure_single

  !> Calculate resudial reduced Helmholtz and differentials
  !!
  !! \param T - Temperature [K]
  !! \param v - Specific volume [m3/mol]
  !! \param n - Mole numbers [mol]
  subroutine Fres_single(nc,seos,T,V,n,F,F_T,F_V,F_n,F_TT,&
       F_TV,F_VV,F_Tn,F_Vn,F_nn)
    use eosdata
    use mbwr_additional, only: MBWR_Fres
    implicit none
    integer, intent(in) :: nc
    real, intent(in) :: T,V,n(nc)
    ! Output.
    class(single_eos), intent(inout) :: seos
    real, optional, intent(out) :: F,F_T,F_V,F_n(nc)
    real, optional, intent(out) :: F_TT,F_TV,F_Tn(nc),F_VV,F_Vn(nc),F_nn(nc,nc)
    ! Locals
    real :: v_l
    real, pointer :: F_n_p, F_nn_p, F_Tn_p, F_Vn_p
    real, target :: F_n_l, F_nn_l, F_Tn_l, F_Vn_l
    !---------------------------------------------------------------------
    if (present(F_n)) then
      F_n_p => F_n_l
    else
      F_n_p => NULL()
    endif
    if (present(F_nn)) then
      F_nn_p => F_nn_l
    else
      F_nn_p => NULL()
    endif
    if (present(F_Tn)) then
      F_Tn_p => F_Tn_l
    else
      F_Tn_p => NULL()
    endif
    if (present(F_Vn)) then
      F_Vn_p => F_Vn_l
    else
      F_Vn_p => NULL()
    endif
    v_l = V/n(1)

    !-------- Specific for each equation of state ------------------------
    Choice_EoS: select case (seos%subeosidx) ! choose the Equation of State
    case (meosMbwr19, meosMbwr32)
      call MBWR_Fres(seos%mbwr_meos(1), T, V, n(1), F, &
           F_T, F_v, F_TT, F_Tv, F_vv, F_n_p, F_Tn_p, F_vn_p, F_nn_p)
      if (present(F_Vn)) F_Vn = F_Vn_l
      if (present(F_n)) F_n = F_n_l
      if (present(F_nn)) F_nn = F_nn_l
      if (present(F_Tn)) F_Tn = F_Tn_l
    case (meosNist)
      call seos%nist(1)%meos%alphaResDerivs_Tv(T, v_l, F, &
           F_T, F_v, F_TT, F_Tv, F_vv, F_n_p, F_Tn_p, F_vn_p, F_nn_p)
      ! F = n*alpha(tau,delta)
      if (present(F)) F = F*n(1)
      if (present(F_T)) F_T = F_T*n(1)
      if (present(F_TT)) F_TT = F_TT*n(1)
      if (present(F_VV)) F_VV = F_VV/n(1)
      if (present(F_n)) then
        F_n = F_n_l
      endif
      if (present(F_nn)) then
        F_nn = F_nn_l/n(1)
      endif
      if (present(F_Tn)) then
        F_Tn = F_Tn_l
      endif
      if (present(F_Vn)) then
        F_Vn = F_Vn_l/n(1)
      endif
    end select Choice_EoS

  end subroutine Fres_single

  !> Calculate reduced ideal Helmholtz energy and differentials
  !!
  !! \param T - Temprature [K]
  !! \param v - Specific volume [m3/mol]
  !! \param n - Mole numbers [mol]
  subroutine Fid_single(nc,comp,seos,T,V,n,F,F_T,F_V,F_n,F_TT,&
       F_TV,F_VV,F_Tn,F_Vn,F_nn)
    use eosdata
    use compdata, only: gendata_pointer
    use ideal, only: Fideal_mix_SI
    implicit none
    integer, intent(in) :: nc
    type (gendata_pointer), intent(in), dimension(nc) :: comp
    real, intent(in) :: T,V,n(nc)
    ! Output.
    class(single_eos), intent(inout) :: seos
    real, optional, intent(out) :: F,F_T,F_V,F_n(nc)
    real, optional, intent(out) :: F_TT,F_TV,F_Tn(nc),F_VV,F_Vn(nc),F_nn(nc,nc)
    ! Locals
    real :: v_l
    real, pointer :: F_n_p, F_nn_p, F_Tn_p, F_Vn_p
    real, target :: F_n_l, F_nn_l, F_Tn_l, F_Vn_l
    !---------------------------------------------------------------------
    if ( seos%subeosidx == meosNist) then
      if (present(F_n)) then
        F_n_p => F_n_l
      else
        F_n_p => NULL()
      endif
      if (present(F_nn)) then
        F_nn_p => F_nn_l
      else
        F_nn_p => NULL()
      endif
      if (present(F_Tn)) then
        F_Tn_p => F_Tn_l
      else
        F_Tn_p => NULL()
      endif
      if (present(F_Vn)) then
        F_Vn_p => F_Vn_l
      else
        F_Vn_p => NULL()
      endif
    endif
    v_l = V/sum(n)

    !-------- Specific for each equation of state ------------------------
    Choice_EoS: select case (seos%subeosidx) ! choose the Equation of State
    case (meosMbwr19, meosMbwr32)
      call Fideal_mix_SI(nc, comp, T, V, n, Fid=F, Fid_T=F_T, Fid_v=F_V, &
           Fid_n=F_n, Fid_TT=F_TT, Fid_vv=F_vv, Fid_nn=F_nn, Fid_Tv=F_TV,&
           Fid_vn=F_Vn, Fid_Tn=F_Tn)
    case (meosNist)
      call seos%nist(1)%meos%alphaIdDerivs_Tv(T,v_l,alp=F,alp_T=F_T,&
           alp_v=F_V,alp_TT=F_TT,alp_Tv=F_TV,alp_vv=F_VV, &
           alp_n=F_n_p, alp_Tn=F_Tn_p, alp_vn=F_vn_p, alp_nn=F_nn_p)
      ! F = n*alpha(tau,delta)
      if (present(F)) F = F*n(1)
      if (present(F_T)) F_T = F_T*n(1)
      if (present(F_TT)) F_TT = F_TT*n(1)
      if (present(F_VV)) F_VV = F_VV/n(1)
      if (present(F_n)) then
        F_n = F_n_l
      endif
      if (present(F_nn)) then
        F_nn = F_nn_l/n(1)
      endif
      if (present(F_Tn)) then
        F_Tn = F_Tn_l
      endif
      if (present(F_Vn)) then
        F_Vn = F_Vn_l/n(1)
      endif
    end select Choice_EoS

  end subroutine Fid_single

end module single_component

max_stars_repo_name: SINTEF/Thermopack

max_issues_repo_name: SINTEF/Thermopack

max_forks_repo_name: SINTEF/Thermopack


text:
      SUBROUTINE SEARCH(XPARAM,ALPHA,SIG,NVAR,GMIN,OKC,OKF, FUNCT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION XPARAM(*), SIG(*)
************************************************************************
*
* SEARCH PERFORMS A LINE SEARCH FOR POWSQ. IT MINIMIZES THE NORM OF
*        THE GRADIENT VECTOR IN THE DIRECTION SIG.
*
* ON INPUT  XPARAM = CURRENT POINT IN NVAR DIMENSIONAL SPACE.
*           ALPHA  = STEP SIZE (IN FACT ALPHA IS CALCULATED IN SEARCH).
*           SIG    = SEARCH DIRECTION VECTOR.
*           NVAR   = NUMBER OF PARAMETERS IN SIG (& XPARAM)
*
* ON OUTPUT XPARAM = PARAMETERS OF MINIMUM.
*           ALPHA  = DISTANCE TO MINIMUM.
*           GMIN   = GRADIENT NORM AT MINIMUM.
*           OKC    = EXITED BEFORE COUNTS WERE EXCEEDED (LOGICAL)
*           OKF    = FUNCTION WAS IMPROVED.
************************************************************************
      COMMON /SIGMA1/ GNEXT, AMIN, ANEXT
      COMMON /SIGMA2/  GNEXT1(MAXPAR), GMIN1(MAXPAR)
      COMMON/KEYWRD/ KEYWRD
      DIMENSION GRAD(MAXPAR),XREF(MAXPAR), GREF(MAXPAR), XMIN1(MAXPAR)
      CHARACTER*80 KEYWRD
      LOGICAL FIRST, DEBUG, OKC, OKF, NOPR
      DATA  FIRST /.TRUE./
      IF( FIRST ) THEN
         FIRST = .FALSE.
C
C    TOLG   = CRITERION FOR EXIT BY RELATIVE CHANGE IN GRADIENT.
C
         DEBUG=(INDEX(KEYWRD,'LINMIN') .NE. 0)
         NOPR=( .NOT. DEBUG)
         LOOKS=0
         OKF=.TRUE.
         TINY=0.1D0
         TOLERG=0.02D0
         IF(NVAR .EQ. 1)TIMY=0.D0
         G=100.D0
         XMAXM=2.D0
         ALPHA=0.1D0
      ENDIF
      ANEXT1=0.D0
      DO 10 I=1,NVAR
         GREF(I)  =GMIN1(I)
         GNEXT1(I)=GMIN1(I)
         XMIN1(I) =XPARAM(I)
   10 XREF(I)  =XPARAM(I)
      IF(ABS(ALPHA) .GT. 0.2)ALPHA=SIGN(0.2D0,ALPHA)
      IF(DEBUG) THEN
         WRITE(6,'('' SEARCH DIRECTION VECTOR'')')
         WRITE(6,'(6F12.6)')(SIG(I),I=1,NVAR)
         WRITE(6,'('' INITIAL GRADIENT VECTOR'')')
         WRITE(6,'(6F12.6)')(GMIN1(I),I=1,NVAR)
      ENDIF
      GB=DOT(GMIN1,GREF,NVAR)
      IF(DEBUG) WRITE(6,'('' GRADIENT AT START OF SEARCH:'',F16.6)')
     1SQRT(GB)
      GSTORE=GB
      AMIN=0.D0
      GMINN=1.D9
C
C
      TA=0.D0
      GA=GB
      GB=1.D9
      ITRYS=0
      GOTO 30
   20 SUM=GA/(GA-GB)
      ITRYS=ITRYS+1
      IF(ABS(SUM) .GT. 3.D0) SUM=SIGN(3.D0,SUM)
      ALPHA=(TB-TA)*SUM+TA
C
C         XPARAM IS THE GEOMETRY OF THE PREDICTED MINIMUM ALONG THE LINE
C
   30 CONTINUE
      DO 40 I=1,NVAR
   40 XPARAM(I)=XREF(I)+ALPHA*SIG(I)
C
C         CALCULATE GRADIENT NORM AND GRADIENTS AT THE PREDICTED MINIMUM
C
      IF(ITRYS.EQ.1)THEN
      DO 50 I=1,NVAR
   50 GRAD(I)=0.D0
      ENDIF
      CALL COMPFG (XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
      LOOKS=LOOKS+1
C
C          G IS THE PROJECTION OF THE GRADIENT ALONG SIG.
C
      G=DOT(GREF,GRAD,NVAR)
      GTOT=SQRT(DOT(GRAD,GRAD,NVAR))
      IF( .NOT. NOPR)
     1WRITE(6,'('' LOOKS'',I3,'' ALPHA ='',F12.6,'' GRADIENT'',F12.3,
     2'' G  ='',F16.6)')
     3LOOKS,ALPHA,SQRT(DOT(GRAD,GRAD,NVAR)),G
      IF(GTOT .LT. GMINN) THEN
         GMINN=GTOT
         IF(ABS(AMIN-ALPHA) .GT.1.D-2) THEN
*
* WE CAN MOVE ANEXT TO A POINT NEAR, BUT NOT TOO NEAR, AMIN, SO THAT THE
* SECOND DERIVATIVESWILLBEREALISTIC(D2E/DX2=(GNEXT1-GMIN1)/(ANEXT-AMIN))
*
            ANEXT=AMIN
            DO 60 I=1,NVAR
   60       GNEXT1(I)=GMIN1(I)
         ENDIF
         AMIN=ALPHA
         DO 70 I=1,NVAR
            IF(GMINN.LT.GMIN) XMIN1(I)=XPARAM(I)
   70    GMIN1(I)=GRAD(I)
         IF(GMIN.GT.GMINN)GMIN=GMINN
      ENDIF
      IF(ITRYS .GT. 8) GOTO 80
      IF (ABS(G/GSTORE).LT.TINY .OR. ABS(G) .LT. TOLERG) GO TO 80
      IF(ABS(G) .LT. MAX(ABS(GA),ABS(GB)) .OR.
     1     GA*GB .GT. 0.D0 .AND. G*GA .LT. 0.D0) THEN
C
C   G IS AN IMPROVEMENT ON GA OR GB.
C
         IF(ABS(GB) .LT. ABS(GA))THEN
            TA=ALPHA
            GA=G
            GO TO 20
         ELSE
            TB=ALPHA
            GB=G
            GO TO 20
         ENDIF
      ELSE
C#         WRITE(6,'(//10X,'' FAILED IN SEARCH, SEARCH CONTINUING'')')
         GOTO 80
      ENDIF
   80 LNSTOP=4
      GMINN=SQRT(DOT(GMIN1,GMIN1,NVAR))
      DO 90 I=1,NVAR
   90 XPARAM(I)=XMIN1(I)
      IF(DEBUG) THEN
         WRITE(6,'('' AT EXIT FROM SEARCH'')')
         WRITE(6,'('' XPARAM'',6F12.6)')(XPARAM(I),I=1,NVAR)
         WRITE(6,'('' GNEXT1'',6F12.6)')(GNEXT1(I),I=1,NVAR)
         WRITE(6,'('' GMIN1 '',6F12.6)')(GMIN1(I),I=1,NVAR)
         WRITE(6,'('' AMIN, ANEXT,GMIN'',4F12.6)')
     1    AMIN,ANEXT,GMIN
      ENDIF
      IF(GMINN.GT.GMIN)THEN
      DO 45 I=1,NVAR
  45  XPARAM(I)=XREF(I)
      ENDIF
      RETURN
C
      END

ext: f

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
Wilson Sauthoff hails from Oxnard Beach, California.

Wilson is a Environmental Science and Management major emphasizing Climate Change and Air Quality at UC Davis.

Wilson is a current ASUCD Vice Chair of the Environmental Policy and Planning Commission.

Wilson is the Bodega Marine Laboratory peer advisor and recruiter.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!----------------------------------------------------------------------
!     Redoing the local bifurcation of the model analysed by
!     Martin Boer, ch. 4 of thesis.
!     Cycle-to-cycle connection 3D food chain
!
!     Homoclinic, but set up also for heteroclinic
!
!     U(1-3) = limit cycle
!     U(4-6) = eigenfunction (ef) stable
!     U(7-9) = limit cycle again ("second cycle")
!     U(10-12) = ef unstable
!     U(13-15) = connecting orbit
!
!     George van Voorn, June 20th 2007
!     Bart Oldeman, December 2008
!----------------------------------------------------------------------

  SUBROUTINE RHS(N,U,PAR,F,DFDU,IJAC)
! ---------- ---

    IMPLICIT NONE
    INTEGER N
    DOUBLE PRECISION U(N),PAR(*),F(N),DFDU(N,N)
    LOGICAL IJAC
    DOUBLE PRECISION f1,f2
    DOUBLE PRECISION d1,d2,x,y,z,a1,a2,b1,b2
!     Declaring variables
    x=U(1)
    y=U(2)
    z=U(3)
!     Parameter values
    d1=PAR(1)
    d2=PAR(2)

    a1 = 5.d0
    a2 = 0.1d0
    b1 = 3.d0
    b2 = 2.d0

    f1=a1/(1d0 + b1*x)
    f2=a2/(1d0 + b2*y)
!     Formulas
    F(1)= x*(1d0 - x) - f1*x*y
    F(2)= f1*x*y - d1*y - f2*y*z
    F(3)= f2*y*z - d2*z

    IF(.NOT.IJAC)RETURN

!     Jacobian elements
    DFDU(1,1)=1d0-2*x-f1*y+b1*x*y*(f1**2)/a1
    DFDU(1,2)=-f1*x
    DFDU(1,3)=0.0d0
    DFDU(2,1)=f1*y-b1*x*y*(f1**2)/a1
    DFDU(2,2)=f1*x-d1-f2*z+b2*y*z*(f2**2)/a2
    DFDU(2,3)=-f2*y
    DFDU(3,1)=0.0d0
    DFDU(3,2)=f2*z-b2*y*z*(f2**2)/a2
    DFDU(3,3)=f2*y-d2

  END SUBROUTINE RHS

  SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
! ---------- ---

    IMPLICIT NONE
    INTEGER NDIM,ICP(*),IJAC
    DOUBLE PRECISION U(NDIM),PAR(*),F(NDIM),DFDU(NDIM,*),DFDP(NDIM,*)
    DOUBLE PRECISION A(3,3)

    CALL RHS(3,U,PAR,F,A,NDIM/=3)
    IF(NDIM==3)RETURN

    F(1:3) = PAR(11)*F(1:3)

    IF(NDIM==6)THEN
!     Variational equations
!     PAR(11) = cycle period
!     PAR(12) = log(FM)
       F(4:6) = PAR(11)*MATMUL(A(:,:),U(4:6))-PAR(12)*U(4:6)
       RETURN
    ENDIF

!     Adjoint variational equations
!     PAR(11) = cycle period
!     PAR(4) = log(FM)
    F(4:6) = -PAR(11)*MATMUL(TRANSPOSE(A(:,:)),U(4:6))-PAR(4)*U(4:6)

!     Functions limit cycle doubled U(7-9)
!     We pretend this is a second different cycle
    CALL RHS(3,U(7:9),PAR,F(7:9),A,.TRUE.)
    F(7:9) = PAR(6)*F(7:9)
!     Period=PAR(6), obviously equal to PAR(11), U(10-12)
!     PAR(5) = log(FM)
    F(10:12) = -PAR(6)*MATMUL(TRANSPOSE(A(:,:)),U(10:12))-PAR(5)*U(10:12)

    IF(NDIM==12)RETURN

!     Connection rescaled, PAR(7) = Tc
    CALL RHS(3,U(13:15),PAR,F(13:15),A,.FALSE.)
    F(13:15) = PAR(7)*F(13:15)

  END SUBROUTINE FUNC

  SUBROUTINE STPNT(NDIM,U,PAR)
! ---------- ---

    IMPLICIT NONE
    INTEGER NDIM
    DOUBLE PRECISION, INTENT(INOUT) :: U(NDIM), PAR(*)

    DOUBLE PRECISION FLOQ,epsilon

    IF(NDIM==12)THEN
       ! extending from NDIM=6 to NDIM=12
       U(7:12)=U(1:6)
       RETURN
    ELSEIF(NDIM==15)THEN
       ! extending from NDIM=12 to NDIM=15
       epsilon=PAR(13)
       U(13:15)=PAR(14:16) + epsilon*PAR(17:19)
       RETURN
    ENDIF

!     0.16<par1=d1<0.32,0.0075<par2=d2<0.015

    PAR(1)=0.5d0
    PAR(2)=0.0125d0
    U(1)=0.7415816238d0
    U(2)=0.1666666667d0
    U(3)=8.664398854d0

!     Guess for FM
    FLOQ=-10d0
    PAR(4)=FLOQ
    PAR(5)=FLOQ
    PAR(12)=FLOQ
    PAR(13)=0.01d0 ! epsilon

!     Eigenvector norms
    PAR(10)=0.

    PAR(25)=0.
    PAR(26)=0.

  END SUBROUTINE STPNT

  SUBROUTINE BCND(NDIM,PAR,ICP,NBC,U0,U1,FB,IJAC,DBC)
! ---------- ----

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NDIM, ICP(*), NBC, IJAC
    DOUBLE PRECISION, INTENT(IN) :: PAR(*), U0(NDIM), U1(NDIM)
    DOUBLE PRECISION, INTENT(OUT) :: FB(NBC)
    DOUBLE PRECISION, INTENT(INOUT) :: DBC(NBC,*)

    DOUBLE PRECISION g1(3),h1(3),h2(3),A(3,3)
    DOUBLE PRECISION epsilon

!     Periodicity boundary conditions
!     -------------------------------

    FB(1:3) = U0(1:3) - U1(1:3)
    FB(4:6) = U1(4:6) - U0(4:6)

    IF(NBC==7)THEN
       FB(7) = DOT_PRODUCT(U0(4:6),U0(4:6)) - PAR(10)
       RETURN
    ENDIF

    FB(7) = DOT_PRODUCT(U0(4:6),U0(4:6)) - PAR(25)

    ! Double

    FB(8:10) = U0(7:9) - U1(7:9)
    FB(11:13) = U1(10:12) - U0(10:12)

    FB(14) = DOT_PRODUCT(U0(10:12),U0(10:12)) - PAR(26)

    IF(NBC==14)RETURN

    ! Cycle 1 end-point connection
    CALL RHS(3,U0(1:3),PAR,g1,A,.FALSE.)
    ! Displacement from the cycle; end connection
    h1(1:3)=U1(13:15)-U0(1:3)

    ! Eq.(12g)/(17h) - h11
    FB(15) = DOT_PRODUCT(h1(1:3),g1(1:3)) - PAR(21)

    ! cycle 2 starting point connection
    ! For coherence we use U(7-9)
    CALL RHS(3,U0(7:9),PAR,g1,A,.FALSE.)
    ! Displacement from the cycle; start connection
    h2(1:3)=U0(13:15)-U0(7:9)

    ! Eq.(12h)/(17i) - h12
    FB(16) = DOT_PRODUCT(h2(1:3),g1(1:3)) - PAR(22)
    ! Eq.(12i)/(17j) - h21
    FB(17) = DOT_PRODUCT(U0(4:6),h1(1:3)) - PAR(23)
    ! Eq.(12j)/(17k) - h22
    FB(18) = DOT_PRODUCT(U0(10:12),h2(1:3)) - PAR(24)
    ! Eq.(12k)
    FB(19) = DOT_PRODUCT(h2(1:3),h2(1:3)) - PAR(8)
  END SUBROUTINE BCND

  SUBROUTINE ICND(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,FI,IJAC,DINT)
! ---------- ----

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NDIM, ICP(*), NINT, IJAC
    DOUBLE PRECISION, INTENT(IN) :: PAR(*)
    DOUBLE PRECISION, INTENT(IN) :: U(NDIM),UOLD(NDIM),UDOT(NDIM),UPOLD(NDIM)
    DOUBLE PRECISION, INTENT(OUT) :: FI(NINT)
    DOUBLE PRECISION, INTENT(INOUT) :: DINT(NINT,*)

!   Integral phase condition
!   ------------------------
    FI(1) = U(1)*UPOLD(1)+U(2)*UPOLD(2)+U(3)*UPOLD(3)
    IF(NINT==1)RETURN
    FI(2) = U(7)*UPOLD(7)+U(8)*UPOLD(8)+U(9)*UPOLD(9)

  END SUBROUTINE ICND

  SUBROUTINE FOPT
  END SUBROUTINE FOPT

  SUBROUTINE PVLS(NDIM,U,PAR)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NDIM
    DOUBLE PRECISION, INTENT(IN) :: U(NDIM)
    DOUBLE PRECISION, INTENT(INOUT) :: PAR(*)
    DOUBLE PRECISION GETP

    DOUBLE PRECISION u0(12),g1(3),g2(3),h1(3),h2(3),A(3,3),b(3),epsilon
    INTEGER I

    IF(NDIM==6)THEN
       DO I=1,6
          PAR(13+I)=GETP('BV0',I,U)
       ENDDO
       PAR(6)=PAR(11)
    ELSEIF(NDIM==12)THEN
      DO I=1,12
         u0(I) = GETP('BV0',I,U)
      ENDDO

      ! f' at base point cycle 1
      CALL RHS(3,u0(1:3),PAR,g1,A,.FALSE.)
      ! f' at base point cycle 2
      CALL RHS(3,u0(7:9),PAR,g2,A,.FALSE.)
      ! start & end of initial point
      epsilon=PAR(13)
      b(1:3) = PAR(14:16)+epsilon*PAR(17:19)

      ! Displacement from the cycle at end connection
      h1(1:3)=b(1:3)-u0(1:3)
      ! Displacement from the cycle at start connection
      h2(1:3)=b(1:3)-u0(7:9)

!     Homotopy parameter values at start continuation
      ! Eq.(12h)/(17i)
      PAR(21) = DOT_PRODUCT(h1(1:3),g1(1:3))
      ! Eq.(12i)/(17j)
      PAR(22) = DOT_PRODUCT(h2(1:3),g2(1:3))
      ! Eq.(12j)/(17k)
      PAR(23) = DOT_PRODUCT(u0(4:6),h1(1:3))
      ! Eq.(12k)/(17l)
      PAR(24) = DOT_PRODUCT(u0(10:12),h2(1:3))

      ! for BC(19)/Eq.(12k)
      PAR(8) = DOT_PRODUCT(h2(1:3),h2(1:3))
   ELSEIF(NDIM==15)THEN
      PAR(27)=GETP('MIN',3,U)
      PAR(28)=GETP('MIN',15,U)
   ENDIF

  END SUBROUTINE PVLS


text:
      DOUBLE PRECISION FUNCTION DCABS1( Z )
*
*  -- Automatically Tuned Linear Algebra Software (ATLAS)
*     (C) Copyright 2000 All Rights Reserved
*
*  -- ATLAS routine -- F77 Interface -- Version 3.9.24 -- December 25, 2000
*
*  Author         : Antoine P. Petitet
*  Originally developed at the University of Tennessee,
*  Innovative Computing Laboratory,  Knoxville TN, 37996-1301, USA.
*
*  ---------------------------------------------------------------------
*
*  -- Copyright notice and Licensing terms:
*
*  Redistribution  and  use in  source and binary forms, with or without
*  modification, are  permitted provided  that the following  conditions
*  are met:
*
*  1. Redistributions  of  source  code  must retain the above copyright
*     notice, this list of conditions and the following disclaimer.
*  2. Redistributions in binary form must reproduce  the above copyright
*     notice,  this list of conditions, and the  following disclaimer in
*     the documentation and/or other materials provided with the distri-
*     bution.
*  3. The name of the University,  the ATLAS group,  or the names of its
*     contributors  may not be used to endorse or promote products deri-
*     ved from this software without specific written permission.
*
*  -- Disclaimer:
*
*  THIS  SOFTWARE  IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
*  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,  INCLUDING,  BUT NOT
*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
*  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
*  OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,  INDIRECT, INCIDENTAL, SPE-
*  CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
*  TO,  PROCUREMENT  OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
*  OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
*  RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT  (IN-
*  CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
*  THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*  ---------------------------------------------------------------------
*
*     .. Scalar Arguments ..
      DOUBLE COMPLEX     Z
*     ..
*
*  Purpose
*  =======
*
*  DCABS1 returns |real( z )| + |imag( z )|  where z is a double complex
*  number.
*
*  Arguments
*  =========
*
*  Z       (input)                       COMPLEX*16
*          On entry, the complex number z.
*
*  =====================================================================
*
*     .. Local Scalars ..
      DOUBLE COMPLEX     ZZ
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   T( 2 )
*     ..
*     .. Equivalence Statements ..
      EQUIVALENCE        ( ZZ, T( 1 ) )
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          DABS
*     ..
*     .. Executable Statements ..
*
      ZZ = Z
      DCABS1 = DABS( T( 1 ) ) + DABS( T( 2 ) )
*
      RETURN
*
*     End of DCABS1
*
      END

max_stars_repo_name: kevleyski/math-atlas

max_issues_repo_name: kevleyski/math-atlas

max_forks_repo_name: kevleyski/math-atlas


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C $Log: vector.f,v $
C Revision 1.3  2009/03/25 12:36:49  gdsjaar
C Add copyright and license notice to all files.
C Permission to assert copyright has been granted; blot is now open source, BSD
C
C Revision 1.2  1994/07/21 15:28:07  gdsjaar
C Moved more commons into includes.
C
c Revision 1.1  1994/04/07  20:17:31  gdsjaar
c Initial checkin of ACCESS/graphics/blotII2
c
c Revision 1.3  1991/06/25  16:10:02  gdsjaar
c Fixed? problem with calls to ugrcol -- changed
c call ugrcol(idelb(ielb),...) to call ugrcol(ielb,...)
c
c Revision 1.2  1990/12/14  08:59:31  gdsjaar
c Added RCS Id and Log to all files
c
C=======================================================================
      SUBROUTINE VECTOR (ISYTYP, VAR, NUM, LENF, NLNKF, HIDENE,
     &   XNE, YNE, ZNE, IN2ELB, ISVOK, VECMAX, BLKCOL, IDELB, *)
C=======================================================================

C   --*** VECTOR *** (DETOUR) Plot nodal or element vector
C   --   Written by Amy Gilkey - revised 03/10/88
C   --   D. P. Flanagan, 12/08/83
C   --
C   --VECTOR draws a vector representing 2 or 3 variables for each node
C   --or element.  It processes each node or element by element block.
C   --Only nodes or elements of selected element blocks are drawn.
C   --
C   --Parameters:
C   --   ISYTYP - IN - the vector type (as in MODTYP of /DETOPT/)
C   --   VAR - IN - the vector variable values
C   --   NUM - IN - the number of nodes (if nodal) or faces (if element)
C   --   LENF - IN - the cumulative face counts by element block
C   --      (only if element)
C   --   NLNKF - IN - the number of nodes per face (only if element)
C   --   HIDENE(i) - IN - true iff node i or face i is hidden (3D only)
C   --   XNE, YNE, ZNE - IN - the nodal coordinates (if nodal)
C   --      or the face center coordinates (if element)
C   --   IN2ELB - IN - the element block for each node; (only if nodal)
C   --      <0 if not in any selected element block
C   --      =0 if in more than one selected element block
C   --   ISVOK - IN - ISVOK(i) is true iff the vector variables are defined
C   --      for element block i (only if element)
C   --   VECMAX - IN - the maximum vector variable value, scaled
C   --   BLKCOL - IN - the user selected colors of the element blocks.
C   --                    BLKCOL(0) = 1 if the user defined material
C   --                                colors should be used in mesh plots.
C   --                              = -1 if program selected colors should
C   --                                be used.
C   --                    BLKCOL(i) = the user selected color of element
C   --                               block i:
C   --                                  -2 - no color selected by user.
C   --                                  -1 - black
C   --                                   0 - white
C   --                                   1 - red
C   --                                   2 - green
C   --                                   3 - yellow
C   --                                   4 - blue
C   --                                   5 - cyan
C   --                                   6 - magenta
C   --   * - return statement if cancel function active
C   --   * - return statement if the cancel function is active
C   --
C   --Common Variables:
C   --   Uses NELBLK of /DBNUMS/
C   --   Uses IS3DIM, NUMNPF of /D3NUMS/
C   --   Uses VECSCL of /ETCOPT/
C   --   Uses ROTMAT of /ROTOPT/
C   --   Uses DTW, VWSCL of /DEVDAT/

      include 'dbnums.blk'
      include 'd3nums.blk'
      include 'etcopt.blk'
      include 'rotopt.blk'
      include 'devdat.blk'

      CHARACTER*(*) ISYTYP
      REAL VAR(NUM,*)
      INTEGER LENF(0:NELBLK)
      INTEGER NLNKF(NELBLK)
      LOGICAL HIDENE(NUM)
      REAL XNE(NUM), YNE(NUM), ZNE(NUM)
      INTEGER IN2ELB(NUMNPF)
      LOGICAL ISVOK(NELBLK)
      REAL VECMAX
      INTEGER BLKCOL(0:NELBLK)
      INTEGER IDELB(*)

      LOGICAL GRABRT
      REAL ZERO3(3)
      SAVE ZERO3

      DATA ZERO3 / 0.0, 0.0, 0.0 /

      IF (VECSCL .EQ. 0.0) THEN
         VSCL = 1.0
      ELSE
         VSCL = VECSCL * 0.05 * VECMAX
      END IF
      ASCL = VWSCL

      IF (ISYTYP .EQ. 'NODE') THEN

         DO 110 IELB = 0, NELBLK

C         --Set the vector color
c            IF (IELB .GT. 0) THEN
c               ITEMP = IDELB(IELB)
c            ELSE
               ITEMP = IELB
c            ENDIF
            CALL UGRCOL (ITEMP, BLKCOL)

            DO 100 INP = 1, NUMNPF

               IF (IS3DIM) THEN
                  IF (HIDENE(INP)) GOTO 100
               END IF

               IF (IN2ELB(INP) .EQ. IELB) THEN
                  IF (GRABRT ()) RETURN 1

C               --Call vector routine

                  IF (.NOT. IS3DIM) THEN
                     CALL VEC (IS3DIM, XNE(INP), YNE(INP), 0.0,
     &                  VAR(INP,1), VAR(INP,2), 0.0, VSCL, ASCL)
                  ELSE
                     XVAR = VAR(INP,1)
                     YVAR = VAR(INP,2)
                     ZVAR = VAR(INP,3)
                     CALL BL_ROTATE (1, 1, ROTMAT, ZERO3,
     &                  XVAR, YVAR, ZVAR, XVAR, YVAR, ZVAR)
                     CALL VEC (IS3DIM, XNE(INP), YNE(INP), ZNE(INP),
     &                  XVAR, YVAR, ZVAR, VSCL, ASCL)
                  END IF
               END IF

  100       CONTINUE

            CALL PLTFLU
  110    CONTINUE

      ELSE IF (ISYTYP .EQ. 'ELEMENT') THEN

         DO 130 IELB = 1, NELBLK
            IF (ISVOK(IELB)) THEN

C            --Set the vector color
c               CALL UGRCOL (IDELB(IELB), BLKCOL)
               CALL UGRCOL (IELB, BLKCOL)

               DO 120 IFAC = LENF(IELB-1)+1, LENF(IELB)
                  IF (IS3DIM) THEN
                     IF (HIDENE(IFAC)) GOTO 120
                  END IF

C               --Call vector routine

                  IF (GRABRT ()) RETURN 1

                  IF (.NOT. IS3DIM) THEN
                     CALL VEC (IS3DIM, XNE(IFAC), YNE(IFAC), 0.0,
     &                  VAR(IFAC,1), VAR(IFAC,2), 0.0, VSCL, ASCL)
                  ELSE
                     XVAR = VAR(IFAC,1)
                     YVAR = VAR(IFAC,2)
                     ZVAR = VAR(IFAC,3)
                     CALL BL_ROTATE (1, 1, ROTMAT, ZERO3,
     &                  XVAR, YVAR, ZVAR, XVAR, YVAR, ZVAR)
                     CALL VEC (IS3DIM, XNE(IFAC), YNE(IFAC), ZNE(IFAC),
     &                  XVAR, YVAR, ZVAR, VSCL, ASCL)
                  END IF

  120          CONTINUE

               CALL PLTFLU
            END IF
  130    CONTINUE

      ELSE IF ((ISYTYP .EQ. 'SIGMAX') .OR. (ISYTYP .EQ. 'SIGMIN')) THEN

         DEG180 = 2.0 * ATAN (1.0)

         DO 150 IELB = 1, NELBLK
            IF (ISVOK(IELB)) THEN

C            --Set the vector color
c               CALL UGRCOL (IDELB(IELB), BLKCOL)
               CALL UGRCOL (IELB, BLKCOL)

               DO 140 IFAC = LENF(IELB-1)+1, LENF(IELB)
                  IF (IS3DIM) THEN
                     IF (HIDENE(IFAC)) GOTO 140
                  END IF

                  SIGXX = VAR(IFAC,1)
                  SIGYY = VAR(IFAC,2)
                  TAUXY = VAR(IFAC,3)
                  Q1 = 0.5 * (SIGXX + SIGYY)
                  Q2 = 0.5 * (SIGXX - SIGYY)
                  Q3 = SQRT (Q2 * Q2 + TAUXY * TAUXY)
                  SIGMAX = Q1 + Q3
                  SIGMIN = Q1 - Q3
                  IF ((TAUXY .EQ. 0.0) .AND. (SIGMAX-SIGYY .EQ. 0.0))
     &               GOTO 140
                  THETA = ATAN2 (TAUXY, SIGMAX - SIGYY)
                  IF (ISYTYP .EQ. 'SIGMAX') THEN
                     XSIG = SIGMAX * COS (THETA)
                     YSIG = SIGMAX * SIN (THETA)
                  ELSE
                     XSIG = SIGMIN * COS (THETA + DEG180)
                     YSIG = SIGMIN * SIN (THETA + DEG180)
                  END IF

C               --Call vector routine

                  IF (GRABRT ()) RETURN 1

                  IF (.NOT. IS3DIM) THEN
                     CALL VEC (IS3DIM, XNE(IFAC), YNE(IFAC), 0.0,
     &                  XSIG, YSIG, 0.0, VSCL, ASCL)
                  ELSE
                     ZSIG = 0.0
                     CALL BL_ROTATE (1, 1, ROTMAT, ZERO3,
     &                  XSIG, YSIG, ZSIG, XSIG, YSIG, ZSIG)
                     CALL VEC (IS3DIM, XNE(IFAC), YNE(IFAC), ZNE(IFAC),
     &                  XSIG, YSIG, ZNE(IFAC), VSCL, ASCL)
                  END IF

  140          CONTINUE

               CALL PLTFLU
            END IF
  150    CONTINUE
      END IF

      RETURN
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


SUBROUTINE DGGSVD_F95( A, B, ALPHA, BETA, K, L, U, V, Q, IWORK, INFO )
!
!  -- LAPACK95 interface driver routine (version 3.0) --
!     UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK
!     September, 2000
!
!  .. USE STATEMENTS ..
   USE LA_PRECISION, ONLY: WP => DP
   USE LA_AUXMOD, ONLY: ERINFO
   USE F77_LAPACK, ONLY: GGSVD_F77 => LA_GGSVD
!  .. IMPLICIT STATEMENT ..
   IMPLICIT NONE
!  .. SCALAR ARGUMENTS ..
   INTEGER, INTENT(OUT), OPTIONAL :: INFO, K, L
!  .. ARRAY ARGUMENTS ..
   REAL(WP), INTENT(INOUT) :: A(:,:), B(:,:)
   REAL(WP), INTENT(OUT) :: ALPHA(:), BETA(:)
   REAL(WP), INTENT(OUT), OPTIONAL, TARGET :: U(:,:), V(:,:), Q(:,:)
   INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IWORK(:)
!----------------------------------------------------------------------
!
! Purpose
! =======
!
!      LA_GGSVD computes the generalized singular values and, optionally,
! the transformation matrices from the generalized singular value
! decomposition (GSVD) of a real or complex matrix pair (A,B), where A
! is m by n and B is p by n. The GSVD of (A,B) is written
!       A = U * SIGMA1(0, R)*Q^H , B = V * SIGMA2(0, R)*Q^H
! where U , V and Q are orthogonal (unitary) matrices of dimensions m by m,
! p by p and n by n, respectively. Let l be the rank of B and r the rank of
! the (m + p) * n matrix ( A )
!                        ( B )
! , and let k = r-l. Then SIGMA1 and SIGMA2 are m*(k + l) and p * (k + l)
! "diagonal" matrices, respectively, and R is a (k + l) * (k + l)
! nonsingular triangular matrix. The detailed structure of SIGMA1 ,SIGMA2
! and R depends on the sign of (m - k - l) as follows:
!       The case m-k-l>=0:
!
!                               k   l
!                      k      ( I   0 )
!        SIGMA1 =      l      ( 0   C )
!                    m-k-l    ( 0   0 )
!
!
!                           k   l
! 	SIGMA2 =    l   ( 0   S )
! 	          p - l ( 0   S )
!
!
!                          n-k-l   k     l
!          (0, R) =   k   (  0    R11   R12  )
!  	              l   (  0     0    R22  )
!
! where C^2 + S^2 = I . We define
! alpha(1)=alpha(2)=...=alpha(k) = 1, alpha(k+i)=c(i i), i=1,2,...,l
! beta(1) = beta(2) = ... = beta(k) = 0, beta(k+i) = s(i i), i=1,2,...,l
!
! The case m-k-l < 0:
!
!                                k    m-k    k+l-m
!            SIGMA1 =    k    (  I     0       0   )
!   	                m-k   (  0     C       0   )
!
!
! 		               k    m-k     k+l-m
!                        m-k   ( 0     S        0  )
!            SIGMA2 =   k+l-m  ( 0     0        I  )
!  	                 p-l   ( 0     0        0  )
!
!
!                                  n-k-l    k     m-k   k+l-m
! 		  	  k      (   0     R11    R12    R13  )
!              (0,R) =   m-k     (   0      0     R22    R23  )
! 	                k+l-m    (   0      0      0     R33  )
!
! where C^2 + S^2 = I . We define
!  alpha(1)=alpha(2)=...=alpha(k) = 1, alpha(k+i)=c(i i), i =1,2,...,m-k,
!  alpha(m+1) = alpha(m+2)=...= alpha(k+l) = 0
!  beta(1)=beta(2)= ... =beta(k)=0, beta(k+i)=s(i i), i=1,2,...,m-k,
!  beta(m+1) = beta(m+2) = ... = beta(k+l) = 1
!
! In both cases the generalized singular values of the pair (A,B) are the
! ratios
!  sigma(i) = alpha(i)/beta(i), i = 1,2, ... ,k+l
!
! The first k singular values are infinite. The finite singular values
! are real and nonnegative.
!     LA_GGSVD computes the real (nonnegative) scalars alpha(i), beta(i),
! i=1,2,..., k+l , the matrix R, and, optionally, the transformation
! matrices U , V and Q.
!
! =========
!
!      SUBROUTINE LA_GGSVD( A, B, ALPHA, BETA, K=k, L=l, &
!                       U=u, V=v, Q=q, IWORK=iwork, INFO=info )
!          <type>(<wp>), INTENT(INOUT) :: A(:,:), B(:,:)
!          REAL(<wp>), INTENT(OUT) :: ALPHA(:), BETA(:)
!          INTEGER, INTENT(OUT), OPTIONAL :: K, L
!          <type>(<wp>), INTENT(OUT), OPTIONAL :: U(:,:), V(:,:), Q(:,:)
!          INTEGER, INTENT(IN), OPTIONAL :: IWORK(:)
!          INTEGER, INTENT(OUT), OPTIONAL :: INFO
!       where
!          <type> ::= REAL | COMPLEX
!          <wp>   ::= KIND(1.0) | KIND(1.0D0)
!
! Arguments
! =========
!
! A      (input/output) REAL or COMPLEX array, shape (:,:) with
!        size(A,1) = m and size(A,2) = n.
!        On entry, the matrix A.
!        On exit, A contains the triangular matrix R, or part of R, as
!        follows:
!        If m-k-l >= 0, then R is stored in A(1:k+l,n-k-l+1:n).
!        If m-k-l < 0, then the matrix
!                  ( R11     R12    R13 )
! 		 (  0      R22    R23 )
!        is stored in A(1:m,n-k-l+1:n).
! B      (input/output) REAL or COMPLEX array, shape (:,:) with
!        size(B,1) = p and size(B,2) = n.
!        On entry, the matrix B.
!        On exit, if m-k-l < 0, then R33 is stored in
!        B(m-k+1:l,n+m-k-l+1:n).
! ALPHA  (output) REAL array, shape (:) with size(ALPHA) = n
!        The real scalars alpha(i) , i = 1, 2,..., k+l.
! BETA   (output) REAL array, shape (:) with size(BETA) = n.
!        The real scalars beta(i) , i = 1, 2, ..., k+l.
!        Note: The generalized singular values of the pair (A,B) are
!        sigma(i) = ALPHA(i)/BETA(i), i = 1, 2, ...,  k+l.
!        If k + l < n, then ALPHA(k+l+1:n) = BETA(k+l+1:n) = 0.
! K, L   Optional (output) INTEGER.
!        The dimension parameters k and l.
! U      Optional (output) REAL or COMPLEX square array, shape (:,:) with
!        size(U,1) = m.
!        The matrix U .
! V      Optional (output) REAL or COMPLEX square array, shape (:,:) with
!        size(V,1) = p.
!        The matrix V .
! Q      Optional (output) REAL or COMPLEX square array, shape (:,:) with
!        size(Q,1) = n.
!        The matrix Q.
! IWORK  Optional (output) INTEGER array, shape(:) with size(IWORK) = n.
!        IWORK contains sorting information. More precisely, the loop
!             for i = k + 1, min(m, k + l)
!                   swap ALPHA(i) and ALPHA(IWORK(i))
!             end
!        will sort ALPHA so that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(n).
! INFO   Optional (output) INTEGER.
!        = 0: successful exit.
!        < 0: if INFO = -i, the i-th argument had an illegal value.
!        > 0: if INFO = 1, the algorithm failed to converge.
!        If INFO is not present and an error occurs, then the program is
!        terminated with an error message.
!----------------------------------------------------------------------
!  .. LOCAL PARAMETERS ..
   CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_GGSVD'
!  .. LOCAL SCALARS ..
   CHARACTER(LEN=1) :: LJOBU, LJOBV, LJOBQ
   INTEGER :: M, N, P, LINFO, ISTAT, ISTAT1, S1U, S2U, S1V, S2V, &
              S1Q, S2Q, LK, LL
!  .. LOCAL ARRAYS ..
   REAL(WP), TARGET :: LLU(1,1), LLV(1,1), LLQ(1,1)
   INTEGER, POINTER :: LIWORK(:)
   REAL(WP), POINTER :: WORK(:)
!  .. INTRINSIC FUNCTIONS ..
   INTRINSIC MAX, PRESENT, SIZE
!  .. EXECUTABLE STATEMENTS ..
   LINFO = 0; ISTAT = 0; M = SIZE(A,1); N = SIZE(A,2); P = SIZE(B,1)
   IF( PRESENT(U) )THEN; S1U = SIZE(U,1); S2U = SIZE(U,2); LJOBU = 'U'
   ELSE; S1U = 1; S2U = 1; LJOBU = 'N'; END IF
   IF( PRESENT(V) )THEN; S1V = SIZE(V,1); S2V = SIZE(V,2); LJOBV = 'V'
   ELSE; S1V = 1; S2V = 1; LJOBV = 'N'; END IF
   IF( PRESENT(Q) )THEN; S1Q = SIZE(Q,1); S2Q = SIZE(Q,2); LJOBQ = 'Q'
   ELSE; S1Q = 1; S2Q = 1; LJOBQ = 'N'; END IF
!  .. TEST THE ARGUMENTS
   IF( M < 0 .OR. N < N )THEN; LINFO = -1
   ELSE IF( SIZE(B,2) /= N .OR. P < 0 ) THEN; LINFO = -2
   ELSE IF( SIZE( ALPHA ) /= N )THEN; LINFO = -3
   ELSE IF( SIZE( BETA ) /= N )THEN; LINFO = -4
   ELSE IF( PRESENT(U) .AND. ( S1U /= MAX(1,M) .OR.  S2U /= M ) )THEN; LINFO = -7
   ELSE IF( PRESENT(V) .AND. ( S1V /= MAX(1,P) .OR.  S2V /= P ) )THEN; LINFO = -8
   ELSE IF( PRESENT(Q) .AND. ( S1Q /= MAX(1,N) .OR.  S2Q /= N ) )THEN; LINFO = -9
   ELSE
    IF (PRESENT(IWORK)) THEN
      LIWORK => IWORK
    ELSE
      ALLOCATE( LIWORK(MAX(1,N)), STAT=ISTAT )
    ENDIF
      IF( ISTAT == 0 ) THEN
         ALLOCATE( WORK( MAX(1,MAX(3*N,M,P)+N) ), STAT=ISTAT)
      END IF
      IF( ISTAT == 0 ) THEN
         IF( PRESENT(U) )THEN
           IF( PRESENT(V) )THEN
             IF( PRESENT(Q) )THEN
                CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), &
                         B, MAX(1,P), ALPHA, BETA, U, MAX(1,S1U), &
                        V, MAX(1,S1V), Q, MAX(1,S1Q), WORK, LIWORK, LINFO )
	      ELSE
	        CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), &
                         B, MAX(1,P), ALPHA, BETA, U, MAX(1,S1U), &
                        V, MAX(1,S1V), LLQ, MAX(1,S1Q), WORK, LIWORK, LINFO )
	      ENDIF
	     ELSE
	      IF( PRESENT(Q) )THEN
	        CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), &
                         B, MAX(1,P), ALPHA, BETA, U, MAX(1,S1U), &
                        LLV, MAX(1,S1V), Q, MAX(1,S1Q), WORK, LIWORK, LINFO )
	      ELSE
	        CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), &
                         B, MAX(1,P), ALPHA, BETA, U, MAX(1,S1U), &
                        LLV, MAX(1,S1V), LLQ, MAX(1,S1Q), WORK, LIWORK, LINFO )
	      ENDIF
	    ENDIF
	   ELSE
	     IF( PRESENT(V) )THEN
	       IF( PRESENT(Q) )THEN
	         CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), &
                         B, MAX(1,P), ALPHA, BETA, LLU, MAX(1,S1U), &
                        V, MAX(1,S1V), Q, MAX(1,S1Q), WORK, LIWORK, LINFO )
	       ELSE
	         CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), &
                         B, MAX(1,P), ALPHA, BETA, LLU, MAX(1,S1U), &
                        V, MAX(1,S1V), LLQ, MAX(1,S1Q), WORK, LIWORK, LINFO )
		ENDIF
	      ELSE
	       IF( PRESENT(Q) )THEN
	         CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), &
                         B, MAX(1,P), ALPHA, BETA, LLU, MAX(1,S1U), &
                        LLV, MAX(1,S1V), Q, MAX(1,S1Q), WORK, LIWORK, LINFO )
		ELSE
		 CALL GGSVD_F77( LJOBU, LJOBV, LJOBQ, M, N, P, LK, LL, A, MAX(1,M), &
                         B, MAX(1,P), ALPHA, BETA, LLU, MAX(1,S1U), &
                        LLV, MAX(1,S1V), LLQ, MAX(1,S1Q), WORK, LIWORK, LINFO )
		ENDIF
	       ENDIF
	      ENDIF
         IF( PRESENT(K) ) K = LK
         IF( PRESENT(L) ) L = LL
      ELSE; LINFO = -100; ENDIF
      DEALLOCATE(WORK, STAT=ISTAT1)
      IF (.NOT. PRESENT(IWORK)) DEALLOCATE(LIWORK, STAT=ISTAT1)
   ENDIF
   CALL ERINFO(LINFO,SRNAME,INFO,ISTAT)
END SUBROUTINE DGGSVD_F95


text:
c###esmod.for
      SUBROUTINE ESMOD
C--------------------------------
C
C     THIS ROUTINE FINDS SPORADIC E MODE AND LOSS INFORMATION
C     (SAME AS SUBROUTINE REGMOD FOR THE OTHER LAYERS)
C
C AC(3) IS ABSORPTION LOSS FACTOR, SEE LUFFY.
C BC(3) IS DENOMINATOR OF LOSS FACTOR.
c
      COMMON/ANOIS/ATNU,ATNY,CC,TM,RCNSE,DU,DL,SIGM,SXGU,SXGL,KJ,JK
      COMMON / CONTRL / IELECT(3), KTOUT(12), MONTHS(12), SUNSP(12),
     A IANTOU, ICARD, INTEG, IRED, ISOUT, ISPROC, ISSN, ITYPE, JDASH,
     B JFREQ, JLONG, KCARD, KRUN, MAPIN, MAXNAM, MONOLD, MOREM, MORES,
     C NUMNAM, NUPROC, MAXMET, MSPEC, M100
      COMMON/SIGD/ DSL,ASM,DSU,AGLAT,DSLF,ASMF,DSUF,ACAV,FEAV,AFE,BFE
     A ,HNU,HTLOSS,XNUZ,XVE
      COMMON/MUFS/EMUF(24),F1MUF(24),F2MUF(24),ESMUF(24),ALLMUF(24),FOT
     A(24),XLUF(24),HPF(24),ANGMUF(24),MODMUF,SIGL(4),SIGU(4),DELMUF(4)
     B ,HPMUF(4),HTMUF(4),FVMUF(4),AFMUF(4),NHOPMF(4),YFOT(4),YHPF(4)
     C ,YMUF(4)
      COMMON /CON /D2R, DCL, GAMA, PI, PI2, PIO2, R2D, RZ, VOFL
      COMMON /DON /ALATD, AMIN, AMIND, BTR, BTRD, DLONG, DMP, ERTR, GCD,
     1 GCDKM, PMP, PWR, TLAT, TLATD, TLONG, TLONGD, RSN, SIGTR, RLAT,
     2 RLATD,RLONG,RLONGD,BRTD,FLUX,ULAT,ULATD,ULONG,ULONGD,SSN,D90R,
     3 D50R,D10R,D90S,D50S,D10S
      COMMON /RON /CLAT(5), CLONG(5), GLAT(5), RD(5), FI(3,5), YI(3,5),
     1HI(3,5), HPRIM(30,5), HTRUE(30,5), FVERT(30,5),KM,KFX, AFAC(30,5),
     2HTR(50,3), FNSQ(50,3)
      COMMON / ZON / ABPS(7), CREL(7), EFF(7), FLDST(7), GRLOS(7),
     1 HN(7), HP(7), PROB(7), RELY(7), RGAIN(7), SIGPOW(7), SN(7),
     2 SPRO(7), TGAIN(7), TIMED(7), TLOSS(7), B(7), FSLOS(7), ADV(7),
     3 OBF(7),NMODE(7),TLLOW(7),TLHGH(7)
      COMMON /TON /ADJ, ADS, GNOS, GOT, REL, SL, SLS
     1, SPR, SU, SUS, XNOISE, ZNOISE, NF
      COMMON/FRQ/FREA(13),FREL(29),FREQ,JMODE,ITXRCP(2)
      COMMON /ES /FS (3, 5), HS (5)
      COMMON /GEOG /GYZ (5), RAT (5), GMDIP (5), CLCK (5), ABIY (5), ART
     1IC (5), SIGPAT (5), EPSPAT (5)
      COMMON / RAYS / ANG(40), IFOB(40,30,5), NANG
C ES MODES
C     ONE SAMPLE MODEL
      K=1
C.....SELECT SAMPLE AREA - ALL MODES WILL BE AT LEAST THIS GOOD
      DO 95 IS = 1,KM
      IF( FS(2,K) - FS(2,IS) ) 95,95,90
   90 K= IS
   95 CONTINUE
C.....MAPS ARE ES OR E - DO NOT USE AT LOW END
      FSDEAD = IFOB(1,3,JMODE)
      FSDEAD = FSDEAD /1000.
      FSDEAD = AMIN1(FSDEAD,3.)
C.....ONLY 2 ES MODES - PRESET ARRAYS
      DO 97 IHT = 4,5
      OBF(IHT)  = 1000.
      ADV(IHT)  = 1000.
      FSLOS(IHT)  = 1000.
      TLOSS(IHT)  = 1000.
      ABPS (IHT)  = 1000.
      EFF  (IHT)  = 0.0
      GRLOS(IHT)  = 1000.
      RGAIN(IHT)  = 0.0
      TGAIN(IHT)  = 0.0
      HN   (IHT)  = -1.
      PROB (IHT)  = 0.001
      CREL (IHT)  = -1000.
      RELY (IHT)  = .001
      SPRO (IHT)  = 0.001
      FLDST(IHT)  = -1000.
      SIGPOW(IHT)  = -1000.
      SN(IHT)     = -1000.
      HP(IHT)     =  -1.
      B (IHT)     =  -1.
      NMODE(IHT)  =   5
      TLLOW(IHT)  =  10.
      TIMED(IHT) = -1.
      TLHGH(IHT)  = 10.
   97 CONTINUE
      IF(FREQ.le.FSDEAD) return
      IF(FS(2,K).le.0.) return
C.....VIRTUAL HEIGHTS SAME AS TRUE HEIGHTS
      SDMAX = 2.*RZ*(PIO2 - ASIN(1./(1.+ HS(K)/RZ) ) )
      IHSRT = GCDKM/SDMAX + 1.
      IHSTP = 2
C.....DO ONLY 2 Es HOPS
      IF(IHSRT .gt.2) return
      IH =3
C.....ABSORPTION LOSS - F2 MODE
      AC = 677.2 * ACAV/ ( (FREQ + GYZ(K)) ** 1.98  + 10.2 )
      DO 150 IHOP = IHSRT, IHSTP
      IH = IH + 1
      GP = IHOP
      GHOP = GCD / GP
      THET = 0.5 * GHOP
      TANS=SIN(THET)/(1.-COS(THET)+HS(K)/RZ)
      PSI=ATAN(TANS)
      SECS=1./COS(PSI)
      SFVMOD = FREQ/SECS
C.....MUF AT THIS HOP DISTANCE
      ESD=FS(2,K)*SECS
      DEL=PIO2-THET-PSI
      CDEL = COS (DEL)
      ADEL = DEL * R2D
      IF(ADEL.lt.AMIND)  go to 150
C.....GROUP PATH
      PATH = 2. * SIN (.5 * GHOP) * (RZ + HS (K)) / CDEL
      HOP = IHOP
      PATH = HOP * PATH
C.....FREE SPACE LOSS
      SFLOS = 32.45 + 20.*ALOG10(PATH * FREQ)
      SINP = RZ * CDEL / (RZ + 100.)
      COSP = SQRT(1. - SINP * SINP)
      SABPS = AC/COSP
      ADX =0.0
C.....CHECK IF Es LAYER < E LAYER (IF YES REMOVE E LAYER BENDING ABOVE HS)
      if(FI(1,K).gt.SFVMOD) ADX=AFE + BFE*ALOG(SFVMOD/FI(1,K) )
      SABPS = SABPS + ADX
C.....PROBABILITY OF REFLECTION
      DUMMY = YMUF(4)
      PROS = PRBMUF(FREQ,ESD,DUMMY,4)
C.....LIMIT LOSS
      PROS = AMIN1(PROS,.90)
      REF = 8.9136 * PROS ** (-0.7)
C.....LOWER DECILE FOR LOSS USE FOT
      ESD  = FS(1,K) * SECS
      DUMMY = YFOT(4)
      PS = PRBMUF(FREQ,ESD,DUMMY,4)
      PS = AMIN1(PS,.9)
      REFL = 8.9136 * PS ** (-0.7)
C.....UPPER DECILE -   HPF
      ESD =  FS(3,K)* SECS
      DUMMY = YHPF(4)
      PS = PRBMUF(FREQ,ESD,DUMMY,4)
      PS  = AMIN1(PS,.9)
      REFU = 8.9136 * PS ** (-0.7)
C.....LOWER DECILE WITH AURORAL
      TLLOW(IH) = DSL + HOP * (REFL - REF)
C.....UPPER DECILE WITH AURORAL
      TLHGH(IH) = DSU + HOP * (REF - REFU)
C.....MAXIMUM OBSERVED IS 25.
      TLLOW(IH) = AMIN1(TLLOW(IH),25.)
      TLHGH(IH) = AMIN1(TLHGH(IH),25.)
      Y = 0.0
C.....GROUND LOSS
      DO 130 IG = 1,KM
      IGX = -IG
      CALL GAIN(IGX,DEL,FREQ,YG,GEFF)
  130 Y = Y + YG
      XKM = KM
      SGRLOS = Y / XKM
C.....ANTENNAS
C125PC      if(MSPEC.ne.125)then
        CALL GAIN(1,DEL,FREQ,STGAIN,STEFF)
        CALL GAIN(2,DEL,FREQ,SRGAIN,DUMMY)
        EFF(IH) = DUMMY
C125PC      else
c.......set gains and EFF to 0 dB or unity
C125PC        stgain=0.
C125PC        srgain=0.
C125PC        EFF(ih)=0.
C125PC      ENDIF
      xtlos = SFLOS + HOP*(SABPS + REF + ADX ) + (HOP -1.) * SGRLOS
     A    - SRGAIN  - STGAIN  + ASM
C.....FIELD STRENGTH
      FLDST(IH) = 107.2 + PWRDB(FREQ) + 20.*ALOG10(FREQ)-XTLOS-SRGAIN
C.....MEDIAN SIGNAL POWER
      SIGPOW(IH) = PWRDB(FREQ) - XTLOS
      DUMMY = YMUF(4)
      PROS = PRBMUF(FREQ,DUMMY,DUMMY,4)
      IHT = IH
C.....OBSCURATION FOR F LAYER NOT USED
      OBF(IHT)    =  + 8.9136*PROS ** (-0.7)
C.....DEVIATION LOSS
      ADV(IHT)    = 0.0
C.....FREE SPACE
      FSLOS(IHT)  =  SFLOS
C.....MEDIAN TRANSMISSION LOSS
      TLOSS(IHT)  =  XTLOS
C.....ABSORPTION LOSS
      ABPS (IHT)  =  AC/COSP + ADX
C.....GROUND LOSS
      GRLOS(IHT)  =  SGRLOS
C.....ANTENNA GAINS
      RGAIN(IHT)  = SRGAIN
      TGAIN(IHT)  = STGAIN
      HN(  IHT)  = HOP
C.....SIGNAL POWER
      SIGPOW(IHT) = PWRDB(FREQ) - XTLOS
C.....SIGNAL-TO-NOISE RATIO
C------------    CHANGED 9/24/91  (LONG PATH RCVR EFF CORRECTION)  FJR
C------------            SN    (IHT) = SIGPOW(IHT) - RCNSE
      SN    (IHT) = SIGPOW(IHT) - RCNSE - EFF(IHT)
C.....HEIGHT
      HP    (IHT) = HS(K)
C.....RADIATION ANGLE
      B     (IHT) =  ADEL
C.....MODE
      NMODE (IHT) =  4
C.....F. DAYS
      PROB  (IHT) =  PROS ** IHOP
C.....TIME DELAY
      TIMED (IHT) =  PATH/VOFL
  150 CONTINUE
      RETURN
      END
C--------------------------------

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


text:
      ! -------------------------------------------------------------
      !
      ! Copyright (C) 2021 The Simons Foundation
      !
      ! Author: Jason Kaye
      !
      ! -------------------------------------------------------------
      !
      ! libdlr is licensed under the Apache License, Version 2.0 (the
      ! "License"); you may not use this file except in compliance with
      ! the License.  You may obtain a copy of the License at
      !
      !     http://www.apache.org/licenses/LICENSE-2.0
      !
      ! Unless required by applicable law or agreed to in writing,
      ! software distributed under the License is distributed on an "AS
      ! IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
      ! express or implied.  See the License for the specific language
      ! governing permissions and limitations under the License.
      !
      ! -------------------------------------------------------------

      program sc_it

      ! Demonstration of discrete Lehmann representation for imaginary
      ! time Green's function with semi-circular density of states,
      ! rho(omega) = sqrt(1-omega^2), using imaginary time sampling.
      !
      ! The DLR expansion is formed from samples at the DLR imaginary
      ! time nodes, and then evaluated both in imaginary time and
      ! Matsubara frequency domains, where its accuracy is measured.
      !
      ! The Green's function is sampled by evaluating the Green's
      ! function using the Lehmann representation with the known
      ! spectral density. The integral in the Lehmann representation is
      ! computed by high-order numerical integration. This reference
      ! Green's function is also used to test the accuracy of the DLR
      ! expansion.

      implicit none
      integer ntst_it,ntst_mf
      real *8 lambda,eps,beta

      ! Input parameters

      lambda = 1000.0d0 ! DLR high energy cutoff
      eps = 1.0d-14     ! DLR error tolerance

      beta = 1000.0d0   ! Inverse temperature

      ntst_it = 2000    ! # imaginary time test points
      ntst_mf = 2000    ! Matsubara frequency test point cutoff

      ! Main subroutine

      call sc_it_main(lambda,eps,ntst_it,ntst_mf,beta)

      end program sc_it

      subroutine sc_it_main(lambda,eps,ntst_it,ntst_mf,beta)

      implicit none
      integer ntst_it,ntst_mf
      real *8 lambda,eps,beta

      integer npg,npo,pg,i,j,r
      integer, allocatable :: it2cfp(:),mf_tst(:)
      real *8 one
      real *8, allocatable :: it2cf(:,:),dlrit(:),dlrrf(:),g(:),gc(:)
      real *8, allocatable :: xgl(:),wgl(:),xgj(:),wgj(:),pbpg(:)
      real *8, allocatable :: it_tst(:),gtst_it(:),gtrue_it(:)
      complex *16, allocatable :: gtst_mf(:),gtrue_mf(:)

      one = 1.0d0

      ! Get DLR frequencies, imaginary time grid

      r = 500 ! Upper bound on DLR rank

      allocate(dlrrf(r),dlrit(r))

      call dlr_it_build(lambda,eps,r,dlrrf,dlrit)

      ! Get imaginary time values -> DLR coefficients matrix (LU form)

      allocate(it2cf(r,r),it2cfp(r))

      call dlr_it2cf_init(r,dlrrf,dlrit,it2cf,it2cfp)

      ! Initialize reference Green's function evaluator (note: this is
      ! not a DLR-related operation; rather, we are initializing a
      ! special-purpose subroutine to evaluate the Green's function with
      ! a semi-circular spectral density. The initialization and
      ! evaluation subroutines are defined at the bottom of this file.)

      pg = 24
      npg = max(ceiling(log(lambda)/log(2.0d0)),1)

      allocate(xgl(pg),wgl(pg),xgj(pg),wgj(pg),pbpg(2*npg+1))

      call gfun_init(pg,npg,pbpg,xgl,wgl,xgj,wgj)

      ! Sample G at imaginary time nodes

      allocate(g(r))

      do i=1,r

        call gfun_it(pg,npg,pbpg,xgl,wgl,xgj,wgj,beta,dlrit(i),&
          g(i))

      enddo

      ! Get DLR coefficients of G

      allocate(gc(r))

      call dlr_it2cf(r,1,it2cf,it2cfp,g,gc)

      ! Get imaginary time test points in relative format

      allocate(it_tst(ntst_it))

      call eqpts_rel(ntst_it,it_tst)

      ! Evaluate DLR on imaginary time test grid

      allocate(gtst_it(ntst_it))

      do i=1,ntst_it

        call dlr_it_eval(r,1,dlrrf,gc,it_tst(i),gtst_it(i))

      enddo

      ! Get Matsubara frequency test points

      allocate(mf_tst(2*ntst_mf+1))

      do i=1,2*ntst_mf+1
        mf_tst(i) = -ntst_mf+i-1
      enddo

      ! Evaluate DLR on Matsubara frequency test grid

      allocate(gtst_mf(2*ntst_mf+1))

      do i=1,2*ntst_mf+1

        call dlr_mf_eval(r,1,dlrrf,-1,gc,mf_tst(i),gtst_mf(i))

      enddo

      ! Evaluate true Green's function on imaginary time test grid

      allocate(gtrue_it(ntst_it))

      do i=1,ntst_it

        call gfun_it(pg,npg,pbpg,xgl,wgl,xgj,wgj,beta,it_tst(i),&
          gtrue_it(i))

      enddo

      ! Evaluate true Green's function on Matsubara frequency test grid

      allocate(gtrue_mf(2*ntst_mf+1))

      do i=1,2*ntst_mf+1

        call gfun_mf(pg,npg,pbpg,xgl,wgl,xgj,wgj,beta,mf_tst(i),&
          gtrue_mf(i))

      enddo

      ! Output errors

      write(6,*) ''
      write(6,*) 'lambda = ',lambda
      write(6,*) 'epsilon = ',eps
      write(6,*) 'beta = ',beta
      write(6,*) ''
      write(6,*) 'DLR rank = ',r
      write(6,*) ''
      write(6,*) 'Imag time max error = ',maxval(abs(gtst_it-gtrue_it))
      write(6,*) 'Mats freq max error = ',maxval(abs(gtst_mf-gtrue_mf))
      write(6,*) ''

      end subroutine sc_it_main

      subroutine gfun_init(n,np,pbp,xgl,wgl,xgj,wgj)

      ! Initialize subroutines gfun_it and gfun_mf, which evaluate
      ! Green's function with semi-circular density at a imaginary time
      ! and Matsubara frequency points, respectively

      implicit none
      integer n,np
      real *8 pbp(2*np+1),xgl(n),wgl(n),xgj(n),wgj(n)

      integer i
      real *8 one

      one = 1.0d0

      ! Gauss-Legendre and Gauss-Jacobi quadrature nodes and weights

      call cdgqf(n,1,0.0d0,0.0d0,xgl,wgl)
      call cdgqf(n,4,0.5d0,0.0d0,xgj,wgj)

      ! Panel endpoints for composite quadrature rule

      pbp(np+1) = 0*one
      do i=1,np
        pbp(np+i+1) = one/2**(np-i)
      enddo
      pbp(1:np) = -pbp(2*np+1:np+2:-1)

      end subroutine gfun_init

      subroutine gfun_it(n,np,pbp,xgl,wgl,xgj,wgj,beta,t,val)

      ! Evaluate Green's function with semi-circular density at an
      ! imaginary time point using the Lehmann representation. Integral
      ! is computed by high-order composite Gauss-Legendre quadrature,
      ! with Gauss-Jacobi quadrature at endpoint panels to treat square
      ! root singularities.

      implicit none
      integer n,np
      real *8 pbp(2*np+1),xgl(n),wgl(n),xgj(n),wgj(n),beta,t,val

      integer ii,jj
      real *8 one,a,b,x,tt
      real *8, external :: kfunf

      one = 1.0d0

      ! Treat 0.5<tau<1, stored in relative format, by symmetry

      tt = abs(t)

      ! Composite Gauss quadrature

      val = 0.0d0
      do ii=2,2*np-1
        a = pbp(ii)
        b = pbp(ii+1)
        do jj=1,n
          x = a+(b-a)*(xgl(jj)+one)/2
          val = val + (b-a)/2*wgl(jj)*kfunf(tt,beta*x)*&
            sqrt(one-x**2)
        enddo
      enddo

      a = one/2
      b = one
      do jj=1,n
        x = a+(b-a)*(xgj(jj)+one)/2
        val = val + ((b-a)/2)**(1.5d0)*wgj(jj)*&
          kfunf(tt,beta*x)*sqrt(one+x)
      enddo

      a = -one
      b = -one/2
      do jj=1,n
        x = a+(b-a)*(-xgj(n-jj+1)+one)/2
        val = val + ((b-a)/2)**(1.5d0)*wgj(n-jj+1)*&
          kfunf(tt,beta*x)*sqrt(one-x)
      enddo

      end subroutine gfun_it

      subroutine gfun_mf(n,np,pbp,xgl,wgl,xgj,wgj,beta,m,val)

      ! Evaluate Green's function with semi-circular density at a
      ! Matsubara frequency point using the Lehmann representation. Integral
      ! is computed by high-order composite Gauss-Legendre quadrature,
      ! with Gauss-Jacobi quadrature at endpoint panels to treat square
      ! root singularities.

      implicit none
      integer n,np,m
      real *8 pbp(2*np+1),xgl(n),wgl(n),xgj(n),wgj(n),beta
      complex *16 val

      integer ii,jj
      real *8 one,a,b,x
      complex *16, external :: kfunmf

      one = 1.0d0

      val = 0.0d0
      do ii=2,2*np-1
        a = pbp(ii)
        b = pbp(ii+1)
        do jj=1,n
          x = a+(b-a)*(xgl(jj)+one)/2
          val = val + (b-a)/2*wgl(jj)*kfunmf(2*m+1,beta*x)*&
            sqrt(one-x**2)
        enddo
      enddo

      a = one/2
      b = one
      do jj=1,n
        x = a+(b-a)*(xgj(jj)+one)/2
        val = val + ((b-a)/2)**(1.5d0)*wgj(jj)*&
          kfunmf(2*m+1,beta*x)*sqrt(one+x)
      enddo

      a = -one
      b = -one/2
      do jj=1,n
        x = a+(b-a)*(-xgj(n-jj+1)+one)/2
        val = val + ((b-a)/2)**(1.5d0)*wgj(n-jj+1)*&
          kfunmf(2*m+1,beta*x)*sqrt(one-x)
      enddo

      end subroutine gfun_mf

max_stars_repo_name: jasonkaye/libdlr

max_issues_repo_name: jasonkaye/libdlr

max_forks_repo_name: jasonkaye/libdlr


text:
! { dg-do compile }
!
! PR fortran/31559
! Do not allow assigning to external functions
!
! Contributed by Steve Kargl <sgk@troutmask.apl.washington.edu>
!
module mod
  implicit none
contains
  integer function bar()
    bar = 4
  end function bar

  subroutine a()
   implicit none
   real :: fun
   external fun
   interface
     function funget(a)
       integer :: a
     end function
     subroutine sub()
     end subroutine sub
   end interface
   sub = 'a'  ! { dg-error "Expected VARIABLE" }
   fun = 4.4  ! { dg-error "Expected VARIABLE" }
   funget = 4 ! { dg-error "is not a VALUE" }
   bar = 5    ! { dg-error "is not a VALUE" }
  end subroutine a
end module mod

end

max_stars_repo_name: vidkidz/crossbridge

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


module FOO
end module FOO
module bar
use foo
end module bar


text: The http://telperion.info Order of Telperion is a mystical knighthood, based in Davis and not at all connected to the local fraternity called Alpha Gamma Omega. All that is done by the Order is secret, except for that which is not secret; regardless, all is done in the appropriate manner, which, incidentally, is secret. Sorry. None of Telperions members have ever been a part of Campus Crusade for Chaos and Confusion and certainly do not recommend that you Steal Bainer Hall this Building.


text:
University of California Agriculture and Natural Resources (ANR) is a division of the UC Office of the President headquartered in Oakland. While ANR is not specifically a UC Davis department, it is affiliated with the campus, specifically the College of Agricultural and Environmental Sciences and the School of Veterinary Medicine. ANR is also affiliated with UC Berkeley and UC Riverside. ANR operates two chief programs: Cooperative Extension and the Agricultural Experiment Station, which are mandated by the fact that the University of California is a land grant university.

Image(Dooley.jpg, right, ANR VP, Daniel Dooley) ANR is led by the Vice President, Daniel Dooley, who provides administrative leadership for the University of Californias programs in Agriculture and Natural Resources.

Prior to 1998, ANR was known as the Division of Agriculture and Natural Resources (DANR), hence the commonly used names of the DANR Building and DANR Lab (now officially the ANR Building and ANR Lab) on campus.

UC ANR has units in several buildings on campus, including
  ANR Building on 225 Hopkins Road
  ANR Lab in Hoagland Hall
  http://anrcs.ucdavis.edu/ Communication Services on 1850 Research Park Drive
  4H Center for Youth Development in Hart Hall
  http://aic.ucdavis.edu/ Agricultural Issues Center in Hunt Hall
  http://efnep.ucdavis.edu/ Expanded Food and Nutrition Education Program (EFNEP) in Meyer Hall
  http://www.ipm.ucdavis.edu/ Statewide Integrated Pest Management Program  in Robbins Hall Annex

They also have a http://ucanr.org website.


text:
Job Board Job hunting can be a royal pain in the youknowwhat.  This page may help streamline that process.  Note that many places that have online applications also have paper applications available instore.
 Retail
 Bookstores
Bring in resume:
  Avid Reader
  OffCampus Books
  Newsbeat

Other bookstores (please indicate the type of application process for any of these, if you know):

  Bizarro World
  Davis Christian Bookroom

 Bicycle Shops and Repair
Apply Online:
  Kens Bike & Ski (http://kensbikeski.com/page.cfm?pageID222 Apply Here)
Print application, then submit instore:
  B&L Bike Shop (http://blbikeshop.com/merchant/504/files/EmploymentApplication.pdf Print Here)
Apply in person:
  APEX Cycles  Come talk to Aaron and see if you are a fit for the shop. 12 week probationary period required

Other Bicycle Shops (please indicate the type of application process for any of these, if you know):

  Davis Bike Exchange
  Freewheeler
  Wheelworks

 Mobile Phones Cell Phone Retailers
Other cell phone retailers (please indicate the type of application process for any of these, if you know):

  Absolute Cellular
  AT&T Wireless
  Heron Technologiesis now Davis Computer
  Mobile Connections
  Parrot Cellular

 Children ChildOriented
Apply online:
  Gap Gap Kids (http://www.gapinc.com/public/Careers/car_jobsearch.shtml Apply Here)

Other childoriented stores (please indicate the type of application process for any of these, if you know):

  Mother and Baby Source
  Teach Your Children

 Clothing Stores
Apply Online:
  The Gap (http://www.gapinc.com/public/Careers/car_jobsearch.shtml Apply Here)
Pick up application instore:
  The Wardrobe
Other clothing stores (please indicate the type of application process for any of these, if you know):

  Ground Zero
  Nina and Tom
  Preeti Girl
  Pinkadot
  Renew Denim
  Riki

 Other
Apply Online:
  Petco (http://www.petco.com/petco_Page_PC_retailstoreopps_Nav_117.aspx Apply Here) Online only
  Cost Plus
  Office Max
Print application, then submit instore:
  Dollar Tree (http://www.dollartree.com/downloads/career/Sales_Associate_Application.pdf Print Here)

Pick up an application instore:
  Ace Hardware
  Davis Home Trends
  The Paint Chip

Submit resume:

  Ink Monkey Graphics

Other places (please indicate the type of application process for any of these, if you know):

  Bath and Body Works
  Big 5 Sporting Goods
  H2O to go
  Davis Academy Driving School
  Davis Shoe Shop & Pedorthic
  Fleet Feet
  Hibbert Lumber
  Nestware
  Rivers to Reef
  Tibet Nepal
  Valley Wine Company
  Volleys Tennis Shop

 Food and Restaurants
 Casual Dining
Apply online:
  Noahs Bagels (http://www.noahs.com/#/careers/ Apply Here)
Print application, then submit instore:
  Beach Hut Deli (http://www.beachhutdeli.com/ Print here)
  Dos Coyotes (http://www.doscoyotes.net/dosapp.pdf Print here)
Pick up application instore:
  Village Bakery
Submit resume:
  Delta of Venus
Other casual dining (please indicate the type of application process for any of these, if you know):

  Ali Baba
  Chuys Taqueria
  El Mariachi
  El Toro Bravo
  Golden Sun
  Hot Dogger
  Konditorei
  Mr. Pickles Sandwich Shop
  Taqueria Davis
  Taqueria El Burrito
  Taqueria Guadalajara
  The Posh Bagel
  The Graduate

 Coffee Shops
Apply online:

  Starbucks (http://www.starbucks.com/aboutus/job_search_pop.asp Apply for a Corporate Position Here)
Print application, then submit instore:
  Common Grounds (http://www.commongroundsdavis.com/files/Application.pdf Print Here)
  Peets Coffee (http://www.peets.com/careers/job_postings.asp?cm_refooter_career_text&cm_spcareer_footer_text Print Here)
  Starbucks (http://media.starbucks.com.edgesuite.net/dotcom/media/pdf/Starbucks_Application.pdf Print Retail Application Here)
Pick up application instore:
  Cloud Forest Cafe
  Mocha Joes
Bring in resume:
  Mishkas

Other coffee shops (please indicate the type of application process for any of these, if you know):

  3rd & U Cafe
  Barista Brew

 Dinein
Submit resume online
  Cafe Bernardo (http://www.paragarys.com/go/prg/careers/submityourresume/index.cfm Submit here, but you can also apply inperson)

Print application, then submit instore:
  Bistro 33 (http://www.bistro33.com/bistro33_davis/SRO_Job_App.pdf Print here)
  Caffe Italia (http://www.dancingtomato.com/employment.html Print here)
  Plutos (http://www.plutosfreshfood.com/documents/PlutosEmpApp_10.19.06.pdf Print here)
  Sushi Unlimited (http://www.sushiunlimited.net/empApp.pdf Print Here)

Pick up application instore:
  Applebees
Submit resume instore:
  Osteria Fasulo
  Tucos Wine Market and Cafe

Other dinein restaurants (please indicate the type of application process for any of these, if you know):

  Black Bear Diner
  Burgers and Brew
  Café Méditerranée
  Cindys
  Crepeville
  Davis Noodle City
  Davis Oshio Cafe
  Davis Sushi Buffet Japanese Restaurant
  Ding How
  Farmers Kitchen Cafe
  Fuji Chef
  Froggys
  G Street Pub now G St. Wunderbar
  Great Wall of China II
  Hoa Viet
  Hometown Chinese Food
  House of Chang
  Hunan
  IHOP
  Jade Garden
  Jusco
  Kathmandu Kitchen
  KetMoRee
  Little Prague
  Moshi Moshi
  Nobu Hiro
  Old Teahouse
  Plainfield Station
  Rajas
  Red Orchid Restaurant
  Seasons
  Sams Mediterranean
  Shanghai Town
  Silver Dragon
  Sophias Thai Kitchen
  Sudwerk
  Sunrise Restaurant
  Symposium
  Thai Nakorn
  Thai Recipes
  Mustard Seed The Mustard Seed
  Wildhorse Grill
  Wok of Flame
  Zen Toro

 Fast Food
Apply online:
  Burger King (http://www.bk.com/companyinfo/careers/rmp.aspx Apply Here)
  Chipotle (http://www4.recruitingcenter.net/clients/chipotle/publicjobs/controller.cfm?jbactionApplyToJob Apply Here)
  Del Taco (http://www.deltaco.com/page42.html Apply Here)
  McDonalds (http://www.mcstate.com/careers/ Apply Here)
  Jamba Juice (http://www.jambajuice.com/jobs_careers#/job_listings/ Apply Here)
  Panda Express (http://www.pandacareers.com/index.asp Apply Here)
  Taco Bell (http://www.yumcareers.com/jobs_framesetup.html Apply Here)
Print application, then return instore:
  InnOut (http://www.innout.com/inoapp.pdf Print Here)
  Ohana Hawaiian BBQ (http://www.ohanabbqcorp.com/pdf/ohib_jobapp.pdf Print Here)
Pick up application instore:
  Carls Jr.
  Chipotle
  Habit Burger
Other fast food restaurants (please indicate the type of application process for any of these, if you know):

  Dairy Queen
  Jack in the Box
  KFC
  Quickly (http://www.quicklyusa.com/jobopportunity.html Print Here)
  Redrum Burger

 Pizza
Apply online:
  Papa Murphys (http://www.selfmgmt.com/clients/papamurphys/ Apply Here)
  Pizza Guys (https://www.pizzaguys.com/home/jobs.cfm Apply Here)
  Round Table Pizza (http://www.roundtablepizza.com/RTP/HI/ Apply Here)
Print application, then return instore:
  Cenarios Pizza (http://www.cenariospizza.com/ Print Here, click downloads and scroll down)
  Woodstocks Pizza (http://woodstocksdavis.com/pg/Jobs.htm Print Here)
Pick up an application instore:
  Lamppost Pizza
Email resume:
  Uncle Vitos Slice of N.Y. (email to jobs@paesanos.biz)
Other pizza places (please indicate the type of application process for any of these, if you know):

  Dominos Pizza
  Little Caesers
  Original Steves

 Sandwiches
Apply online:
  Subway (http://www.subway.com/subwayroot/AboutSubway/Employment/Local/index.aspx Apply Here)
Print application, then return instore:
  Beach Hut Deli (http://www.beachhutdeli.com/BHD_emp_ap.pdf Print Here)
Pick up an application instore:
Other sandwich places (please indicate the type of application process for any of these, if you know):

  Quiznos
  Togos
  Zias Delicatessen

 Ice Cream and Dessert
Print application, then return instore:
  Cultive Frozen Yogurt (http://cultiveyogurt.com/images/cultive_application.pdf Print Here)
Pick up an application instore:
  Baskin Robbins
Email Resume:
  Yolo Berry Yogurt (email to yoloberry@gmail.com)

Other ice cream and dessert places (please indicate the type of application process for any of these, if you know):

  Ciocolat
  Fluffy Donuts
  IceKrimski cafe
  Let them eat Cake
  SugarPlum
  The Candy House of Davis
  Tutti Frutti Frozen Yogurt

 Grocery, Convenience and Drug Stores
Apply Online:
  7Eleven (http://www.7eleven.com/Careers/SalesAssociate/tabid/143/Default.aspx Apply Here)
  Nugget Market (https://jobsearch.unicru.com/JLohome.aspx?rscid{8eb07c836aff42b0e543ed3a6af31c61}&RPCR Apply Here) Online only
  Safeway (http://www.safeway.com/Employment/#iframetop Apply Here) Online applications only
  Savemart (http://www.savemart.com/storelevel.php Apply Here)
Print application, then submit instore:
  Davis Food Coop (http://www.davisfood.coop/employment.html Print Here)
Pick up application instore:
  Circle K
  Rite Aid
Other Grocery, Convenience and Drug Stores (please indicate the type of application process for any of these, if you know):

  AM/PM
  CVS Pharmacy
  El Macero Pharmacy
  Fast and Easy Mart
  Kims Mart
  Quick Shop Market
  The Olive Drive Market

 Services
 Copies
Apply Online:
  FedEx Office (formerly known as Kinkos) (http://fedex.hodesiq.com/careers/job_search.aspx Apply Here)
Pick up an application instore:
  Davis Copy Shop
Other copy places (please indicate the type of application process for any of these, if you know):

  Copyland

 Gyms and Fitness
Other gyms (please indicate the type of application process for any of these, if you know):

  Davis Athletic Club
  Davis Swim and Fitness Club
  Peak Performance
  Physical Edge

Shipping
Pick up application instore:
  Parcel Dispatch PDQ
Submit Resume:
  PostMarks
Other shipping places (please indicate the type of application process for any of these, if you know): UPS Store

 Tutoring
Apply Online:
  Kaplan (http://www.kaptest.com/Teach_for_Kaplan/Overview/CS_teach_about.html Apply Here)
  The Princeton Review (http://www.princetonreview.com/teachfortheprincetonreview.aspx?uidbadge Apply Here)
  Club Z! InHome Tutoring (http://www.clubztutoring.com/tutorform.php Apply or inquire here)
 Print application, then submit instore:
  OneonOne Tutoring and Educational Mapping (http://www.oneononedavis.com/index5.html Print Here)
Submit resume:
  Partners in Learning
Other tutoring venues (please indicate the type of application process for any of these, if you know):

  Kumon Math & Reading Centers
  LINK Tutoring Services

 Research

  Evaluating Grain Products Study Another way to earn money (and score free tasty food) is by participating in studies like this one, where food is given out free, and participants are paid for their time!
  Dietary Fiber Study Free food for 3 weeks and participants are paid handsomely for their time!

 Entertainment
 Theaters

  Regal Cinemas Davis Holiday 6
  Regal Cinemas Davis Stadium 5
  Varsity Theater

 Music

Submit resume:
  Armadillo Music
Other music places (please indicate the type of application process for any of these, if you know): Watermelon Music

 Video Rental
Other video rental places (please indicate the type of application process for any of these, if you know):
  Stonegate Video


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2019 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_bfgs_b2

CONTAINS

! The subroutine bfgs_b2 is an implementation of a BFGS algorithm with
! box constraints. It is motivated by the paper:
! Byrd et al., "A Limited Memory Algorithm for Bound Constrained Optimization",
! SIAM Journal on Scientific Computing 16, 1190 (1995)
! However, it is substantially different from that paper.
!
!                            GM'2019

! Note: lastGradient and lastParameters have to be 0.0 if this is the first call to the subroutine and these
!       vectors don't exist yet. Also historyLength has to be 0 in that case.

SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints,enabledConstraints,parameters,&
                   lastParameters,equalityLinCombi,equalityCriterion,maxHistoryLength,paramCorrections,&
                   gradientCorrections,iStep,mixParam,l_converged,convCrit)

   USE m_constants

   IMPLICIT NONE

   INTEGER,        INTENT(IN)    :: vecLength
   INTEGER,        INTENT(IN)    :: maxHistoryLength
   INTEGER,        INTENT(INOUT) :: iStep ! The index of the bfgs_b iteration
   REAL,           INTENT(IN)    :: mixParam
   REAL,           INTENT(IN)    :: equalityCriterion
   REAL,           INTENT(INOUT)    :: gradient(vecLength)
   REAL,           INTENT(INOUT) :: lastGradient(vecLength)
   REAL,           INTENT(IN)    :: minConstraints(vecLength)
   REAL,           INTENT(IN)    :: maxConstraints(vecLength)
   LOGICAL,        INTENT(IN)    :: enabledConstraints(vecLength)
   REAL,           INTENT(INOUT) :: parameters(vecLength)
   REAL,           INTENT(INOUT) :: lastParameters(vecLength)
   REAL,           INTENT(IN)    :: equalityLinCombi(vecLength)
   REAL,           INTENT(INOUT) :: paramCorrections(vecLength,maxHistoryLength)
   REAL,           INTENT(INOUT) :: gradientCorrections(vecLength,maxHistoryLength)
   LOGICAL,        INTENT(OUT)   :: l_converged
   REAL,           INTENT(IN)    :: convCrit

   INTEGER                :: i, j, stepIndex, numFreeParams, numFixedParams
   INTEGER                :: historyLength
   REAL, PARAMETER        :: eps = 1.0e-10
   REAL, PARAMETER        :: largeEps = 1.0e-6
   REAL                   :: temp, norm, subspaceEqualityCrit, maxDisp, dampingFactor
   REAL                   :: ddot

   LOGICAL                :: converged

   REAL                   :: lastGradCorrection(vecLength)
   REAL                   :: lastParamCorrection(vecLength)
   REAL                   :: prelimNextPoint(vecLength)
   REAL                   :: prelimGradient(vecLength)
   REAL                   :: tempVecA(vecLength)
   REAL                   :: subspaceGrad(vecLength)
   REAL                   :: subspaceNextPoint(vecLength)
   REAL                   :: subspaceLinCombi(vecLength)
   REAL                   :: diffPoint(vecLength)
   REAL                   :: diffGrad(vecLength)

   LOGICAL                :: freeParams(vecLength)

   REAL, ALLOCATABLE      :: hessMat(:,:), invHessMat(:,:)
   REAL, ALLOCATABLE      :: tempMatB(:,:), tempMatC(:,:)
   REAL, ALLOCATABLE      :: subspaceBasis(:,:)
   REAL, ALLOCATABLE      :: subspaceInvHessMat(:,:)

   WRITE(*,*) 'bfgs_b2 (input): parameters(:vecLength): ', parameters(:vecLength)
   WRITE(*,*) 'bfgs_b2 (input): gradient(:vecLength): ', gradient(:vecLength)

   WRITE(5000,*) '========================================================================='
   WRITE(5000,*) 'ENTERED BFGS_B2: '

   !!! Test start
!   DO i = 1, vecLength - 1
!      IF(ABS(gradient(i)).LT.eps) CYCLE
!      gradient(i) = sqrt(ABS(gradient(i)))*2.0*ATAN(gradient(i))/pi_const
!   END DO
   !!! Test end

   l_converged = .FALSE.

   ! 0. Add last data pair to history iff it is not obviously destructive for the positive definiteness of the Hessian matrix

   IF(ANY(lastGradient(:).NE.0.0).OR.ANY(lastParameters(:).NE.0.0)) THEN

      WRITE(2501,'(a,12f10.6)') '    gradient: ', gradient(:)
      WRITE(2501,'(a,12f10.6)') 'lastGradient: ', lastGradient(:)
      WRITE(2502,'(a,12f10.6)') '    parameters: ', parameters(:)
      WRITE(2502,'(a,12f10.6)') 'lastParameters: ', lastParameters(:)

      lastGradCorrection(:) = gradient(:) - lastGradient(:)
      lastParamCorrection(:) = parameters(:) - lastParameters(:)
      temp = ddot(vecLength,lastParamCorrection(:),1,lastGradCorrection(:),1)
      IF (temp.GT.eps) THEN
         ! Add data pair to history
         istep = istep + 1
         stepIndex = MOD(iStep - 1,maxHistoryLength) + 1
         paramCorrections(:,stepIndex) = lastParamCorrection(:)
         gradientCorrections(:,stepIndex) = lastGradCorrection(:)
      END IF

      ! Check for convergence and return if converged
      WRITE(*,*) 'MAXVAL(ABS(lastGradCorrection(:))): ', MAXVAL(ABS(lastGradCorrection(:)))
      WRITE(*,*) 'MAXVAL(ABS(lastParamCorrection(:))): ', MAXVAL(ABS(lastParamCorrection(:)))
      WRITE(*,*) 'MAXVAL(ABS(lastGradCorrection(:vecLength-1))): ', MAXVAL(ABS(lastGradCorrection(:vecLength-1)))
      WRITE(*,*) 'MAXVAL(ABS(lastParamCorrection(:vecLength-1))): ', MAXVAL(ABS(lastParamCorrection(:vecLength-1)))
!      IF(ALL(ABS(lastGradCorrection(:)).LT.convCrit).AND.ALL(ABS(lastParamCorrection(:)).LT.convCrit)) THEN
      IF(ALL(ABS(lastParamCorrection(:vecLength-1)).LT.convCrit)) THEN
         l_converged = .TRUE.
         RETURN
      END IF
   END IF

   ! 1. Construct Hessian matrix and its inverse.

   historyLength = MIN(istep,maxHistoryLength)

   ALLOCATE(hessMat(vecLength,vecLength))
   ALLOCATE(invHessMat(vecLength,vecLength))
   ALLOCATE(tempMatB(vecLength,vecLength))

   hessMat(:,:) = 0.0
   invHessMat(:,:) = 0.0
   tempMatB(:,:) = 0.0

   DO i = 1, vecLength
      hessMat(i,i) = 1.0
      invHessMat(i,i) = 1.0
   END DO

   IF(historyLength.EQ.0) THEN
      DO i = 1, vecLength
         hessMat(i,i) = 1.0 / mixParam
         invHessMat(i,i) = mixParam
      END DO
   END IF

   DO i = 1, historyLength
      stepIndex = MOD(iStep - historyLength + i - 1,maxHistoryLength) + 1

      ! First the Hessian matrix
      ! second correction term
      tempMatB(:,:) = 0.0
      tempVecA(:) = 0.0
      CALL dgemv('N',vecLength,vecLength,1.0,hessMat(:,:),vecLength,paramCorrections(:,stepIndex),1,0.0,tempVecA(:),1)
      norm = ddot(vecLength,paramCorrections(:,stepIndex),1,tempVecA(:),1)
      CALL dgemm('N','N',vecLength,vecLength,1,-1.0/norm,tempVecA(:),vecLength,tempVecA(:),1,0.0,tempMatB,vecLength)

      ! first correction term
      norm = ddot(vecLength,gradientCorrections(:,stepIndex),1,paramCorrections(:,stepIndex),1)
      CALL dgemm('N','N',vecLength,vecLength,1,1.0/norm,gradientCorrections(:,stepIndex),vecLength,gradientCorrections(:,stepIndex),1,1.0,tempMatB,vecLength)
      hessMat(:,:) = hessMat(:,:) + tempMatB(:,:)

      ! Now the inverse of the Hessian matrix
      tempMatB(:,:) = 0.0
      tempVecA(:) = 0.0
      CALL dgemv('N',vecLength,vecLength,1.0,invHessMat(:,:),vecLength,gradientCorrections(:,stepIndex),1,0.0,tempVecA(:),1)
      norm = ddot(vecLength,paramCorrections(:,stepIndex),1,gradientCorrections(:,stepIndex),1)

      ! second correction term
      CALL dgemm('N','N',vecLength,vecLength,1,1.0,tempVecA(:),vecLength,paramCorrections(:,stepIndex),1,0.0,tempMatB,vecLength)

      tempMatB(:,:) = tempMatB(:,:) + TRANSPOSE(tempMatB(:,:))
      tempMatB(:,:) = -tempMatB(:,:) / norm

      ! first correction term
      temp = norm + ddot(vecLength,gradientCorrections(:,stepIndex),1,tempVecA(:),1)
      CALL dgemm('N','N',vecLength,vecLength,1,temp/(norm*norm),paramCorrections(:,stepIndex),vecLength,paramCorrections(:,stepIndex),1,1.0,tempMatB,vecLength)
      invHessMat(:,:) = invHessMat(:,:) + tempMatB(:,:)

      !Check whether invHessMat is the inverse of hessMat (for debugging only)
      tempMatB(:,:) = 0.0
      CALL dgemm('N','N',vecLength,vecLength,vecLength,1.0,hessMat(:,:),vecLength,invHessMat(:,:),vecLength,0.0,tempMatB,vecLength)

   END DO

   ! 2. Determine next point based on the gradient, the inverse of the Hessian, and the box constraints

   freeParams(:) = .TRUE.
   prelimNextPoint(:) = parameters(:)
   prelimGradient(:) = gradient(:)

   WRITE(5300,*) '=========================================================='
   WRITE(5300,*) 'bfgs_b2: Starting next point search: '
   WRITE(5300,*) '=========================================================='

   converged = .FALSE.
   DO WHILE (.NOT.converged)

      WRITE(5300,*) 'Point A: '
      WRITE(5300,*) 'prelimNextPoint: '
      WRITE(5300,'(5f20.10)') prelimNextPoint(:)
      WRITE(5300,*) 'prelimGradient: '
      WRITE(5300,'(5f20.10)') prelimGradient(:)

      ! Determine free parameters (parameters not on boundaries with a gradient pointing outwards)
      numFixedParams = 0
      DO i = 1, vecLength
         IF(enabledConstraints(i)) THEN
            IF((prelimNextPoint(i).LE.minConstraints(i)+eps).AND.(prelimGradient(i).GT.0.0)) THEN
               freeParams(i) = .FALSE.
               numFixedParams = numFixedParams + 1
            END IF
            IF((prelimNextPoint(i).GE.maxConstraints(i)-eps).AND.(prelimGradient(i).LT.0.0)) THEN
               freeParams(i) = .FALSE.
               numFixedParams = numFixedParams + 1
            END IF
         END IF
      END DO
      numFreeParams = vecLength - numFixedParams

      ! Perform basis transformation to obtain inverse Hessian matrix for the free parameter subspace

      ALLOCATE (subspaceBasis(vecLength,numFreeParams))
      subspaceBasis = 0.0
      subspaceGrad = 0.0
      j = 1
      DO i = 1, vecLength
         IF(freeParams(i)) THEN
            subspaceBasis(i,j) = 1.0
            subspaceGrad(j) = prelimGradient(i)
            j = j + 1
         END IF
      END DO

      ALLOCATE (tempMatC(vecLength,numFreeParams))
      ALLOCATE (subspaceInvHessMat(numFreeParams,numFreeParams))
      tempMatC(:,:) = 0.0
      subspaceInvHessMat(:,:) = 0.0

      WRITE(5000,*) '========================================================================='
      WRITE(5000,*) 'prelimNextPoint: '
      WRITE(5000,'(5f20.10)') prelimNextPoint(:)
      WRITE(5000,*) 'prelimGradient: '
      WRITE(5000,'(5f20.10)') prelimGradient(:)
      WRITE(5000,*) 'subspaceBasis: '
      DO i = 1, SIZE(subspaceBasis,2)
         WRITE(5000,'(5f20.10)') subspaceBasis(:,i)
      END DO

      CALL dgemm('N','N',vecLength,numFreeParams,vecLength,1.0,invHessMat(:,:),vecLength,subspaceBasis(:,:),vecLength,0.0,tempMatC(:,:),vecLength)
      CALL dgemm('N','N',numFreeParams,numFreeParams,vecLength,1.0,TRANSPOSE(subspaceBasis(:,:)),numFreeParams,tempMatC(:,:),vecLength,0.0,subspaceInvHessMat(:,:),numFreeParams)

      DEALLOCATE(tempMatC)
      DEALLOCATE(subspaceBasis)

      ! Calculate next point (ignoring constraints)

      tempVecA(:) = 0.0

      WRITE(5000,*) 'numFreeParams: ', numFreeParams
      WRITE(5000,*) 'subspaceInvHessMat: '
      DO i = 1, SIZE(subspaceInvHessMat,2)
         WRITE(5000,'(5f20.10)') subspaceInvHessMat(:,i)
      END DO
      WRITE(5000,*) '-subspaceGrad: '
      WRITE(5000,'(5f20.10)') -subspaceGrad(:numFreeParams)
      WRITE(5000,*) 'tempVecA: '
      WRITE(5000,'(5f20.10)') tempVecA(:numFreeParams)

      CALL dgemv('N',numFreeParams,numFreeParams,1.0,subspaceInvHessMat(:,:),numFreeParams,-subspaceGrad(:numFreeParams),1,0.0,tempVecA(:numFreeParams),1)

      WRITE(5000,*) '----------------------------'
      WRITE(5000,*) 'tempVecA: '
      WRITE(5000,'(5f20.10)') tempVecA(:numFreeParams)
      WRITE(5000,*) 'prelimNextPoint: '
      WRITE(5000,'(5f20.10)') prelimNextPoint(:)

      !!! TEST START
!      temp = SQRT(ddot(numFreeParams,tempVecA(:numFreeParams),1,tempVecA(:numFreeParams),1))
      temp = MAXVAL(ABS(tempVecA(:numFreeParams-1)))
      maxDisp = 0.05
      IF (temp.GT.maxDisp) tempVecA(:numFreeParams) = tempVecA(:numFreeParams) / (temp/maxDisp)
      dampingFactor = 0.5
      !!! TEST END

      j = 1
      DO i = 1, vecLength
         IF(freeParams(i)) THEN
            prelimNextPoint(i) = prelimNextPoint(i) + dampingFactor * tempVecA(j)
!            prelimNextPoint(i) = prelimNextPoint(i) + tempVecA(j)
            j = j + 1
         END IF
      END DO

      WRITE(5000,*) 'prelimNextPoint: '
      WRITE(5000,'(5f20.10)') prelimNextPoint(:)

      DEALLOCATE(subspaceInvHessMat)

      ! Constrain next point to box and check convergence
      converged = .TRUE.
      DO i = 1, vecLength
         IF (enabledConstraints(i).AND.freeParams(i)) THEN
            IF(prelimNextPoint(i).LT.minConstraints(i)-eps) THEN
!               converged = .FALSE.
            END IF
            IF(prelimNextPoint(i).GT.maxConstraints(i)+eps) THEN
!               converged = .FALSE.
            END IF
            IF(prelimNextPoint(i).LT.minConstraints(i)) THEN
               prelimNextPoint(i) = minConstraints(i)
            END IF
            IF(prelimNextPoint(i).GT.maxConstraints(i)) THEN
               prelimNextPoint(i) = maxConstraints(i)
            END IF
         END IF
      END DO

      WRITE(5000,*) 'prelimNextPoint: '
      WRITE(5000,'(5f20.10)') prelimNextPoint(:)
      WRITE(5000,*) 'prelimGradient: '
      WRITE(5000,'(5f20.10)') prelimGradient(:)

      ! If not converged calculate prelimGradient for next point according to quadratic model
      IF(.NOT.converged) THEN
         diffPoint(:) = prelimNextPoint(:) - parameters(:)
         diffGrad(:) = 0.0
         CALL dgemv('N',vecLength,vecLength,1.0,hessMat(:,:),vecLength,diffPoint(:),1,0.0,diffGrad(:),1)
         prelimGradient(:) = prelimGradient(:) + diffGrad(:)
      END IF

      ! Test: Set the gradient entry for the Lagrange multiplier explicitly
      ! (I hope this makes the approach more stable. The here introduced data is not automatically obtained)
      temp = ddot(vecLength,equalityLinCombi(:vecLength),1,prelimNextPoint(:vecLength),1)
      prelimGradient(vecLength) = temp - equalityCriterion   ! derivative with respect to the Lagrange parameter

      WRITE(5000,*) 'prelimGradient: '
      WRITE(5000,'(5f20.10)') prelimGradient(:)

   END DO

   ! 3. Project prelimNextPoint onto the allowed subspace given by the equality constraints
   ! Note: I only consider a single equality constraint so this part becomes simpler

   ! Determine free parameters (parameters not on boundaries with a gradient pointing outwards)

!   WRITE(5200,*) '=========================================================='
!   WRITE(5200,*) 'bfgs_b2: Starting charge constraint ensurance: '
!   WRITE(5200,*) '=========================================================='

   converged = .FALSE.
   DO WHILE (.NOT.converged)

!      WRITE(5200,*) 'Point A: '
!      WRITE(5200,*) 'prelimNextPoint: '
!      WRITE(5200,'(5f20.10)') prelimNextPoint(:)

      freeParams(:) = .TRUE.
      numFixedParams = 0
      DO i = 1, vecLength
         IF(enabledConstraints(i)) THEN
            IF(prelimNextPoint(i).LT.minConstraints(i) - eps) THEN
               freeParams(i) = .FALSE.
               numFixedParams = numFixedParams + 1
            END IF
            IF(prelimNextPoint(i).GT.maxConstraints(i) + eps) THEN
               freeParams(i) = .FALSE.
               numFixedParams = numFixedParams + 1
            END IF
            IF(prelimNextPoint(i).LT.minConstraints(i)) THEN
               prelimNextPoint(i) = minConstraints(i)
            END IF
            IF(prelimNextPoint(i).GT.maxConstraints(i)) THEN
               prelimNextPoint(i) = maxConstraints(i)
            END IF
         END IF
      END DO
      numFreeParams = vecLength - numFixedParams

      subspaceNextPoint(:) = 0.0
      subspaceLinCombi(:) = 0.0
      subspaceEqualityCrit = equalityCriterion
      j = 1
      DO i = 1, vecLength
         IF(freeParams(i)) THEN
            subspaceNextPoint(j) = prelimNextPoint(i)
            subspaceLinCombi(j) = equalityLinCombi(i)
            j = j + 1
         ELSE
            subspaceEqualityCrit = subspaceEqualityCrit - prelimNextPoint(i) * equalityLinCombi(i)
         END IF
      END DO

      temp = ddot(numFreeParams,subspaceLinCombi(:),1,subspaceNextPoint(:),1)
      converged = .TRUE.
      IF(ABS(temp - subspaceEqualityCrit).GT.largeEps) THEN
         converged = .FALSE.

!         WRITE(5100,*) 'subspaceNextPoint(:): '
!         WRITE(5100,*) subspaceNextPoint(:)
!         WRITE(5100,*) 'subspaceEqualityCrit: ', subspaceEqualityCrit
!         WRITE(5100,*) 'temp: ', temp

         subspaceNextPoint(:) = subspaceNextPoint(:) * subspaceEqualityCrit / temp
      END IF
      j = 1
      DO i = 1, vecLength
         IF(freeParams(i)) THEN
            prelimNextPoint(i) = subspaceNextPoint(j)
            j = j + 1
         END IF
      END DO

      ! Also ensure that the equality criterion is fulfilled for the whole space.
      ! (This is only for pathological cases)

      temp = ddot(vecLength,equalityLinCombi(:vecLength),1,prelimNextPoint(:vecLength),1)
      IF(ABS(temp - equalityCriterion).GT.largeEps) THEN
         converged = .FALSE.
         prelimNextPoint(:vecLength) = prelimNextPoint(:vecLength) * equalityCriterion / temp
      END IF

   END DO

   ! 4. Update lastParameters, lastGradient, parameters

   lastParameters(:) = parameters(:)
   lastGradient(:) = gradient(:)
   parameters(:) = prelimNextPoint(:)
!   parameters(:) = 0.25*(prelimNextPoint(:) - lastParameters(:)) + lastParameters(:)

   DEALLOCATE(tempMatB)
   DEALLOCATE(invHessMat)
   DEALLOCATE(hessMat)

   WRITE(*,*) 'bfgs_b2 (output): parameters(:vecLength): ', parameters(:vecLength)

END SUBROUTINE bfgs_b2

END MODULE m_bfgs_b2

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR