!> This module is basically a store for the input parameters that are specified in the namelists \a knobs and \a parameters. In general, the names of the public variables in this module are the same as the name of the input parameter they correspond to.
module run_parameters
  implicit none

  private

  public :: init_run_parameters, finish_run_parameters
  public :: check_run_parameters, wnml_run_parameters
  public :: write_trinity_parameters
  public :: beta, zeff, tite, reset, immediate_reset
  public :: fphi, fapar, fbpar
  public :: code_delt_max, wunits, woutunits, tunits
  public :: nstep, wstar_units, eqzip, margin
  public :: secondary, tertiary, harris
  public :: ieqzip, k0
  public :: delt_option_switch, delt_option_hand
  public :: delt_option_auto, vnm_init
  public :: avail_cpu_time, margin_cpu_time
  public :: rhostar, neo_test
  public :: do_eigsolve
  !> If true and nonlinear_mode is "off", return
  !! simple diffusive estimates of fluxes to trinity
  public :: trinity_linear_fluxes
  !> If true and nonlinear_mode is "off",
  !! and trinity_linear_fluxes is true, return
  !! quasilinear estimates of fluxes to trinity
  public :: trinity_ql_fluxes
  public :: user_comments

  !> If true use old diagnostics. Included for testing
  !! only and will eventually be removed. Please use
  !! the new diagnostics module!
  public :: use_old_diagnostics

  public :: set_overrides

  real :: beta, zeff, tite
  real :: fphi, fapar, fbpar, faperp
  real :: delt, code_delt_max, user_delt_max, margin
  real, dimension (:), allocatable :: wunits, woutunits, tunits
  real, dimension (2) :: vnm_init
  real :: avail_cpu_time, margin_cpu_time
  integer :: nstep
  logical :: reset=.false.
  logical :: immediate_reset
  logical :: wstar_units, eqzip
  logical :: secondary, tertiary, harris
  real :: k0
  integer :: delt_option_switch
  integer, parameter :: delt_option_hand = 1, delt_option_auto = 2
  logical :: initialized = .false.
  logical :: rpexist, knexist
  real :: rhostar
  logical :: neo_test
  logical :: trinity_linear_fluxes, do_eigsolve
  logical :: trinity_ql_fluxes
  logical :: use_old_diagnostics

  character(len=100000) :: user_comments

  integer, allocatable :: ieqzip(:,:)
  integer :: eqzip_option_switch
  integer, parameter :: &
       eqzip_option_none = 1, &
       eqzip_option_secondary = 2, &
       eqzip_option_tertiary = 3, &
       eqzip_option_equilibrium = 4

contains
  subroutine check_run_parameters(report_unit)
    implicit none
    integer, intent(in) :: report_unit
    if (fphi /= 1.) then
       write (report_unit, *)
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('fphi in the knobs namelist = ',e11.4)") fphi
       write (report_unit, fmt="('fphi is a scale factor of all instances of Phi (the electrostatic potential).')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *)
    end if

    if (fapar == 0.) then
       write (report_unit, fmt="('A_parallel will not be included in the calculation.')")
    end if
    if (fapar == 1.) then
       write (report_unit, fmt="('A_parallel will be included in the calculation.')")
    end if
    if (fapar /= 0. .and. fapar /= 1.) then
       write (report_unit, *)
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('fapar in the knobs namelist = ',e11.4)") fapar
       write (report_unit, fmt="('fapar is a scale factor of all instances of A_parallel (the parallel vector potential).')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *)
    end if

    if (fbpar == 0.) then
       write (report_unit, fmt="('B_parallel will not be included in the calculation.')")
    end if
    if (fbpar == 1.) then
       write (report_unit, fmt="('B_parallel will be included in the calculation.')")
    end if
    if (fbpar /= 0. .and. fbpar /= 1.) then
       write (report_unit, *)
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('fbpar in the knobs namelist = ',e11.4)") fbpar
       write (report_unit, fmt="('fbpar is a scale factor of all instances of B_parallel &
           & (the perturbed parallel magnetic field).')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *)
    end if

    if (eqzip) then
       write (report_unit, *)
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('eqzip = T in the knobs namelist.')")
       write (report_unit, fmt="('This freezes some modes in time for a secondary stability analysis.')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *)
       if (secondary) write (report_unit, fmt="('Mode with kx = 0, ky = ky_min fixed in time')")
       if (tertiary)  write (report_unit, fmt="('Mode with ky = 0, kx = kx_min fixed in time')")
    end if

    write (report_unit, *)
    if(immediate_reset)then
       write (report_unit, fmt="('The time step will be reset immediately after cfl violation detected.')")
    else
       write (report_unit, fmt="('The time step will be reset just before the next time step after cfl violation detected.')")
    endif
    write (report_unit, *)
  end subroutine check_run_parameters

  subroutine wnml_run_parameters(unit,electrons,collisions)
    implicit none
    integer, intent(in) :: unit
    logical, intent(in) :: electrons, collisions
    if (rpexist) then
       write (unit, *)
       write (unit, fmt="(' &',a)") "parameters"
       write (unit, fmt="(' beta = ',e17.10)") beta       ! if zero, fapar, fbpar should be zero
       if (collisions) write (unit, fmt="(' zeff = ',e17.10)") zeff
       if (.not. electrons)  write (unit, fmt="(' tite = ',e17.10)") tite
!CMR, 10/2/2011: zip not in this namelist, so removing it!
!       if (zip) write (unit, fmt="(' zip = ',L1)") zip
       write (unit, fmt="(' /')")
    endif
    if (knexist) then
       write (unit, *)
       write (unit, fmt="(' &',a)") "knobs"
       write (unit, fmt="(' fphi   = ',f6.3)") fphi
       write (unit, fmt="(' fapar  = ',f6.3)") fapar
       write (unit, fmt="(' fbpar = ',f6.3)") fbpar
       write (unit, fmt="(' delt = ',e17.10)") delt
       write (unit, fmt="(' nstep = ',i8)") nstep
       write (unit, fmt="(' wstar_units = ',L1)") wstar_units
       if (eqzip) then
          write (unit, fmt="(' eqzip = ',L1)") eqzip
          write (unit, fmt="(' secondary = ',L1)") secondary
          write (unit, fmt="(' tertiary = ',L1)") tertiary
       end if
       write (unit, fmt="(' margin = ',e17.10)") margin
       select case (delt_option_switch)
       case (delt_option_auto)
          write (unit, fmt="(' delt_option = ',a)") '"check_restart"'
       case (delt_option_hand)
          ! nothing
       end select
       write (unit, fmt="(' immediate_reset = ',L1)") immediate_reset
       write (unit, fmt="(' /')")
    endif
  end subroutine wnml_run_parameters

  subroutine init_run_parameters
    use kt_grids, only: init_kt_grids, naky, nakx => ntheta0
    use gs2_time, only: init_delt, user2code
    implicit none

    if (initialized) return
    initialized = .true.

    call read_parameters

    call init_kt_grids
    call init_delt (delt)
    call user2code (user_delt_max, code_delt_max)

    if(.not. allocated(wunits)) allocate (wunits(naky))
    if(.not. allocated(woutunits)) allocate (woutunits(naky))
    if(.not. allocated(tunits)) allocate (tunits(naky))

! omega_* normalization of time:
    call adjust_time_norm

    if(.not.allocated(ieqzip)) allocate(ieqzip(nakx,naky))
    ieqzip(1:nakx,1:naky)=1
    select case (eqzip_option_switch)
    case (eqzip_option_secondary)
       ! suppress evolution of secondary mode
       ieqzip(1,2) = 0
    case (eqzip_option_tertiary)
       ! suppress evolution of tertiary mode
       ieqzip(2,1) = 0
       ieqzip(nakx,1) = 0
    case (eqzip_option_equilibrium)
       ! suppress evolution of 1D equilibrium (x dependent)
       ieqzip(1:nakx,1) = 0
    end select
  end subroutine init_run_parameters

  subroutine read_parameters
    use file_utils, only: input_unit, error_unit, input_unit_exist
    use mp, only: proc0, broadcast
    use gs2_save, only: init_dt, init_vnm
    use text_options, only: text_option, get_option_value
    use kt_grids, only: gryfx
    implicit none
    type (text_option), dimension (4), parameter :: eqzipopts = &
         (/ text_option('none', eqzip_option_none), &
            text_option('secondary', eqzip_option_secondary), &
            text_option('tertiary', eqzip_option_tertiary), &
            text_option('equilibrium', eqzip_option_equilibrium) /)
    character (len=20) :: eqzip_option
    type (text_option), dimension (3), parameter :: deltopts = &
         (/ text_option('default', delt_option_hand), &
            text_option('set_by_hand', delt_option_hand), &
            text_option('check_restart', delt_option_auto) /)
    character(20) :: delt_option
    integer :: ierr, istatus, in_file
    real :: delt_saved
    real, dimension (2) :: vnm_saved

    real :: teti  ! for back-compatibility
    namelist /parameters/ beta, zeff, tite, teti, k0, rhostar, user_comments
    namelist /knobs/ fphi, fapar, fbpar, delt, nstep, wstar_units, eqzip, &
         delt_option, margin, secondary, tertiary, faperp, harris, &
!         avail_cpu_time, eqzip_option, include_lowflow, neo_test
         avail_cpu_time, margin_cpu_time, eqzip_option, neo_test, &
         trinity_linear_fluxes, do_eigsolve, immediate_reset, &
         use_old_diagnostics,  trinity_ql_fluxes

    if (proc0) then
       beta = 0.0
       zeff = 1.0
       tite = 1.0
       teti = -100.0
       rhostar = 3.e-3
       user_comments = ''
       k0 = 1.

       fbpar = -1.0
       faperp = 0.0

       ! EGH added these defaults Feb 2018
       delt = 0.1
       nstep = 100

!       include_lowflow = .false.
       neo_test = .false.
       wstar_units = .false.
       eqzip_option = 'none'
       eqzip = .false.
       secondary = .true.
       tertiary = .false.
       harris = .false.
       delt_option = 'default'
       margin = 0.05
       avail_cpu_time = 1.e10
       margin_cpu_time = 300.
       trinity_linear_fluxes = .false.
       trinity_ql_fluxes = .false.
       do_eigsolve = .false.
       immediate_reset = .true.
       use_old_diagnostics = .false.

       in_file = input_unit_exist("parameters", rpexist)
!       if (rpexist) read (unit=input_unit("parameters"), nml=parameters)
       if (rpexist) read (unit=in_file,nml=parameters)

       in_file = input_unit_exist("knobs", knexist)
!       if (knexist) read (unit=input_unit("knobs"), nml=knobs)
       if (knexist) read (unit=in_file, nml=knobs)

       if (teti /= -100.0) tite = teti

! Allow faperp-style initialization for backwards compatibility.
! Only fbpar is used outside of this subroutine.
       if (fbpar == -1.) then
          fbpar = faperp
       end if

! Override fapar, fbpar and beta if set inconsistently:

       !> If we have zero beta then disable apar and bpar fields
       !! as these contribute nothing to result but slow down calculation.
       if(beta.eq.0) then
          if(((fapar.ne.0) .or. (fbpar.ne.0)).and. proc0) then
             ierr = error_unit()
             write(ierr,'("Warning: Disabling apar and bpar as beta = 0.")')
          endif
          fapar = 0.
          fbpar = 0.
       endif

       !> If we have non-zero beta then alert the user that they have some
       !! of the perturbed magnetic fields disabled. This may be intended
       !! behaviour, so we throw a warning. However, if both are disabled,
       !! we set beta=0. This will make the simulation electrostatic, as
       !! probably intended if both apar and bpar are set to zero.
       if(fapar.eq.0 .and. fbpar.eq.0 .and. beta.ne.0) then
          if(proc0) then
             ierr = error_unit()
             write(ierr,'("Warning: Both fapar and fbpar are zero: setting beta = 0.")')
          endif
          beta = 0.0
       endif
       if((beta.ne.0).and.proc0) then
          ierr = error_unit()
          if(fapar.eq.0) write(ierr,'("Warning: Running with finite beta but fapar=0.")')
          if(fbpar.eq.0) write(ierr,'("Warning: Running with finite beta but fbpar=0.")')
       endif

       if (eqzip) then
          if (secondary .and. tertiary) then
             ierr = error_unit()
             write (ierr, *) 'Forcing secondary = FALSE'
             write (ierr, *) 'because you have chosen tertiary = TRUE'
             secondary = .false.
          end if
          if (secondary .and. harris) then
             ierr = error_unit()
             write (ierr, *) 'Forcing secondary = FALSE'
             write (ierr, *) 'because you have chosen harris = TRUE'
             secondary = .false.
          end if
          if (tertiary .and. harris) then
             ierr = error_unit()
             write (ierr, *) 'Forcing tertiary = FALSE'
             write (ierr, *) 'because you have chosen harris = TRUE'
             tertiary = .false.
          end if
       endif

       ierr = error_unit()
       call get_option_value &
            (delt_option, deltopts, delt_option_switch, ierr, &
            "delt_option in knobs",.true.)

       call get_option_value ( &
            eqzip_option, eqzipopts, eqzip_option_switch, error_unit(), &
            "eqzip_option in knobs",.true.)

!!$       ! eqzip_option replaces eqzip, secondary, tertiary, harris
!!$       if (eqzip .and. eqzip_option_switch == eqzip_option_none) then
!!$          if (harris) then
!!$             eqzip_option_switch = eqzip_option_equilibrium
!!$             write(error_unit(),*) 'eqzip_option is set to equilibrium'
!!$          else
!!$             if (tertiary) then
!!$                eqzip_option_switch = eqzip_option_tertiary
!!$                write(error_unit(),*) 'eqzip_option is set to tertiary'
!!$             else
!!$                eqzip_option_switch = eqzip_option_secondary
!!$                write(error_unit(),*) 'eqzip_option is set to secondary'
!!$             end if
!!$          end if
!!$       else if (eqzip .and. eqzip_option_switch /= eqzip_option_none) then
!!$          write(error_unit(),*) 'eqzip, secondary, tertiary, harris are ignored'
!!$          write(error_unit(),*) 'because eqzip_option exists'
!!$       end if

    end if

    ! FOR GRYFX
    if (gryfx()) then
      delt = sqrt(2.0)*delt/2.0
      nstep = nstep*2
    end if

    call broadcast (delt_option_switch)
    call broadcast (delt)
    call broadcast (beta)
    call broadcast (zeff)
    call broadcast (tite)
    call broadcast (fphi)
    call broadcast (fapar)
    call broadcast (fbpar)
    call broadcast (nstep)
    call broadcast (wstar_units)
    call broadcast (eqzip)
    call broadcast (secondary)
    call broadcast (tertiary)
    call broadcast (harris)
    call broadcast (margin)
    call broadcast (k0)
    call broadcast (avail_cpu_time)
    call broadcast (margin_cpu_time)
    call broadcast (eqzip_option_switch)
!    call broadcast (include_lowflow)
    call broadcast (rhostar)
    call broadcast (neo_test)
    call broadcast (trinity_linear_fluxes)
    call broadcast (trinity_ql_fluxes)
    call broadcast (do_eigsolve)
    call broadcast (immediate_reset)
    call broadcast (use_old_diagnostics)

    user_delt_max = delt

    delt_saved = delt
    if (delt_option_switch == delt_option_auto) then
       vnm_init = 1.0
       call init_vnm (vnm_saved, istatus)
       if (istatus == 0) vnm_init = vnm_saved
       call init_dt (delt_saved, istatus)
       if (istatus == 0) delt  = delt_saved
    endif

  end subroutine read_parameters

  subroutine write_trinity_parameters(trinpars_unit)
    integer, intent(in) :: trinpars_unit
    write(trinpars_unit, "(A15)") '&run_parameters'
    write (trinpars_unit, *) ' beta = ', beta
    write (trinpars_unit, "(A1)") '/'

  end subroutine write_trinity_parameters

  subroutine set_overrides(tstep_ov)
    use overrides, only: timestep_overrides_type
    type(timestep_overrides_type), intent(in) :: tstep_ov
    if (tstep_ov%override_immediate_reset) immediate_reset=tstep_ov%immediate_reset
  end subroutine set_overrides

  subroutine adjust_time_norm
    use file_utils, only: error_unit
    use mp, only: proc0
    use kt_grids, only: aky
    implicit none
!CMR: Sep 2010
! Attempt to understand time normalisation variables, which are arrays(naky)
!    TUNITS: DT(KY)=TUNITS(KY).CODE_DT
!            This is a generally very useful variable to store ky dependent
!            timestep in the code time normalisation.
!            Used to multiply ky independent source terms on RHS of GKE.
!    WUNITS: WUNITS(KY)=AKY(KY)*TUNITS(KY)/2
!            Auxiliary variable.  Used to save compute operations when
!            evaluating source terms on RHS of GKE that are proportional to ky.
!            !! The Mysterious factor 1/2 Explained !!
!            The factor of 1/2 arises because those source terms were first
!            specified using the normalisation Tref=mref vtref^2
! [R Numata et al, "AstroGK: Astrophysical gyrokinetics code", JCP, 2010].
!CMRend
    if (wstar_units) then
       wunits = 1.0
       where (aky /= 0.0)
          tunits = 2.0/aky
       elsewhere
          tunits = 0.0
       end where
       if (any(tunits == 0.0) .and. proc0) then
          write (error_unit(), *) &
               "WARNING: wstar_units=.true. and aky=0.0: garbage results"
          print *, &
               "WARNING: wstar_units=.true. and aky=0.0: garbage results"
       end if
!CMR: Sep 2010
!  Changes to allow wstar_units to be used
!CMRend
    else
       tunits = 1.0
       wunits = aky/2.0
    end if
    woutunits = 1.0/tunits
  end subroutine adjust_time_norm

  subroutine finish_run_parameters
    implicit none

    if (allocated(wunits)) deallocate (wunits, woutunits, tunits)
    if (allocated(ieqzip)) deallocate (ieqzip)

    initialized = .false.
  end subroutine finish_run_parameters
end module run_parameters


text:
      PROGRAM WRFILE
      CHARACTER FNAME*80
      WRITE(*,*) 'Enter Filename'
      READ(*,*) FNAME
C Use of NAME in OPEN is VMS Fortran extension
      OPEN (UNIT=8,NAME=FNAME,ERR=999)
      WRITE(8,100) FNAME
 100  FORMAT(1x,'Sample Output to File ',A80)
      STOP
 999  WRITE(*,*) 'Error opening file'
      STOP 2
      END

ext: f

max_stars_repo_name: kurtsansom/ftnchek

max_issues_repo_name: kurtsansom/ftnchek

max_forks_repo_name: kurtsansom/ftnchek


text:
MODULE constants_NSWC
! Contains the NSWC functions SPMPAR, DPMPAR, EPSLN, DEPSLN, EXPARG & DXPARG
!-----------------------------------------------------------------------
!     WRITTEN using F90 intrinsics by
!        Alan Miller
!        CSIRO Mathematical & Information Sciences
!        CLAYTON, VICTORIA, AUSTRALIA 3169
!     Latest revision - 1 February 1997
!-----------------------------------------------------------------------

IMPLICIT NONE
INTEGER, PARAMETER     :: dp = SELECTED_REAL_KIND(15, 60)

CONTAINS

FUNCTION spmpar (i) RESULT(fn_val)
!-----------------------------------------------------------------------

!     SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR
!     THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
!     I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
!     SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
!     ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN

!        SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,

!        SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,

!        SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
!-----------------------------------------------------------------------

IMPLICIT NONE
INTEGER, INTENT(IN) :: i
REAL                :: fn_val

! Local variable
REAL                :: one = 1.0

SELECT CASE (i)
  CASE (1)
    fn_val = EPSILON(one)
  CASE (2)
    fn_val = TINY(one)
  CASE (3)
    fn_val = HUGE(one)
END SELECT

RETURN
END FUNCTION spmpar

FUNCTION dpmpar (i) RESULT(fn_val)
!-----------------------------------------------------------------------

!     DPMPAR PROVIDES THE DOUBLE PRECISION MACHINE CONSTANTS FOR
!     THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
!     I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
!     DOUBLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
!     ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN

!        DPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,

!        DPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,

!        DPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
!-----------------------------------------------------------------------

IMPLICIT NONE
INTEGER, INTENT(IN) :: i
REAL (dp)           :: fn_val

! Local variable
REAL (dp)    :: one = 1._dp

SELECT CASE (i)
  CASE (1)
    fn_val = EPSILON(one)
  CASE (2)
    fn_val = TINY(one)
  CASE (3)
    fn_val = HUGE(one)
END SELECT

RETURN
END FUNCTION dpmpar

FUNCTION epsln () RESULT(fn_val)
!--------------------------------------------------------------------
!     THE EVALUATION OF LN(EPS) WHERE EPS IS THE SMALLEST NUMBER
!     SUCH THAT 1.0 + EPS .GT. 1.0 .  L IS A DUMMY ARGUMENT.
!--------------------------------------------------------------------
IMPLICIT NONE
REAL                :: fn_val

! Local variable
REAL                :: one = 1.0

fn_val = LOG( EPSILON(one) )
RETURN
END FUNCTION epsln

FUNCTION exparg (l) RESULT(fn_val)
!--------------------------------------------------------------------
!     IF L = 0 THEN  EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
!     EXP(W) CAN BE COMPUTED.
!
!     IF L IS NONZERO THEN  EXPARG(L) = THE LARGEST NEGATIVE W FOR
!     WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO.
!
!     NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED.
!--------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: l
REAL                :: fn_val

! Local variable
REAL                :: one = 1.0

IF (l == 0) THEN
  fn_val = LOG( HUGE(one) )
ELSE
  fn_val = LOG( TINY(one) )
END IF
RETURN
END FUNCTION exparg

FUNCTION depsln () RESULT(fn_val)
!--------------------------------------------------------------------
!     THE EVALUATION OF LN(EPS) WHERE EPS IS THE SMALLEST NUMBER
!     SUCH THAT 1.D0 + EPS .GT. 1.D0 .  L IS A DUMMY ARGUMENT.
!--------------------------------------------------------------------
IMPLICIT NONE
REAL (dp)           :: fn_val

! Local variable
REAL (dp)    :: one = 1._dp

fn_val = LOG( EPSILON(one) )
RETURN
END FUNCTION depsln

FUNCTION dxparg (l) RESULT(fn_val)
!--------------------------------------------------------------------
!     IF L = 0 THEN  DXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
!     DEXP(W) CAN BE COMPUTED.
!
!     IF L IS NONZERO THEN  DXPARG(L) = THE LARGEST NEGATIVE W FOR
!     WHICH THE COMPUTED VALUE OF DEXP(W) IS NONZERO.
!
!     NOTE... ONLY AN APPROXIMATE VALUE FOR DXPARG(L) IS NEEDED.
!--------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: l
REAL (dp)           :: fn_val

! Local variable
REAL (dp)    :: one = 1._dp

IF (l == 0) THEN
  fn_val = LOG( HUGE(one) )
ELSE
  fn_val = LOG( TINY(one) )
END IF
RETURN
END FUNCTION dxparg

END MODULE constants_NSWC

SUBROUTINE smplx (a, b0, c, ka, m, n0, ind, ibasis, x, z, iter, mxiter,   &
                  numle, numge, bi, rerr)
!-----------------------------------------------------------------------
!     SIMPLEX PROCEDURE FOR SOLVING LINEAR PROGRAMMING PROBLEMS
!-----------------------------------------------------------------------
! Finds non-negative x's to maximize:
!     c1.x1 + ... + cn.xn
!
! Subject to the constraints:
!     a11.x1 + ... + a1n.xn {<=,=,>=} b1
!              ...
!     am1.x1 + ... + amn.xn {<=,=,>=} bm
!
! The constraints must be ordered so that the <= constraints, if any, come
! first, then any >= constraints, with the equality constraints last.
!
! Arguments:
! a(ka,n0)   INPUT   coefficients in the constraints.
! b0(m)      INPUT   right-hand side of the constraints.
! c(n0)      INPUT   vector of `costs' in the objective function.
! ka         INPUT   first dimension of array a.
! m          INPUT   dimension of array b0.
! n0         INPUT   2nd dimension of array a.
! ind        IN/OUT  If IND = 0 on input, the routine selects its own initial
!                    basis, else if IND = 1 then the indices of the initial
!                    basis should be in IBASIS.
!                    On output:
!                    IND = 0 the problem was solved
!                    IND = 1 the problem has no solution
!                    IND = 2 MXITER iterations were performed; more needed
!                    IND = 3 sufficient accuracy could not be maintained to
!                            solve the problem
!                    IND = 4 the problem has an unbounded solution
!                    IND = 5 input error detected
!                    IND = 6 the solution may have been obtained
! ibasis(m)  IN/OUT  indices of the variables in the basis.
! x()        OUTPUT  Dimension must be >= n + numle + numge.
!                    If IND = 0 or 6, it contains the values of the original,
!                    slack and surplus variables.
! z          OUTPUT  If IND = 0 or 6, contains the value of the objective.
! iter       OUTPUT  number of iterations used.
! mxiter     INPUT   maximum number of iterations.
! numle      INPUT   number of <= constraints.
! numge      INPUT   number of >= constraints.
! bi(m,m)    OUTPUT  the inverse of the basis matrix.
! rerr       OUTPUT  the estimated relative error achieved.
! N.B. The last two arguments in the NSWC routine have been eliminated.
!      These were work spaces.   They are now declared as automatic arrays.
!-----------------------------------------------------------------------
!     WRITTEN BY ALFRED H. MORRIS JR.
!        NAVAL SURFACE WEAPONS CENTER
!        DAHLGREN, VIRGINIA
!------------------------
!     INITIAL VERSION  DEC  1977
!     LAST UPDATE      OCT  1990
!-----------------------------------------------------------------------
!     Converted using F90 intrinsics by
!        Alan Miller
!        CSIRO Mathematical & Information Sciences
!        CLAYTON, VICTORIA, AUSTRALIA 3169
!     Latest revision - 5 February 1997
!-----------------------------------------------------------------------

USE constants_NSWC
IMPLICIT NONE

INTEGER, INTENT(IN)                    :: ka, m, n0, mxiter, numle, numge
INTEGER, INTENT(IN OUT)                :: ind
INTEGER, DIMENSION(:), INTENT(IN OUT)  :: ibasis
INTEGER, INTENT(OUT)                   :: iter
REAL (dp), DIMENSION(:,:), INTENT(IN)  :: a
REAL (dp), DIMENSION(:), INTENT(IN)    :: b0, c
REAL (dp), INTENT(OUT)                 :: z, rerr
REAL (dp), DIMENSION(:), INTENT(OUT)   :: x
REAL (dp), DIMENSION(:,:), INTENT(OUT) :: bi

INTERFACE
  SUBROUTINE smplx1 (a, b0, c, ka, m, n0, ind, ibasis, r, z, iter, mxiter,  &
                     eps0, rerrmn, rerrmx, rerr, numle, numge, bi)
    USE constants_NSWC
    IMPLICIT NONE
    INTEGER, INTENT(IN)                    :: ka, m, n0, mxiter, numle, numge
    INTEGER, INTENT(OUT)                   :: ind, iter
    REAL (dp), DIMENSION(:,:), INTENT(IN)  :: a
    REAL (dp), DIMENSION(:), INTENT(IN)    :: b0, c
    REAL (dp), DIMENSION(:,:), INTENT(OUT) :: bi
    REAL (dp), DIMENSION(:), INTENT(OUT)   :: r
    REAL (dp), INTENT(OUT)                 :: z, rerr
    REAL (dp), INTENT(IN)                  :: eps0, rerrmn, rerrmx
    INTEGER, DIMENSION(:), INTENT(IN OUT)  :: ibasis
  END SUBROUTINE smplx1
END INTERFACE

!------------------------
!     DIMENSION X(N0+NUMLE+NUMGE)
!------------------------

!     ********** EPS0 IS A MACHINE DEPENDENT PARAMETER. ASSIGN EPS0
!                THE VALUE U WHERE U IS THE SMALLEST POSITIVE FLOATING
!                POINT NUMBER SUCH THAT 1.0 + U .GT. 1.0.

!     Local variables
REAL (dp) :: eps0, rerrmn, rerrmx

eps0 = dpmpar(1)

!------------------------
rerrmn = 10.0_dp*eps0
rerrmx = 1.d-4
IF (eps0 < 1.d-13) rerrmx = 1.d-5

CALL smplx1(a, b0, c, ka, m, n0, ind, ibasis, x, z, iter, mxiter, eps0,   &
        rerrmn, rerrmx, rerr, numle, numge, bi)
RETURN
END SUBROUTINE smplx

SUBROUTINE smplx1 (a, b0, c, ka, m, n0, ind, ibasis, r, z, iter, mxiter,  &
                   eps0, rerrmn, rerrmx, rerr, numle, numge, bi)
!----------------------
!     NSTEP = 1   ELIMINATE THE NEGATIVE VARIABLES
!     NSTEP = 2   PHASE 1 OF THE SIMPLEX ALGORITHM
!     NSTEP = 3   PHASE 2 OF THE SIMPLEX ALGORITHM
!----------------------
!     MXITER = THE MAXIMUM NUMBER OF ITERATIONS PERMITTED
!     ITER = THE NUMBER OF THE CURRENT ITERATION
!     ICOUNT = THE NUMBER OF ITERATIONS SINCE THE LAST INVERSION
!----------------------
!     NUMLE = THE NUMBER OF .LE. CONSTRAINTS
!     NUMGE = THE NUMBER OF .GE. CONSTRAINTS
!----------------------
!     THE ROUTINE ASSUMES THAT THE .LE. CONSTRAINTS PRECEDE THE .GE.
!     CONSTRAINTS AND THAT THE .EQ. CONSTRAINTS COME LAST. THERE ARE
!     M CONSTRAINTS. X(N0+I) IS THE SLACK, SURPLUS, OR ARTIFICIAL
!     VARIABLE FOR THE I-TH CONSTRAINT (I=1, ..., M).
!----------------------
!     N0 = THE NUMBER OF ORGINAL VARIABLES
!     NS = THE NUMBER OF ORGINAL AND SLACK VARIABLES
!     N  = THE NUMBER OF ORGINAL, SLACK, AND SURPLUS VARIABLES
!     NUM = THE TOTAL NUMBER OF VARIABLES
!----------------------
!     RERRMN = THE SMALLEST RELATIVE ERROR TOLERANCE USED
!     RERRMX = THE LARGEST RELATIVE ERROR TOLERACE USED
!     RERR   = THE ESTIMATED CURRENT RELATIVE ERROR
!----------------------
!     ASSUME THAT
!         B0 = (B0(1), ..., B0(M))
!         C  = (C(1), ..., C(N0))
!         Z  = C(1)*X(1)+...+C(N0)*X(N0)
!     THE PROBLEM IS TO MAXIMIZE Z SUBJECT TO
!         AX(LE,EQ,GE)B0
!         X.GE.0
!----------------------
!     ON INPUT IND CAN HAVE THE VALUES
!         IND = 0   NO BEGINNING BASIS IS PROVIDED BY THE USER
!         IND = 1   THE ARRAY IBASIS HAS BEEN SET BY THE USER
!     ON OUTPUT IND IS ASSIGNED ONE OF THE VALUES
!         IND = 0   Z WAS SUCCESSFULLY MAXIMIZED
!         IND = 1   THE PROBLEM HAS NO FEASIBLE SOLUTION
!         IND = 2   MXITER ITERATIONS WERE PERFORMED
!         IND = 3   SUFFICIENT ACCURACY CANNOT BE MAINTAINED
!         IND = 4   THE PROBLEM HAS AN UNBOUNDED SOLUTION
!         IND = 5   THERE IS AN INPUT ERROR
!         IND = 6   Z WAS POSSIBLY MAXIMIZED
!----------------------
!     BASIS IS AN INTEGER ARRAY OF DIMENSION N0+M. FOR J.LE.N
!         BASIS(J) = 1  IF X(J) IS A BASIC VARIABLE
!         BASIS(J) = 0  IF X(J) IS NOT A BASIC VARIABLE
!     IF THE BASIC VARIABLES ARE X(I1), ..., X(IM) THEN
!         IBASIS = (I1, ..., IM)
!     ALSO XB(1), ..., XB(M) ARE THE CORRESPONDING VALUES OF THE
!     BASIC VARIABLES.
!----------------------
!     BI IS AN MXM ARRAY CONTAINING THE INVERSE OF THE BASIS MATRIX.
!----------------------
!     R IS AN ARRAY OF DIMENSION N. ON OUTPUT R CONTAINS THE CURRENT
!     VALUE OF X. DURING COMPUTATION R NORMALLY CONTAINS THE REDUCED
!     COSTS USED FOR THE SELECTION OF THE VARIABLE TO BE MADE BASIC.
!----------------------
USE constants_NSWC
IMPLICIT NONE
INTEGER, INTENT(IN)                    :: ka, m, n0, mxiter, numle, numge
INTEGER, INTENT(OUT)                   :: ind, iter
REAL (dp), DIMENSION(:,:), INTENT(IN)  :: a
REAL (dp), DIMENSION(:), INTENT(IN)    :: b0, c
REAL (dp), DIMENSION(:,:), INTENT(OUT) :: bi
REAL (dp), DIMENSION(:), INTENT(OUT)   :: r
REAL (dp), INTENT(OUT)                 :: z, rerr
REAL (dp), INTENT(IN)                  :: eps0, rerrmn, rerrmx
INTEGER, DIMENSION(:), INTENT(IN OUT)  :: ibasis
!----------------------

INTERFACE
  SUBROUTINE crout1(a, ka, n, iend, indx, temp, ierr)
    USE constants_NSWC
    IMPLICIT NONE
    INTEGER, INTENT(IN)                     :: ka, n, iend
    INTEGER, INTENT(OUT)                    :: ierr
    REAL (dp), DIMENSION(:), INTENT(IN OUT) :: a, temp
    INTEGER, DIMENSION(:), INTENT(IN OUT)   :: indx
  END SUBROUTINE crout1
END INTERFACE

!     Local variables
INTEGER   :: i, ibeg, icount, iend, ierr, ii, il, imin, iout, ip, j, jj,  &
             jmin, jp, k, ki, kj, l, ll, lrow, m0, mcheck, ms, n, npos,   &
             nrow, ns, nstep, num, bflag
REAL (dp) :: xb(m), y(m)
INTEGER   :: basis(m+n0), indx(m)
REAL (dp) :: amax, binorm, bmax, bmin, bnorm, cmin, const, eps, epsi, ratio, &
             rerr1, rmin, rtol, s, sgn, t, tol, total, w, xmax, zero = 0._dp, &
             dsum, dsump, dsumn, dt

!     ****** XMAX IS A MACHINE DEPENDENT CONSTANT. XMAX IS THE
!            LARGEST POSITIVE FLOATING POINT NUMBER.

xmax = dpmpar(3)

!----------------------
iter = 0
icount = 0
mcheck = MIN(5, 1 + m/15)
z = zero

!                CHECK FOR INPUT ERRORS

ms = numle + numge
IF (m < 2 .OR. n0 < 2 .OR. ms > m .OR. ka < m)  &
GO TO 12
DO i = 1, m
  IF (b0(i) < zero) GO TO 12
  xb(i) = zero
END DO
rtol = xmax
DO i = 1, n0
  IF (c(i) /= zero) rtol = MIN(ABS(c(i)), rtol)
END DO
rtol = rerrmx*rtol
GO TO 20

12 ind = 5
RETURN

!     FORMATION OF THE IBASIS AND BASIS ARRAYS.
!     (IF IND = 1  THEN THE IBASIS ARRAY IS DEFINED BY THE USER.)

20 ns = n0 + numle
n = ns + numge
IF (ind == 0) GO TO 30
num = n
DO i = 1, m
  IF (ibasis(i) > n) num = num + 1
END DO
GO TO 32
22 IF (ind == 0) GO TO 590
ind = 0

30 num = n0 + m
DO i = 1, m
  ibasis(i) = n0 + i
END DO
32 bflag = 0
basis(1:n) = 0
DO i = 1, m
  ki = ibasis(i)
  basis(ki) = 1
END DO
IF (ind == 1) GO TO 100

!          CALCULATION OF XB AND BI WHEN IND = 0

rerr = rerrmn
DO j = 1, m
  xb(j) = b0(j)
  bi(1:m,j) = zero
  bi(j,j) = 1.0_dp
END DO
IF (numge == 0) GO TO 630
jmin = numle + 1
DO j = jmin, ms
  xb(j) = -xb(j)
  bi(j,j) = -1.0_dp
END DO
GO TO 601

!                  REORDER THE BASIS

100 ibeg = 1
iend = m
DO i = 1, m
  IF (ibasis(i) > n0) THEN
    indx(ibeg) = ibasis(i)
    ibeg = ibeg + 1
  ELSE
    indx(iend) = ibasis(i)
    iend = iend - 1
  END IF
END DO
IF (iend == m) GO TO 22
ibasis(1:m) = indx(1:m)

!            REINVERSION OF THE BASIS MATRIX

DO j = 1, m
  kj = ibasis(j)
  IF (kj <= n0) GO TO 110
  IF (kj <= ns) GO TO 120
  IF (kj <= n) GO TO 130
  GO TO 120

  110 bi(1:m,j) = a(1:m,kj)
  CYCLE

  120 l = kj - n0
  bi(1:m,j) = zero
  bi(l,j) = 1.0_dp
  CYCLE

  130 l = kj - n0
  bi(1:m,j) = zero
  bi(l,j) = -1.0_dp
END DO

icount = 0
CALL crout1 (bi(1:,1), m, m, iend, indx, y, ierr)
IF (ierr /= 0) GO TO 580

!         CHECK THE ACCURACY OF BI AND RESET RERR

bnorm = zero
DO j = 1, m
  kj = ibasis(j)
  IF (kj <= n0) GO TO 140
  total = 1.0_dp
  GO TO 142
  140 total = SUM( ABS(a(1:m,kj)) )
  142 bnorm = MAX(bnorm, total)
END DO

binorm = zero
DO j = 1, m
  total = SUM( ABS(bi(1:m,j)) )
  binorm = MAX(binorm, total)
END DO
rerr = MAX(rerrmn, eps0*bnorm*binorm)
IF (rerr > 1.d-2) GO TO 580
bflag = 0

!                 RECALCULATION OF XB

DO i = 1, m
  dsump = zero
  dsumn = zero
  DO l = 1, m
    dt = bi(i,l)*b0(l)
    IF (dt > zero) THEN
      dsump = dsump + dt
    ELSE
      dsumn = dsumn + dt
    END IF
  END DO
  xb(i) = dsump + dsumn
  s = dsump
  t = dsumn
  tol = rerrmx*MAX(s, -t)
  IF (ABS(xb(i)) < tol) xb(i) = zero
END DO
GO TO 601

!     FIND THE NEXT VECTOR A(--, JP) TO BE INSERTED INTO THE BASIS

200 jp = 0
rmin = zero
IF (nstep == 3) rmin = -rtol
DO j = 1, n0
  IF (basis(j) /= 0) CYCLE
  IF (r(j) >= rmin) CYCLE
  jp = j
  rmin = r(j)
END DO
IF (n0 == n) GO TO 203
jmin = n0 + 1
rmin = rmin*1.1_dp
DO j = jmin, n
  IF (basis(j) /= 0) CYCLE
  IF (r(j) >= rmin) CYCLE
  jp = j
  rmin = r(j)
END DO
203 IF (jp /= 0) GO TO 300
IF (nstep < 2) THEN
   GO TO 800
ELSE IF (nstep == 2) THEN
  GO TO 230
ELSE
  GO TO 250
END IF

!     INSERT THE VALUES OF THE ORGINAL, SLACK, AND SURPLUS
!             VARIABLES INTO R, THEN TERMINATE.

220 r(1:n) = zero
DO i = 1, m
  ki = ibasis(i)
  IF (ki <= n) r(ki) = xb(i)
END DO
RETURN

!             COMPLETION OF THE NSTEP = 2 CASE

230 DO i = 1, m
  IF (ibasis(i) <= n) CYCLE
  IF (xb(i) > zero) GO TO 800
END DO
GO TO 680

240 IF (icount >= 5) GO TO 100
ind = 1
GO TO 220

!             COMPLETION OF THE NSTEP = 3 CASE

250 IF (rerr > 1.d-2) GO TO 251
ind = 0
GO TO 800
251 IF (icount >= 5) GO TO 100
ind = 6
GO TO 800

!     IF MXITER ITERATIONS HAVE NOT BEEN PERFORMED THEN BEGIN THE
!     NEXT ITERATION. COMPUTE THE JP-TH COLUMN OF BI*A AND STORE IT IN Y.

300 IF (iter < mxiter) GO TO 301
ind = 2
GO TO 220
301 iter = iter + 1
icount = icount + 1
IF (jp > ns) GO TO 330
IF (jp > n0) GO TO 320

nrow = 0
amax = zero
DO i = 1, m
  IF (a(i,jp) == zero) CYCLE
  nrow = nrow + 1
  indx(nrow) = i
  amax = MAX(ABS(a(i,jp)), amax)
END DO
IF (nrow /= 0) GO TO 310
ind = 4
GO TO 220

310 rerr1 = rerrmx*amax
DO i = 1, m
  dsum = zero
  DO ll = 1, nrow
    l = indx(ll)
    dsum = dsum + bi(i,l)*a(l,jp)
  END DO
  y(i) = dsum
  IF (ABS(y(i)) >= 5.d-3) CYCLE
  bmax = zero
  DO l = 1, m
    bmax = MAX(ABS(bi(i,l)), bmax)
  END DO
  tol = rerr1*bmax
  IF (ABS(y(i)) < tol) y(i) = zero
END DO
GO TO 350

320 l = jp - n0
y(1:m) = bi(1:m,l)
GO TO 350

330 l = jp - n0
y(1:m) = -bi(1:m,l)

350 DO i = 1, m
  IF (y(i) /= zero) GO TO 360
END DO
r(jp) = zero
iter = iter - 1
icount = icount - 1
GO TO 200

360 IF (nstep == 2) GO TO 430
IF (nstep > 2) GO TO 440

!     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC FOR THE NSTEP = 1 CASE

npos = 0
ip = 0
eps = zero
epsi = xmax
DO i = 1, m
  IF (xb(i) < zero .OR. y(i) <= zero) CYCLE
  ratio = xb(i)/y(i)
  IF (ratio < epsi) THEN
    epsi = ratio
    npos = 1
    indx(1) = i
    CYCLE
  ELSE IF (ratio > epsi) THEN
    CYCLE
  END IF
  npos = npos + 1
  indx(npos) = i
END DO
IF (npos == 0) GO TO 420
IF (epsi == zero) GO TO 460

DO i = 1, m
  IF (xb(i) >= zero .OR. y(i) >= zero) CYCLE
  ratio = xb(i)/y(i)
  IF (ratio > epsi) CYCLE
  IF (ratio < eps) CYCLE
  eps = ratio
  ip = i
END DO
IF (ip /= 0) GO TO 500
GO TO 460

420 DO i = 1, m
  IF (xb(i) >= zero .OR. y(i) >= zero) CYCLE
  ratio = xb(i)/y(i)
  IF (ratio < eps) CYCLE
  eps = ratio
  ip = i
END DO
GO TO 500

!     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC FOR THE NSTEP = 2 CASE

430 npos = 0
epsi = xmax
DO i = 1, m
  IF (y(i) <= zero) CYCLE
  ratio = xb(i)/y(i)
  IF (ratio < epsi) THEN
    epsi = ratio
    npos = 1
    indx(1) = i
  ELSE IF (ratio > epsi) THEN
    CYCLE
  END IF
  npos = npos + 1
  indx(npos) = i
END DO
GO TO 450

!     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC FOR THE NSTEP = 3 CASE

440 npos = 0
epsi = xmax
DO i = 1, m
  IF (y(i) < zero) THEN
    IF (ibasis(i) <= n) CYCLE
    ip = i
    GO TO 500
  ELSE IF (y(i) > zero) THEN
    ratio = xb(i)/y(i)
    IF (ratio < epsi) THEN
      epsi = ratio
      npos = 1
      indx(1) = i
    ELSE IF (ratio > epsi) THEN
      CYCLE
    END IF
    npos = npos + 1
    indx(npos) = i
  END IF
END DO

450 IF (npos /= 0) GO TO 460
IF (icount >= 5) GO TO 100
ind = 4
GO TO 220

!              TIE BREAKING PROCEDURE

460 ip = indx(1)
IF (npos == 1) GO TO 500
ip = 0
bmin = xmax
cmin = xmax
DO ii = 1, npos
  i = indx(ii)
  l = ibasis(i)
  IF (l > n0) GO TO 461
  IF (c(l) <= zero) cmin = MIN(zero, cmin)
  IF (c(l) > cmin) CYCLE
  imin = i
  cmin = c(l)
  CYCLE
  461 IF (l <= n) GO TO 462
  ip = i
  GO TO 500
  462 lrow = l - n0
  s = b0(lrow)
  IF (lrow <= numle) THEN
    IF (s > bmin) CYCLE
    ip = i
    bmin = s
  ELSE
    s = -s
    bmin = MIN(zero, bmin)
    IF (s > bmin) CYCLE
    ip = i
    bmin = s
  END IF
END DO
IF (cmin <= zero .OR. ip == 0) ip = imin

!               TRANSFORMATION OF XB

500 IF (xb(ip) == zero) GO TO 510
const = xb(ip)/y(ip)
DO i = 1, m
  s = xb(i)
  xb(i) = xb(i) - const*y(i)
  IF (xb(i) >= zero) CYCLE
  IF (s >= zero .OR. xb(i) >= rerrmx*s) xb(i) = zero
END DO
xb(ip) = const

!               TRANSFORMATION OF BI

510 DO j = 1, m
  IF (bi(ip,j) == zero) CYCLE
  const = bi(ip,j)/y(ip)
  bi(1:m,j) = bi(1:m,j) - const*y(1:m)
  bi(ip,j) = const
END DO

!             UPDATING IBASIS AND BASIS

iout = ibasis(ip)
ibasis(ip) = jp
basis(iout) = 0
basis(jp) = 1
IF (iout > n) num = num - 1

!        CHECK THE ACCURACY OF BI AND RESET RERR

IF (rerr > 1.d-2) GO TO 530
k = 0
DO j = 1, m
  kj = ibasis(j)
  IF (kj > n0) CYCLE
  total = DOT_PRODUCT( bi(j,1:m), a(1:m, kj) )
  rerr = MAX(rerr, ABS(1.0_dp - total))
  k = k + 1
  IF (k >= mcheck) GO TO 522
END DO
522 IF (rerr <= 1.d-2) GO TO 600

!        THE ACCURACY CRITERIA ARE NOT SATISFIED

530 IF (icount < 5) GO TO 600
bflag = 1
GO TO 100

580 IF (iter == 0) GO TO 12
IF (bflag == 0) GO TO 590
bflag = 0
DO ip = 1, m
  IF (jp == ibasis(ip)) GO TO 582
END DO
582 ibasis(ip) = iout
basis(jp) = 0
basis(iout) = 1
IF (iout > n) num = num + 1
GO TO 100

590 ind = 3
GO TO 220

!       SET UP THE R ARRAY FOR THE NSTEP = 1 CASE

600 IF (nstep == 2) GO TO 630
IF (nstep > 2) GO TO 700
601 DO j = 1, m
  IF (xb(j) < zero) GO TO 610
END DO
GO TO 630

610 nstep = 1
m0 = 0
DO l = 1, m
  IF (xb(l) >= zero) CYCLE
  m0 = m0 + 1
  indx(m0) = l
END DO

DO j = 1, m
  dsump = zero
  dsumn = zero
  DO ll = 1, m0
    l = indx(ll)
    IF (bi(l,j) < zero) THEN
      dsumn = dsumn + bi(l,j)
    ELSE IF (bi(l,j) > zero) THEN
      dsump = dsump + bi(l,j)
    END IF
  END DO
  y(j) = dsump + dsumn
  s = dsump
  t = dsumn
  tol = rerrmx*MAX(s, -t)
  IF (ABS(y(j)) < tol) y(j) = zero
END DO
GO TO 650

!       SET UP THE R ARRAY FOR THE NSTEP = 2 CASE

630 IF (n == num) GO TO 680
nstep = 2
m0 = 0
DO l = 1, m
  IF (ibasis(l) <= n) CYCLE
  m0 = m0 + 1
  indx(m0) = l
END DO

DO j = 1, m
  dsump = zero
  dsumn = zero
  DO ll = 1, m0
    l = indx(ll)
    IF (bi(l,j) < zero) THEN
      dsumn = dsumn + bi(l,j)
    ELSE IF (bi(l,j) > zero) THEN
      dsump = dsump + bi(l,j)
    END IF
  END DO
  y(j) = -(dsump + dsumn)
  s = dsump
  t = dsumn
  tol = rerrmx*MAX(s, -t)
  IF (ABS(y(j)) < tol) y(j) = zero
END DO

650 DO j = 1, n0
  IF (basis(j) == 0) THEN
    r(j) = DOT_PRODUCT( y(1:m), a(1:m,j) )
  ELSE
    r(j) = zero
  END IF
END DO

660 IF (n0 == ns) GO TO 670
jmin = n0 + 1
DO j = jmin, ns
  r(j) = zero
  IF (basis(j) /= 0) CYCLE
  jj = j - n0
  r(j) = y(jj)
END DO

670 IF (ns == n) GO TO 200
jmin = ns + 1
DO j = jmin, n
  r(j) = zero
  IF (basis(j) /= 0) CYCLE
  jj = j - n0
  r(j) = -y(jj)
END DO
GO TO 200

!      SET UP A NEW R ARRAY FOR THE NSTEP = 3 CASE

680 nstep = 3
DO j = 1, m
  dsum = zero
  DO l = 1, m
    il = ibasis(l)
    IF (il <= n0) dsum = dsum + c(il)*bi(l,j)
  END DO
  y(j) = dsum
END DO

DO j = 1, n0
  r(j) = zero
  IF (basis(j) /= 0) CYCLE
  dsum = -c(j) + DOT_PRODUCT( y(1:m), a(1:m,j) )
  r(j) = dsum
  IF (r(j) >= zero) CYCLE
  tol = rerrmx*ABS(c(j))
  IF (ABS(r(j)) < tol) r(j) = zero
END DO
GO TO 660

!       UPDATE THE R ARRAY FOR THE NSTEP = 3 CASE

700 const = r(jp)
DO j = 1, n0
  IF (basis(j) /= 0) THEN
    r(j) = zero
  ELSE
    total = DOT_PRODUCT( bi(ip,1:m), a(1:m,j) )
    r(j) = r(j) - const*total
    IF (r(j) >= zero) CYCLE
    tol = rerrmx*ABS(c(j))
    IF (ABS(r(j)) < tol) r(j) = zero
  END IF
END DO

IF (n0 == ns) GO TO 720
jmin = n0 + 1
DO j = jmin, ns
  IF (basis(j) /= 0) THEN
    r(j) = zero
  ELSE
    jj = j - n0
    r(j) = r(j) - const*bi(ip,jj)
  END IF
END DO

720 IF (ns == n) GO TO 200
jmin = ns + 1
DO j = jmin, n
  IF (basis(j) /= 0) THEN
    r(j) = zero
  ELSE
    jj = j - n0
    r(j) = r(j) + const*bi(ip,jj)
  END IF
END DO
GO TO 200
!-----------------------------------------------------------------------
!               REFINE XB AND STORE THE RESULT IN Y
!-----------------------------------------------------------------------
800 y(1:m) = zero

m0 = 0
DO j = 1, m
  kj = ibasis(j)
  IF (kj <= n0) GO TO 810
  IF (kj <= ns) GO TO 820
  IF (kj <= n) GO TO 830
  GO TO 820

  810 m0 = m0 + 1
  indx(m0) = j
  CYCLE

  820 l = kj - n0
  y(l) = xb(j)
  CYCLE

  830 l = kj - n0
  y(l) = -xb(j)
END DO

IF (m0 == 0) THEN
  r(1:m) = b0(1:m) - y(1:m)
ELSE
  DO i = 1, m
    dsum = y(i)
    DO jj = 1, m0
      j = indx(jj)
      kj = ibasis(j)
      dsum = dsum + a(i,kj)*xb(j)
    END DO
    r(i) = b0(i) - dsum
  END DO
END IF

rerr1 = MIN(rerrmx, rerr)
DO i = 1, m
  y(i) = zero
  IF (xb(i) < zero) THEN
    sgn = -1.0_dp
    dsump = zero
    dsumn = xb(i)
  ELSE IF (xb(i) > zero) THEN
    sgn = 1.0_dp
    dsump = xb(i)
    dsumn = zero
  ELSE
    CYCLE
  END IF
  DO l = 1, m
    dt = bi(i,l)*r(l)
    IF (dt > zero) THEN
      dsump = dsump + dt
    ELSE
      dsumn = dsumn + dt
    END IF
  END DO
  w = dsump + dsumn
  IF (w == zero) CYCLE
  IF (sgn /= SIGN(1.0_dp, w)) CYCLE
  s = dsump
  t = dsumn
  tol = rerr1*MAX(s, -t)
  IF (ABS(w) > tol) y(i) = w
END DO
IF (nstep == 2) GO TO 870
IF (nstep > 2) GO TO 880

!         CHECK THE REFINEMENT (NSTEP = 1)

DO i = 1, m
  IF (y(i) >= zero) GO TO 861
  IF (y(i) < -rerrmx) GO TO 240
  y(i) = zero
  861 xb(i) = y(i)
END DO
GO TO 630

!         CHECK THE REFINEMENT (NSTEP = 2)

870 DO i = 1, m
  IF (ibasis(i) <= n) GO TO 871
  IF (y(i) > rerrmx) GO TO 240
  y(i) = zero
  871 xb(i) = y(i)
END DO
GO TO 680

!              COMPUTE Z  (NSTEP = 3)

880 dsum = zero
DO i = 1, m
  ki = ibasis(i)
  IF (ki > n0) GO TO 881
  dsum = dsum + c(ki)*y(i)
  881 xb(i) = y(i)
END DO
z = dsum
GO TO 220
END SUBROUTINE smplx1

SUBROUTINE crout1(a, ka, n, iend, indx, temp, ierr)
!     ******************************************************************
!     CROUT PROCEDURE FOR INVERTING MATRICES
!     ******************************************************************
!     A IS A MATRIX OF ORDER N WHERE N IS GREATER THAN OR EQUAL TO 1.
!     THE INVERSE OF A IS COMPUTED AND STORED IN A.

!     KA = LENGTH OF THE COLUMNS OF THE ARRAY A

!     IEND MAY BE 0, 1, ..., N-1.  IT IS ASSUMED THAT EACH OF THE FIRST
!     IEND COLUMNS OF THE MATRIX A CONTAINS ONLY ONE NONZERO ELEMENT
!     AND THAT THE NONZERO ELEMENT IS 1 OR -1.

!     indx IS AN ARRAY OF DIMENSION N-1 OR LARGER THAT IS USED BY THE
!     ROUTINE FOR KEEPING TRACK OF THE ROW INTERCHANGES THAT ARE MADE.

!     TEMP IS A TEMPORARY STORAGE ARRAY.

!     IERR REPORTS THE STATUS OF THE RESULTS. IF A IS NONSINGULAR THEN
!     THE INVERSE OF A IS COMPUTED AND IERR=0. OTHERWISE IF A IS FOUND
!     TO BE SINGULAR THEN IERR=1 AND THE ROUTINE TERMINATES.
!     --------------------
USE constants_NSWC
IMPLICIT NONE
INTEGER, INTENT(IN)                     :: ka, n, iend
INTEGER, INTENT(OUT)                    :: ierr
REAL (dp), DIMENSION(:), INTENT(IN OUT) :: a, temp
INTEGER, DIMENSION(:), INTENT(IN OUT)   :: indx

!     Local variables
INTEGER   :: i, ibeg, ij, ik, il, j, jcol, jj, k, kcol, kcount, kj,  &
             kj0, kk, kl, km1, kp1, l, lj, lj0, lk, lmin, maxdim, mcol,  &
             ncol, nk, nm1, nmj, nmk, nn
REAL (dp) :: zero = 0._dp
REAL (dp) :: dsum, c, pmin, s

maxdim = ka*n
mcol = iend*ka
IF (iend == 0) GO TO 100

!           PROCESS THE FIRST IEND COLUMNS OF A

kcol = 0
DO k = 1, iend
  kk = kcol + k
  nk = kcol + n
  DO lk = kk, nk
    IF (a(lk) < zero) GO TO 20
    IF (a(lk) > zero) GO TO 30
  END DO
  GO TO 300

  20 l = lk - kcol
  lj0 = mcol + l
  DO lj = lj0, maxdim, ka
    a(lj) = -a(lj)
  END DO

  30 l = lk - kcol
  indx(k) = l
  IF (k == l) GO TO 32
  lj = lk
  DO kj = kk, maxdim, ka
    c = a(kj)
    a(kj) = a(lj)
    a(lj) = c
    lj = lj + ka
  END DO
  32 kcol = kcol + ka
END DO

!           PROCESS THE REMAINING COLUMNS OF A

100 nm1 = n - 1
ierr = 0
pmin = zero
ibeg = iend + 1
IF (ibeg == n) GO TO 190

k = ibeg
km1 = iend
kp1 = k + 1
kcol = mcol
kk = kcol + k
DO kcount = ibeg, nm1

!     SEARCH FOR THE K-TH PIVOT ELEMENT (K=IBEG, ..., N-1)

  l = k
  s = ABS(a(kk))
  DO i = kp1, n
    ik = kcol + i
    c = ABS(a(ik))
    IF (s >= c) CYCLE
    l = i
    s = c
  END DO

  IF (k > ibeg .AND. s >= pmin) GO TO 120
  pmin = s
  IF (s == zero) GO TO 300

!              INTERCHANGING ROWS K AND L

  120 indx(k) = l
  IF (k == l) GO TO 130
  kj0 = mcol + k
  lj  = mcol + l
  DO kj = kj0, maxdim, ka
    c = a(kj)
    a(kj) = a(lj)
    a(lj) = c
    lj = lj + ka
  END DO

!       COMPUTE THE K-TH ROW OF U (K=IBEG, ..., N-1)

  130 c = a(kk)
  IF (k > ibeg) GO TO 140
  kj0 = kk + ka
  DO kj = kj0, maxdim, ka
    a(kj) = a(kj)/c
  END DO
  GO TO 160

  140 kl = mcol + k
  DO l = ibeg, km1
    temp(l) = a(kl)
    kl = kl + ka
  END DO

  kj0 = kk + ka
  DO kj = kj0, maxdim, ka
    jcol = kj - k
    dsum = -a(kj)
    DO l = ibeg, km1
      lj = jcol + l
      dsum = dsum + temp(l)*a(lj)
    END DO
    a(kj) = -dsum/c
  END DO

!      COMPUTE THE K-TH COLUMN OF L (K=IBEG+1, ..., N)

  160 km1 = k
  k = kp1
  kp1 = k + 1
  kcol = kcol + ka
  kk = kcol + k
  DO l = ibeg, km1
    lk = kcol + l
    temp(l) = a(lk)
  END DO

  DO i = k, n
    il = mcol + i
    dsum = zero
    DO l = ibeg, km1
      dsum = dsum + a(il)*temp(l)
      il = il + ka
    END DO
    a(il) = a(il) - dsum
  END DO
END DO

!           CHECK THE N-TH PIVOT ELEMENT

190 ncol = maxdim - ka
nn = ncol + n
c = ABS(a(nn))
IF (c > pmin) GO TO 200
IF (c == zero) GO TO 300

!          REPLACE L WITH THE INVERSE OF L

200 IF (ibeg == n) GO TO 213
jj = mcol + ibeg
i = ka + 1
DO j = ibeg, nm1
  a(jj) = 1.0_dp/a(jj)
  temp(j) = a(jj)
  kj = jj
  DO km1 = j, nm1
    k = km1 + 1
    kj = kj + 1
    dsum = zero
    kl = kj
    DO l = j, km1
      dsum = dsum + a(kl)*temp(l)
      kl = kl + ka
    END DO
    a(kj) = -dsum/a(kl)
    temp(k) = a(kj)
  END DO
  jj = jj + i
END DO
213 a(nn) = 1.0_dp/a(nn)
IF (n == 1) RETURN

!       SOLVE UX = Y WHERE Y IS THE INVERSE OF L

DO nmk = 1, nm1
  k = n - nmk
  lmin = MAX(ibeg, k+1)
  kl = (lmin-1)*ka + k
  DO l = lmin, n
    temp(l) = a(kl)
    a(kl) = zero
    kl = kl + ka
  END DO

  kj0 = mcol + k
  DO kj = kj0, maxdim, ka
    dsum = -a(kj)
    lj = (kj - k) + lmin
    DO l = lmin, n
      dsum = dsum + temp(l)*a(lj)
      lj = lj + 1
    END DO
    a(kj) = -dsum
  END DO
END DO

!                 COLUMN INTERCHANGES

jcol = ncol - ka
DO nmj = 1, nm1
  j = n - nmj
  k = indx(j)
  IF (j == k) GO TO 251
  ij = jcol
  ik = (k-1)*ka
  DO i = 1, n
    ij = ij + 1
    ik = ik + 1
    c = a(ij)
    a(ij) = a(ik)
    a(ik) = c
  END DO
  251 jcol = jcol - ka
END DO
RETURN

!                    ERROR RETURN

300 ierr = 1
RETURN
END SUBROUTINE crout1

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
      program convertcoords4games
      implicit double precision(a-h,o-z)

      character*2 lab
      real*8 c(3)
      real an
      integer alab2anum

c     open(unit=1,file='gamesscoords',status='unknown')

1     continue
      read(5,*,end=2)lab,(c(k),k=1,3)
      n=alab2anum(lab)
      an=float(n)
c     write(1,100)lab,an,(c(k),k=1,3)
      write(6,100)lab,an,(c(k),k=1,3)
100   format(a2,f6.1,3f13.6)
      go to 1
2     continue
c     close(unit=1)
      end

      integer function alab2anum(alab)

c given the atom label this returns the atomic number

      character*2 :: alab,atomic_labels(110)

      data atomic_labels/'H','He','Li','Be','B','C','N','O','F','Ne',
     .'Na','Mg','Al','Si','P','S','Cl','Ar','K','Ca','Sc','Ti','V','Cr',
     . 'Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb',
     . 'Sr','Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn',
     . 'Sb','Te','I','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu',
     . 'Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','W','Re','Os',
     . 'Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac',
     .   'Th','Pa','U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No'
     . ,'Lr','Rf','Db','Sg','Bh','Hs','Mt','Ds'/

      integer n

      do n=1,110
        if(atomic_labels(n).eq.alab)then
      alab2anum=n
      return
        endif
      enddo
      alab2anum=0

      return
      end

max_stars_repo_name: mickcollins/SMFA

max_issues_repo_name: mickcollins/SMFA

max_forks_repo_name: mickcollins/SMFAPAC


text:
About Melissa

Melissa loves giving massages to help people. She believes that it challenges the body, mind, and spirit while staying centered and focused on the needs of the person or animal at hand. She states It is amazing to feel a muscle change shape and soften under my fingertips, and exciting for a client to express how great they feel after a session. Melissas mission for her massage practice is to help families live a pain and stress free life physically, mentally, and spiritually. One of her deepest passions is helping families during childbearing years. Her vision for teaching infant massage to parents and caregivers is to help babies become calmer, smarter, sweeter individuals, who grow up to see the value of massage in their lives.  She hopes to set a positive example of living a balanced life by taking care of her health, her family and friends, and everyone who comes into her practice.

Melissas focus is providing a comforting cozy space so that when someone walks into the studio, they feel they are in a healing sanctuary. She communicates with each individual upfront about what they want her to focus on, and strives to be centered and intuitive enough to feel the muscle or scar tissue that needs work. She tries to sense the perfect amount of pressure for allowing each muscle fiber to melt back into a normal perfect physical state, and bring blood flow back into the area. Melissa strongly encourages any feedback during a session. Melissas goal is for her clients to come in and allow her to help reduce their stress and relieve their muscle pain.

Training

Melissa became certified in massage therapy in 2004 and completed CPR training in Davis. She is licensed for massage in the city of Davis and has a massage practitioners license through the state of California (lic.# 18681). She is insured through ABMP.  In July of 2005 she traveled to The Boulder College of Massage and was trained in canine massage therapy.  In October of 2005 she was trained in equine massage therapy through Healing Arts Institute.  Melissa finished a 35 hour certification program called Pregnancy Massage for the Childbearing year from Stephanie Halderman at Healing Arts Institute in September 2010, and belongs to the American Pregnancy Massage Association. In 2012 she became a Certified Infant Massage Instructor CIMI® through the International Loving Touch of Foundation, a national accredited program. I am continually taking continued education courses to help my clients.

Background

Melissa has been a veterinary technician since 1996. In 2002 she was introduced to alternative medicine, and saw the speed and ease of recovery from surgery. Her dog, Sasha, had her anterior cruciate ligament repaired. Pain medications, physical therapy, acupuncture, and massage were used during her recovery. Melissa saw decreased drug use, decreased recovery time, and stress reduction. She realized through her training and also receiving massage on a monthly basis that routine massage is part of the physical and mental equation of being healthy. Massage brings balance to the body, along with eating healthy, exercising, drinking enough water, and getting enough rest.

Website Testimonals

I would highly recommend Melissa as a massage therapist. I started running last year and have noticed a huge difference in recovery time after getting a Sports Massage from her. She listens to you and takes the time to understand your specific issue and/or needs. I will not go to anyone else because she is the best. I met her at the 2009 Turkey Trot and have been hooked since then. She is caring, patient and gives the absolute BEST massages. Try her and you will agree! — Arlene

I highly recommend a massage with Melissa. She was fantastic at listening and adjusting to what I wanted. The atmosphere was great as well (super soft sheets and a fountain). Ive tried a few places around Davis and this was definitely my favorite  and actually one of the best massages Ive had anywhere. — JLW

I just got a massage from Melissa and it was wonderful! I was a little unsure after not seeing any comments here, but the ability to book the appointment online was what encouraged me to check this place out. Melissa asked good questions to see what kind of massage I was looking for and made me very comfortable. She was professional and her technique was amazing. I dont think Ive had any better massage before. Highly recommended! — Adrienne W.

20140416 22:54:50 nbsp  I started getting pregnancy massages from Melissa in my third trimester after I started having significant sciatica pain. After only two sessions, the pain was absolutely gone. That was all it took to convince me that she was going to be key to helping me get through the pain and discomfort of my third trimester. Ive ended up booking weekly appointments and it has helped me tremendously with lower back pain, swelling, and plantar fasciitis. She always tailors my massages to focus on current trouble spots which is really beneficial. If you are pregnant and feeling discomfort, book an appointment now!  Users/darcid4

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    Copyright(C) 2008-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
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
C=======================================================================
      SUBROUTINE STORE (ISTEP, TYP, IBEGIN, IEND, NWRDS,
     &   NUMELB, IDELB, ISEVOK, VISELB, MAXNE, VARVAL, MERR)
C=======================================================================

C   --*** STORE *** (ALGEBRA) Store database variable data
C   --   Written by Amy Gilkey - revised 08/16/88
C   --
C   --STORE reads each database variable and stores it, if it is wanted,
C   --into the proper location of VARVAL.  The /VAR../ variables between
C   --IBEGIN and IEND are assumed to be ordered on IDVAR, with no
C   --repetitions.
C   --
C   --Parameters:
C   --   ISTEP  - IN  - the time step number
C   --   TYP    - IN  - the variable types
C   --   IBEGIN - IN  - the starting /VAR../ index of the variables
C   --   IEND   - IN  - the ending /VAR../ index of the variables
C   --   NWRDS  - IN  - the number of words to be read for each variable
C   --   NUMELB - IN  - the number of elements per block
C   --   ISEVOK - IN  - the element block variable truth table;
C   --                  variable i of block j exists iff ISEVOK(j,i)
C   --   MAXNE  - IN  - the VARVAL dimension (max of NUMEL and NUMNP)
C   --   VARVAL - OUT - the returned needed input variables
C   --   MERR   - OUT - error flag C   --
C   --Common Variables:
C   --   Uses IDVAR, ISTVAR of /VAR../
C   --   Uses NDBIN of /DBASE/

      PARAMETER (ICURTM = 1, ILSTTM = 2, IONETM = 3)
      include 'ag_namlen.blk'
      include 'ag_var.blk'
      include 'ag_dbase.blk'
      include 'ag_dbnums.blk'

      CHARACTER TYP
      INTEGER NUMELB(*)
      INTEGER IDELB(*)
      LOGICAL ISEVOK(NELBLK,NVAREL)
      LOGICAL VISELB(NELBLK)
      REAL VARVAL(MAXNE,*)
      INTEGER MERR
      MERR = 0

C     Read and store all global variables (if needed for the current step)
      IF (TYP .EQ. 'G') THEN
         IF (ISTVAR(ICURTM,IBEGIN) .NE. 0) THEN
            NSTO = ISTVAR(ICURTM,IBEGIN)
            ID = IDVAR(IBEGIN)
            call exggv(ndbin, istep, nvargl, varval(id,nsto), ierr)
         END IF
C     Read and store all needed nodal variables for the current step
      ELSE IF (TYP .EQ. 'N') THEN
         DO 100 NVAR = IBEGIN, IEND
            IF (ISTVAR(ICURTM,NVAR) .NE. 0) THEN
               NSTO = IABS (ISTVAR(ICURTM,NVAR))
               call exgnv(ndbin, istep, idvar(nvar), nwrds,
     &                    varval(1,nsto), ierr)
            END IF
  100    CONTINUE
C     Read and store all needed element variables for the current step
      ELSE IF (TYP .EQ. 'E') THEN
         IEL = 1
         DO 130 IELB = 1, NELBLK
            IDEB = IDELB(IELB)
            NELEM = NUMELB(IELB)
            IF (VISELB(IELB)) THEN
              DO 120 NVAR = IBEGIN, IEND
               INVAR = IDVAR(NVAR)
               IF (ISTVAR(ICURTM,NVAR) .NE. 0) THEN
                  NSTO = IABS (ISTVAR(ICURTM,NVAR))
                  IF (ISEVOK(IELB,INVAR)) THEN
                    call exgev(ndbin, istep, INVAR,
     &                IDEB, NELEM, varval(iel,nsto), ierr)
                  ELSE
C                  --Make sure values for undefined elements are zero
                     DO 110 I = IEL, IEL+NUMELB(IELB) - 1
                        VARVAL(I,NSTO) = 0.0
  110                CONTINUE
                  END IF
               END IF
  120         CONTINUE
            END IF
            IEL = IEL + NELEM
  130    CONTINUE
      ELSE
         WRITE(*,10)
     &   'Unknown Type in Subroutine Store: TYP =  ',typ
 10      format(A,A)
         MERR = 1
      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


module iteration
  ! Module to perform the field iteration
  use globdata
  use io, only: print_array
  implicit none
  include "mpif.h"

  integer :: info

  private
  public :: iterate, par_iterate, stitch

contains

  subroutine iterate
    real(mkt) :: fac1, fac2, fac3
    integer :: i, j, n

    fac1 = d*dt
    fac2 = 1_mkd/dx**2; fac3 = 1_mkd/dy**2

    do n = n_start, n_max
      if ( mod(n,2) == 0 ) then
        do j = 2, nx_loc-1
          do i = 2, ny_loc-1
            T(i,j) = fac1*( (T_old(i+1,j) - 2*T_old(i,j) + T_old(i-1,j))*fac2 + &
              (T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1))*fac3 ) + T_old(i,j)
          enddo
        enddo
      else
        do j = 2, nx_loc-1
          do i = 2, ny_loc-1
            T_old(i,j) = fac1*( (T(i+1,j) - 2*T(i,j) + T(i-1,j))*fac2 + &
              (T(i,j+1) - 2*T(i,j) + T(i,j-1))*fac3 ) + T(i,j)
          enddo
        enddo
      endif
    enddo

    ! In case n_max is not even update the T field one last time
    if ( mod(n_max,2) /= 0 ) then
      do j = 2, nx-1
        do i = 2, ny-1
          T(i,j) = fac1*( (T_old(i+1,j) - 2*T_old(i,j) + T_old(i-1,j))*fac2 + &
            (T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1))*fac3 ) + T_old(i,j)
        enddo
      enddo
    endif

  end subroutine iterate

  subroutine par_iterate
    real(mkt) :: fac1, fac2, fac3
    real(mkt), dimension(1:ny_loc-2) :: send_buffer, recv_buffer
    integer :: i, j, n, m, ndest, nsrc, send_column, recv_column
    logical :: even

    fac1 = d*dt
    fac2 = 1_mkd/dx**2; fac3 = 1_mkd/dy**2

    ! Iterate
    do n = n_start, n_max
      if ( mod(n,2) == 0 ) then; even = .true.; else; even = .false.; endif

        if ( even ) then
          do j = 2, nx_loc-1
            do i = 2, ny_loc-1
              T(i,j) = fac1*( (T_old(i+1,j) - 2*T_old(i,j) + T_old(i-1,j))*fac2 + &
                (T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1))*fac3 ) + T_old(i,j)
            enddo
          enddo
        else
          do j = 2, nx_loc-1
            do i = 2, ny_loc-1
              T_old(i,j) = fac1*( (T(i+1,j) - 2*T(i,j) + T(i-1,j))*fac2 + &
                (T(i,j+1) - 2*T(i,j) + T(i,j-1))*fac3 ) + T(i,j)
            enddo
          enddo
        endif

        ! Exhange ghost layers, first right then left
        exchange: do m = 1, 2
          ! Define exchange direction
          if ( m == 1 ) then
            send_column = nx_loc-1
            recv_column = 1
            ndest = rank+1
            nsrc  = rank-1
            if (ndest == pnum) ndest = MPI_Proc_null ! Check if we are at end
            if (nsrc == -1)    nsrc  = MPI_Proc_null ! columns and if so nullify
          else
            send_column = 2
            recv_column = nx_loc
            ndest = rank-1
            nsrc  = rank+1
            if (ndest == -1)   ndest = MPI_Proc_null
            if (nsrc == pnum)  nsrc  = MPI_Proc_null
          endif

          ! Define buffers to be sent
          if ( even ) then
            do i = 1, ny_loc-2
              send_buffer(i) = T(i+1,send_column)
            enddo
          else
            do i = 1, ny_loc-2
              send_buffer(i) = T_old(i+1,send_column)
            enddo
          endif

          ! Send/recv ghost layers
          call MPI_Sendrecv(send_buffer, ny_loc-2, MPI_Double_precision, &
            ndest, 0, recv_buffer, ny_loc-2, MPI_Double_precision, nsrc, 0, &
            MPI_Comm_world, MPI_Status_ignore, info)
          if ( info /= 0 ) print *, 'Sendrecv error'

          if ( nsrc == -2 ) cycle exchange ! Cycle if source is null

          ! Update ghost layers
          if ( even ) then
            do i = 1, ny_loc-2
              T(i+1,recv_column) = recv_buffer(i)
            enddo
          else
            do i = 1, ny_loc-2
              T_old(i+1,recv_column) = recv_buffer(i)
            enddo
          endif
        enddo exchange
      enddo

      ! In case n_max is not even update the T field one last time
      if ( mod(n_max,2) /= 0 ) then
        do j = 2, nx-1
          do i = 2, ny-1
            T(i,j) = fac1*( (T_old(i+1,j) - 2*T_old(i,j) + T_old(i-1,j))*fac2 + &
              (T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1))*fac3 ) + T_old(i,j)
          enddo
        enddo
      endif
    end subroutine par_iterate

    subroutine stitch
      use alloc_mod

      real(8), dimension(:,:), allocatable :: rbuf_T
      integer :: info, irank, cur_nx, cur_ny, i, j, col_count

      if ( rank > 0 ) then
        call MPI_Send(T, size(T,1)*size(T,2), MPI_Double_precision, 0, 0, &
          MPI_Comm_world, info)
      endif

      if ( rank == 0 ) then
        ! Expand the root field array to global size (nx x ny)
        call alloc_field(T, nx, ny, info)

        ! Increment counter denoting the last inserted column (-1 removes g. layer)
        col_count = field_sizes(2)-1

        ! Insert remaining fields
        do irank = 1, pnum-1
          ! Get size of the current array (irank+1 to skip root's index)
          cur_ny = field_sizes((irank+1)*2-1)
          cur_nx = field_sizes((irank+1)*2)

          if ( irank == 1 ) then
            allocate(rbuf_T(cur_ny, cur_nx))
          elseif (cur_ny /= field_sizes(irank*2-1) .or. &
              cur_nx /= field_sizes(irank*2  )) then
              deallocate(rbuf_T)
              allocate(rbuf_T(cur_ny, cur_nx))
            endif

            call MPI_Recv(rbuf_T, cur_ny*cur_nx, MPI_Double_precision, irank, 0,&
              MPI_Comm_world, MPI_Status_ignore, info)

            do j = 2, cur_nx
              do i = 1, cur_ny
                T(i, col_count+j-1) = rbuf_T(i,j)
              enddo
            enddo

            ! Update last column (-2 since internal slices have two ghost layers)
            col_count = col_count+cur_nx-2

          enddo
          deallocate(rbuf_T)
        endif
      end subroutine stitch
    end module iteration


text:
*DECK ORTRAN
      SUBROUTINE ORTRAN (NM, N, LOW, IGH, A, ORT, Z)
C***BEGIN PROLOGUE  ORTRAN
C***PURPOSE  Accumulate orthogonal similarity transformations in the
C            reduction of real general matrix by ORTHES.
C***LIBRARY   SLATEC (EISPACK)
C***CATEGORY  D4C4
C***TYPE      SINGLE PRECISION (ORTRAN-S)
C***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
C***AUTHOR  Smith, B. T., et al.
C***DESCRIPTION
C
C     This subroutine is a translation of the ALGOL procedure ORTRANS,
C     NUM. MATH. 16, 181-204(1970) by Peters and Wilkinson.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C
C     This subroutine accumulates the orthogonal similarity
C     transformations used in the reduction of a REAL GENERAL
C     matrix to upper Hessenberg form by  ORTHES.
C
C     On INPUT
C
C        NM must be set to the row dimension of the two-dimensional
C          array parameters, A and Z, as declared in the calling
C          program dimension statement.  NM is an INTEGER variable.
C
C        N is the order of the matrix A.  N is an INTEGER variable.
C          N must be less than or equal to NM.
C
C        LOW and IGH are two INTEGER variables determined by the
C          balancing subroutine  BALANC.  If  BALANC  has not been
C          used, set LOW=1 and IGH equal to the order of the matrix, N.
C
C        A contains some information about the orthogonal trans-
C          formations used in the reduction to Hessenberg form by
C          ORTHES  in its strict lower triangle.  A is a two-dimensional
C          REAL array, dimensioned A(NM,IGH).
C
C        ORT contains further information about the orthogonal trans-
C          formations used in the reduction by  ORTHES.  Only elements
C          LOW through IGH are used.  ORT is a one-dimensional REAL
C          array, dimensioned ORT(IGH).
C
C     On OUTPUT
C
C        Z contains the transformation matrix produced in the reduction
C          by  ORTHES  to the upper Hessenberg form.  Z is a two-
C          dimensional REAL array, dimensioned Z(NM,N).
C
C        ORT has been used for temporary storage as is not restored.
C
C     Questions and comments should be directed to B. S. Garbow,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     ------------------------------------------------------------------
C
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
C                 system Routines - EISPACK Guide, Springer-Verlag,
C                 1976.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   760101  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  ORTRAN
C
      INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1
      REAL A(NM,*),ORT(*),Z(NM,*)
      REAL G
C
C     .......... INITIALIZE Z TO IDENTITY MATRIX ..........
C***FIRST EXECUTABLE STATEMENT  ORTRAN
      DO 80 I = 1, N
C
         DO 60 J = 1, N
   60    Z(I,J) = 0.0E0
C
         Z(I,I) = 1.0E0
   80 CONTINUE
C
      KL = IGH - LOW - 1
      IF (KL .LT. 1) GO TO 200
C     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
      DO 140 MM = 1, KL
         MP = IGH - MM
         IF (A(MP,MP-1) .EQ. 0.0E0) GO TO 140
         MP1 = MP + 1
C
         DO 100 I = MP1, IGH
  100    ORT(I) = A(I,MP-1)
C
         DO 130 J = MP, IGH
            G = 0.0E0
C
            DO 110 I = MP, IGH
  110       G = G + ORT(I) * Z(I,J)
C     .......... DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES.
C                DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ..........
            G = (G / ORT(MP)) / A(MP,MP-1)
C
            DO 120 I = MP, IGH
  120       Z(I,J) = Z(I,J) + G * ORT(I)
C
  130    CONTINUE
C
  140 CONTINUE
C
  200 RETURN
      END


text:
      PROGRAM Surfvmd
      IMPLICIT NONE
C ********************************************************************
C *                                                                  *
C * This software is an unpublished work containing confidential and *
C * proprietary information of Birkbeck College. Use, disclosure,    *
C * reproduction and transfer of this work without the express       *
C * written consent of Birkbeck College are prohibited. This notice  *
C * must be attached to all copies or extracts of the software.      *
C *                                                                  *
C * (c) 2002 Oliver Smart & Birmingham University                    *
C *                                                                  *
C ********************************************************************
C
C Modification history:
C
C Date  Author          Modification
C 05/02 O.S.S. 		Original version
C
C reads in the

C routine FREDA used here to decode lines include common vbles
C vbles used by freda
C col: 	c*80		the command?
C kl: 	int 		no of words inputed 0 then default
C fn:	r*4(40) 	no.'s as they occur
C kn:     int     	no. of no's
C fl:     c*1(40,6)	words?
      INCLUDE 'FREDAINC'

C
C screen, keyboard
C sph_process now command line driven - comment out all
C interactive question lines with !(inter)
!(cinter)      INTEGER                   NIN
!(cinter)      PARAMETER(                NIN = 5)
!(cinter)C input/output file stream
!(cinter)      INTEGER                   SIN, SOUT
      INTEGER                   NOUT
      PARAMETER(                NOUT= 6)

      CHARACTER*200		FINPUT      ! inputfilename (penultimate argument)
      INTEGER			SIN         ! its stream number
      INTEGER			END_FINPUT  ! last character in filename
      CHARACTER*200		FOUTPUT     ! inputfilename (penultimate argument)
      INTEGER			SOUT        ! its stream number
      INTEGER			END_FOUTPUT ! last character in filename
      CHARACTER*200 		FSPHINPUT
      INTEGER			SIN2         ! its stream number
      INTEGER			END_FSPHINPUT

!(cinter)C filename
!(cinter)      CHARACTER*200             FNAME

C string for text records in hydra file
      CHARACTER*80              STRING

C abort indicator
      LOGICAL                   LABORT

C draw centre lines?
      LOGICAL                   LCEN

C override radius?
      DOUBLE PRECISION		ROVER

C lower limit radius - any smaller than this and ignore on read
      DOUBLE PRECISION		RLOWER
C number of ignored records
      INTEGER			NIGN

C sphere centre storage arrays
C maximum number, actual number, co-ords and radius
C logical indicates whether the record is the last centre
      INTEGER                   SPMAX, SPNO
      PARAMETER(                SPMAX = 500000)
      DOUBLE PRECISION          SPXYZ( 3, SPMAX), SPRAD( SPMAX)
      LOGICAL                   SPLAST( SPMAX)
C capsule option needs to have the second centre
C n.b. SPRAD is the real capsule radius not the effective rad`
      DOUBLE PRECISION          SPSEC( 3, SPMAX)
C store whether point is product of connolly - pick up by -999 	residue number
      LOGICAL			SPCONN( SPMAX)
      ! need effective radius for colouring connolly
      DOUBLE PRECISION		SPEFFR( SPMAX)

! VBLES FOR s/r SPHPEG
! flag to indicate that the sphere is peg accessible
!  0: not peg accessible
!  1: accessible from +ve end
! -1: accessible from -ve end
!  2: in the shell (within distance shell_add of a peg sphere)
! -2: ditto
      INTEGER			SP_PEG_ACC( SPMAX)

! now with  SHELL_ADD need to modify the radius of some spheres
! normally this would be equal to sprad but if sp_peg_acc = +/-2
! then this radius will be lower
      DOUBLE PRECISION		SP_PEG_RAD( SPMAX)

! SCOTDO -> SP_TDO
! SCOENC -> SP_ENC
! SCOCIRCRAD -> SP_CIRCRAD
! n.b. vble may now be changed in calculation if sphere has no intersection
! with plane under consideration
      LOGICAL			SP_TDO( SPMAX)
C s/r coarea needs a store for the east/north coords of
C each point on grid - these are measured relative to the
C initial point
! SCOCIRCRAD: the radius of the point in the plane under consideration
! in the case of calls from concal this we equal scocrad
      DOUBLE PRECISION		SP_ENC( 2, SPMAX),
     &				SP_CIRCRAD( SPMAX)

! working flag for s/r coarea
      LOGICAL			SP_PEG_FLAG( SPMAX)

C loop counts x,y&z; sphere centre count, dots count,
C 2nd sphere centre count, colour Pass count
      INTEGER                   XCOUNT, SCOUNT

C dot density to be used in output
      INTEGER                   DOTDEN

C program uses a sphere of isotropically distributed points
C on a sphere of unit radius. PTMAX is the array bound,
C PTNO is the number of dots (dependent on dotden)
      INTEGER                   PTMAX, PTNO
      PARAMETER(                PTMAX = 10000)
      DOUBLE PRECISION          PTXYZ( 3, PTMAX)

C error indicator
      LOGICAL                   LERR

C variable to define colours
      DOUBLE PRECISION          RCUT(0:3)
C if we are doing colour NPASS=3
C if we are not doing colour NPASS=1
      INTEGER                   NPASS

C if doing a grid surface then check each new potential vertex
C to see if it lies with distance 0.6*rad*dmult (where rad is
C the pore radius at the point) of any line already output
C- stored in linxyz.
C If the new point is outside this range then it is accepted
C it is joined to any vertex within the range rad*dmult the
C straight line is then corrected so that it curves with the local
C pore surface
C
C within a distance DOTD2 (the square of the minimum distance apart).
C VERMAX is the array bound, VERNO is the number of accepted dots.
C fourth number is the radius of the sphere which produced the dot
      INTEGER                   VERMAX, VERNO
      PARAMETER(                VERMAX = 10000)
      DOUBLE PRECISION          VERXYZ( 4, VERMAX)

C also store points taken from output lines
      INTEGER                   LINMAX, LINNO
      PARAMETER(                LINMAX = 30*VERMAX)
      DOUBLE PRECISION          LINXYZ( 4,LINMAX)

C 21/07/94 program changed so that dot is output if it is
C at least dmult*sphere-radius from any other dot
      DOUBLE PRECISION          DMULT

C control flag - if true do a grid surface
      LOGICAL			LGRID

C logical vble capsules on
      LOGICAL			LCAPS

C string for output
      CHARACTER*90		OUTSTR

C doing sos output
      LOGICAL			SOSOUT

      ! vbles for get_rec_commands
      CHARACTER*3000 		REC_COMMANDS_HELP_MESSAGE
      INTEGER			REC_MAX  ! array bound
      PARAMETER(		REC_MAX = 15)
      INTEGER     		REC_TOT  ! actual number of arguements
      CHARACTER*25		REC_COMMANDS(REC_MAX)
      INTEGER			REC_REQUIRE_NUM(REC_MAX)
      INTEGER			REC_REQUIRE_STRING(REC_MAX)
      INTEGER			RCOUNT
      LOGICAL			REC_FOUND( REC_MAX)
      DOUBLE PRECISION		REC_NUM( 10, REC_MAX)
      INTEGER			REC_NUM_NUM( REC_MAX)
      CHARACTER*200		REC_STRING( 10, REC_MAX)
      INTEGER			REC_STRING_NUM( REC_MAX)

      LOGICAL			LBIN		! binary file indicator
      LOGICAL                   OPENRO	! open files readonly function

! vbles for forpeg option - hole channel vector etc.
      DOUBLE PRECISION		CVECT(3)
! unit vectors in plane normal to cvect
! needed to a number of routines - originally h2dmap but now conn
! and addend
! these are normal - first is east direction (x on plane)
!                    and second is north     (y on plane)
      DOUBLE PRECISION		PERPVE(3), PERPVN(3)
! a sampling distance - the distance between planes
      DOUBLE PRECISION		SAMPLE
! option CONN introduced 1/00 -
! CONNR(1) is the connolly surface probe radius (default 1.15Angs)
!          the option is turned off by setting this radius negative
! CONNR(2) is the grid size to be used in the calculation - default 1/2 probe
      DOUBLE PRECISION		CONNR(2)
! the radius at which an end is reached
      DOUBLE PRECISION		ENDRAD

! PEG is on OR NOT
      LOGICAL 			LPEG

! vbles for working peg effect on conductance -
! first number the ratio of conductance with PEG about to water
! second number is the ratio for non-penetrating PEG
      DOUBLE PRECISION		PEGRAT(2)

! probe radius - +ve direction and -ve
      DOUBLE PRECISION		PROBE_PLUS, PROBE_NEG

! (n037) additional shell around peg spheres where the
! conductivity of the solvent is reduced.
      DOUBLE PRECISION		SHELL_ADD

      INTEGER			REC_SPHERE
      DOUBLE PRECISION		TRIA(3),TRIB(3),TRIC(3),
     &				NORMA(3),NORMB(3),NORMC(3)

      INTEGER			TCOUNT
      DOUBLE PRECISION		DIST

! end of decs ******************

!(cinter)C turn on VT codes  with bold after prompt-
!(cinter)      CALL VTCON( .TRUE.)

C initialise the no.of accepted dots to 0
      VERNO = 0

C capsule option defaults to off
      LCAPS = .FALSE.
C initialize
      LERR = .FALSE.

C write greeting
      WRITE( NOUT, '(A)')
     &' *** Program sph_process (replacing sphqpt) ***',
     &' Reads a .sph file produced by hole',
     &' option SPHPDB which contains sphere centre info,',
     &' and processes it to a either dot surface or a number'//
     &                                     ' of other options',
     &' either as a dot or grid representation',
     &' Copyright 1996,1997 by Oliver Smart and Birkbeck College',
     &' Copyright 1998 by Oliver Smart and Birmingham University',
     &' Program modification number 0001 15-11-00'

C write link time of program to screen
      CALL VERTIM( NOUT)

      ! get the controlling command line arguments

      ! establish values for s/r get
      REC_COMMANDS_HELP_MESSAGE =
     &'                     (help message start)\n'//
     &' Usage: '//
     &' surf2vmd  ((options)) '//
     &        ' -sph input_file.sph input_file.tri output_file\n'//
     &' options: \n'//
     &'  -colour (Rlow Rmid): produce a coloured surface\n'//
     &'                       surface will be \n'//
     &'                       red for pore radii < Rlow \n'//
     &'                       green for Rlow < pore radii < Rmid\n'//
     &'                       blue  for pore radii > Rmid\n'//
     &'                       default Rlow = 1.15 angs\n'//
     &'                       default Rmid = 2.30 angs\n'//
     &'                       (for connolly use effective radius\n'//
     &'                        of slab rather than pore radius\n'//
     &'  -color same as -colour for yanks\n'

!      call chrend( rec_commands_help_message, line_end)
!      write( nout, '(a)') rec_commands_help_message(1:line_end)
      REC_TOT = 3
      DO RCOUNT = 1,REC_TOT          ! initialize require store
        REC_REQUIRE_NUM(RCOUNT) = 0
	REC_REQUIRE_STRING(RCOUNT) = 0
      ENDDO
      REC_COMMANDS(1) = 'colour'
      REC_COMMANDS(2) = 'color'
      REC_COMMANDS(3) = 'sph'
      REC_REQUIRE_STRING(3) = 1

      CALL GET_REC_COMMANDS( NOUT, LERR,
     &			     REC_COMMANDS_HELP_MESSAGE,
     &                       REC_TOT, REC_COMMANDS,
     &                       REC_REQUIRE_NUM, REC_REQUIRE_STRING,
     &                       REC_FOUND, REC_NUM, REC_NUM_NUM,
     &			     REC_STRING, REC_STRING_NUM,
     &			     .TRUE., FINPUT, .TRUE., FOUTPUT)
      IF (LERR) GOTO 55555

      ! REC_COMMANDS(1) = 'colour'
      ! REC_COMMANDS(2) = 'color'
      ! if we are doing colour NPASS=3
      ! if we are not doing colour NPASS=1 accept any radii
      NPASS = 1
      RCUT(0) = -1.0
      RCUT(1) = 9999.0
      IF (REC_FOUND(1).OR.REC_FOUND(2)) THEN
        NPASS = 3
        RCUT(1) = 1.15
	RCUT(2) = 2.30
	RCUT(3) = 9999.0 ! default values for colour
	! pick up the radii if specified
	IF (REC_NUM(1,1).GT.-1E10) RCUT(1) = REC_NUM(1,1)
	IF (REC_NUM(1,2).GT.-1E10) RCUT(1) = REC_NUM(1,2)
	IF (REC_NUM(2,1).GT.-1E10) RCUT(2) = REC_NUM(2,1)
	IF (REC_NUM(2,2).GT.-1E10) RCUT(2) = REC_NUM(2,2)
      ENDIF

      ! REC_COMMANDS(3) = 'sph'
      IF (REC_FOUND(3)) THEN
        FSPHINPUT = REC_STRING(1, 3)
      ELSE
        WRITE( NOUT, '(A)')
     &'Error you must specify the .sph file as an input: ',
     &'surf2vmd -sph input.sph input.tri output.vmd_plot'
        LERR = .TRUE.
	GOTO 55555

      ENDIF

      ! end of unpacking arguments

      ! find last character in each filename
      CALL CHREND( FOUTPUT, END_FOUTPUT)
      CALL CHREND( FINPUT, END_FINPUT)
      CALL CHREND( FSPHINPUT, END_FSPHINPUT)

! open the input sph  file as readonly
      IF (.NOT.OPENRO( SIN, FINPUT, NOUT)) THEN
        WRITE( NOUT, '(A,A,A)')
     &' Error cannot open input file: ''', FINPUT(1:END_FINPUT), ''''
        LERR = .TRUE. ! have failed
        GOTO 55555
      ENDIF
! open the input sph  file as readonly
      IF (.NOT.OPENRO( SIN2, FSPHINPUT, NOUT)) THEN
        WRITE( NOUT, '(A,A,A)')
     &' Error cannot open input file: ''', FSPHINPUT(1:END_FSPHINPUT),
     &                                                            ''''
        LERR = .TRUE. ! have failed
        GOTO 55555
      ENDIF

! open the output file as new
! use s/r newop - this keeps old version removes tilda's spaces etc. etc.
! arguments output stream, error indicator,
! reduce_output_variable (0 for full), filename, filestream, binary indicator
      IF (FOUTPUT(1:4).NE.'none') THEN
        LBIN = .FALSE.
        CALL NEWOP( NOUT, LERR, 0,
     &                  FOUTPUT, SOUT, LBIN)
        IF (LERR) THEN
          WRITE( NOUT, '(A,A,A)')
     &' Error cannot open output file: ''', FOUTPUT(1:END_FOUTPUT), ''''
          LERR = .TRUE. ! have failed
          GOTO 55555
        ENDIF
!	WRITE( SOUT, '(A)')
!     &'draw delete all ',   ! header
!     &'draw color green '
      ELSE
        ! no output file - only to be used for -peg option
        SOUT = -1
	WRITE( NOUT, '(A)')
     &'WARNING ''none'' specified as output file - so have not opened',
     &'        one.  This only makes sense with -peg option'
      ENDIF

      WRITE( NOUT, '(A,A,A)')
     &'Have opened ', FINPUT(1:END_FINPUT),   ' for surf tri input ',
     &'        and ', FSPHINPUT(1:END_FSPHINPUT), ' for sph input',
     &'        and ', FOUTPUT(1:END_FOUTPUT), ' for output'

C have opened input and output files
C and got initial conditions start work
      ROVER  = 0.
      RLOWER = -1D06
      CALL SPH_PROCESS_READ( NOUT, LERR, SIN2, ROVER, RLOWER,
     &				   SPMAX, SPNO, SPXYZ, SPRAD, SPLAST,
     &				   SPSEC, SPCONN, SPEFFR, LCAPS,
     & 				   SAMPLE, CVECT, PERPVE, PERPVN,
     &				   CONNR, ENDRAD )
      write( *, *) sprad(1)
      IF (LERR) GOTO 55555
      DO SCOUNT = 1, SPNO
        DO TCOUNT = 1, SPNO
	  IF(SCOUNT.NE.TCOUNT) THEN ! look at the pair if different
	    ! Q: is the sphere scount inside sphere tcount
	    DIST = (SPXYZ(1,SCOUNT)-SPXYZ(1,TCOUNT))**2 +
     &             (SPXYZ(2,SCOUNT)-SPXYZ(2,TCOUNT))**2 +
     &             (SPXYZ(3,SCOUNT)-SPXYZ(3,TCOUNT))**2
            IF (DIST.LT.SPRAD(TCOUNT)**2) THEN
	      ! could be inside - centre of scount is in tcount
	      DIST = SQRT(DIST)
	      DIST = DIST+ABS(SPRAD(SCOUNT))
	      IF (DIST.LE.ABS(SPRAD(TCOUNT))) THEN
	        ! inside - mark by chaging case of radius
		SPRAD(SCOUNT) = -ABS(SPRAD(SCOUNT))
	      ENDIF
	    ENDIF

	  ENDIF
	ENDDO
      ENDDO
      DO SCOUNT = 1, SPNO
        WRITE( SOUT,'(4F10.4)')
     &   SPXYZ(1,SCOUNT), SPXYZ(2,SCOUNT), SPXYZ(3,SCOUNT),
     &   SPRAD(SCOUNT)
      ENDDO
C! now we have the .sph file in memory go through the tri file
C10    CONTINUE
C
C        READ( SIN, '(A)', END=55) COL
C        CALL FREDA(1,80,FL,6)
C        IF (KN.EQ.1) THEN
C 	 REC_SPHERE = FN(1)
C        ELSE
C 	 stop
C        ENDIF

C        READ( SIN, '(A)', END=55) COL
C        CALL FREDA(1,80,FL,6)
C        IF (KN.EQ.6) THEN
C 	 TRIA(1)  = FN(1)
C 	 TRIA(2)  = FN(2)
C 	 TRIA(3)  = FN(3)
C 	 NORMA(1) = FN(4)
C 	 NORMA(2) = FN(5)
C 	 NORMA(3) = FN(6)
C        ELSE
C	 stop
C        ENDIF

C	READ( SIN, '(A)', END=55) COL
C	CALL FREDA(1,80,FL,6)
C	IF (KN.EQ.6) THEN
C	  TRIB(1)  = FN(1)
C	  TRIB(2)  = FN(2)
C	  TRIB(3)  = FN(3)
C	  NORMB(1) = FN(4)
C	  NORMB(2) = FN(5)
C	  NORMB(3) = FN(6)
C	ELSE
C	  stop
C	ENDIF
C	CALL FREDA(1,80,FL,6)
C	READ( SIN, '(A)', END=55) COL
C	IF (KN.EQ.6) THEN
C	  TRIC(1)  = FN(1)
C	  TRIC(2)  = FN(2)
C	  TRIC(3)  = FN(3)
C	  NORMC(1) = FN(4)
C	  NORMC(2) = FN(5)
C	  NORMC(3) = FN(6)
C	ELSE
C	  stop
C	ENDIF

C!	WRITE( SOUT, '(12(A,3F11.5) )')
C!     &'draw trinorm { ', TRIA,  ' }  {', TRIB, '}  {', TRIC, ' }  {',
C!     & NORMA,  ' }  {', NORMB, ' }  {',   NORMC, ' }'
C	WRITE( SOUT, '(12(A,3F11.5) )')
C     &'draw triange { ', TRIA,  ' }  {', TRIB, '}  {', TRIC, ' }'
C      GOTO 10
C55    CONTINUE

55555 CONTINUE
      STOP 'FORTRAN STOP - sphqpt normal completion.'
      END