text:
 Users/achash can use this space to introduce him/herself and establish an identity.

20110531 14:12:07 nbsp Welcome to the Wiki!  It looks like you may be associated with the Arlington Farms management.  Please take a few minutes to read over Welcome to the Wiki/Business Owner (which is in the process of being reworked; if you want, have a look at one of the drafts http://daviswiki.org/Welcome_to_the_Wiki/Business_Owner/Talk/Draft_C here or http://daviswiki.org/Welcome_to_the_Wiki/Business_Owner/Talk/Draft_D here... if you have any input on them, it would be very welcome on the Talk page http://daviswiki.org/Welcome_to_the_Wiki/Business_Owner/Talk here!).

The captions you removed were initially written by one user and were restored after you removed them the first time by another user, who wrote I like those caps.  The Wiki is frequently home to quirkiness, as well as to different points of view on exactly how a page should look.  Frankly, those top two pictures on Arlington Farms look freaky.  Im not quite sure how the halo effect got on there, but the color saturation is insane.  Some mention of the fact that the place doesnt actually look like that is warranted.

If you dont like the current captions, you might try writing new ones.  Theres no guarantee someone wont disagree with the change (or come up with a whole new caption), but trying to improve on content is generally seen as better than outright deleting another users contributions. Users/TomGarberson

20110603 16:58:56 nbsp Yes, Im the Assistant Residential Manager for Arlington Farm. Thank you for explaining the changes with the captions! My manager and I will be working on changing them soon to something slightly different :) Users/achash


text:
*
*
      SUBROUTINE PSLASIZESYEV( JOBZ, N, DESCA, MINSIZE )
*
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ
      INTEGER            MINSIZE, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * )
*     ..
*
*  Purpose
*  =======
*
*  PSLASIZESYEV computes the amount of memory needed by PSSYEV
*  to calculate:
*    1)  Eigenvectors and eigenvalues if JOBZ = 'V'
*    2)  Eigenvalues only if JOBZ = 'N'
*
*  Arguments
*  =========
*
*     NP = the number of rows local to a given process.
*     NQ = the number of columns local to a given process.
*
*  JOBZ     (global input) CHARACTER*1
*           Specifies whether or not to compute the eigenvectors:
*           = 'N':  Compute eigenvalues only.
*           = 'V':  Compute eigenvalues and eigenvectors.
*
*  N        (global input) INTEGER
*           Size of the matrix to be tested.  (global size)
*
*  DESCA    (global input) INTEGER array dimension ( DLEN_ )
*
*  MINSIZE  (global output) INTEGER
*           Workspace required for PSSYEV to:
*           1)  Eigenvectors and eigenvalues if JOBZ = 'V'
*           2)  Eigenvalues only if JOBZ = 'N'
*
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
     $                   MB_, NB_, RSRC_, CSRC_, LLD_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
*     ..
*     .. Local Scalars ..
      LOGICAL            WANTZ
      INTEGER            CONTEXTC, CSRC_A, IACOL, IAROW, ICOFFA, IROFFA,
     $                   LCM, LCMQ, LDC, MQ0, MYCOL, MYPCOLC, MYPROWC,
     $                   MYROW, NB, NN, NP, NP0, NPCOL, NPCOLC, NPROCS,
     $                   NPROW, NPROWC, NQ, NRC, QRMEM, RSRC_A,
     $                   SIZEMQRLEFT, SIZEMQRRIGHT
*     ..
*     .. External Functions ..
*
*
      LOGICAL            LSAME
      INTEGER            ILCM, INDXG2P, NUMROC, SL_GRIDRESHAPE
      EXTERNAL           ILCM, INDXG2P, LSAME, NUMROC, SL_GRIDRESHAPE
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, BLACS_GRIDEXIT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*       This is just to keep ftnchek happy
      IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
     $    RSRC_.LT.0 )RETURN
*
      CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
      NB = DESCA( MB_ )
      N = DESCA( M_ )
      RSRC_A = DESCA( RSRC_ )
      CSRC_A = DESCA( CSRC_ )
      LCM = ILCM( NPROW, NPCOL )
      LCMQ = LCM / NPCOL
      IROFFA = 0
      ICOFFA = 0
      IAROW = INDXG2P( 1, NB, MYROW, RSRC_A, NPROW )
      IACOL = INDXG2P( 1, NB, MYCOL, CSRC_A, NPCOL )
      NP = NUMROC( N+IROFFA, NB, MYROW, IAROW, NPROW )
      NQ = NUMROC( N+ICOFFA, NB, MYCOL, IACOL, NPCOL )
      SIZEMQRLEFT = MAX( ( NB*( NB-1 ) ) / 2, ( NP+NQ )*NB ) + NB*NB
      SIZEMQRRIGHT = MAX( ( NB*( NB-1 ) ) / 2,
     $               ( NQ+MAX( NP+NUMROC( NUMROC( N+ICOFFA, NB, 0, 0,
     $               NPCOL ), NB, 0, 0, LCMQ ), NP ) )*NB ) + NB*NB
      NN = MAX( N, NB, 2 )
      NP0 = NUMROC( NN, NB, 0, 0, NPROW )
      MQ0 = NUMROC( NN, NB, 0, 0, NPCOL )
      NPROCS = NPROW*NPCOL
      WANTZ = LSAME( JOBZ, 'V' )
      LDC = 0
*
*     Create the new context that is used in PSSYEV
*
      IF( WANTZ ) THEN
         CONTEXTC = SL_GRIDRESHAPE( DESCA( CTXT_ ), 0, 1, 1, NPROCS, 1 )
         CALL BLACS_GRIDINFO( CONTEXTC, NPROWC, NPCOLC, MYPROWC,
     $                        MYPCOLC )
         NRC = NUMROC( N, NB, MYPROWC, 0, NPROCS)
         LDC = MAX( 1, NRC )
         CALL BLACS_GRIDEXIT( CONTEXTC )
      END IF

*
*     Compute the total amount of space needed
*
      IF( WANTZ ) THEN
        QRMEM = 5*N + MAX(  2*NP0 +MQ0 + NB*NN, 2*NN-2 ) + N*LDC
        MINSIZE = MAX ( SIZEMQRLEFT, SIZEMQRRIGHT, QRMEM )
      ELSE
         MINSIZE = 5*N + 2*NP0 +MQ0 + NB*NN
      END IF
*
      RETURN
*
*     End of PSLASIZESYEV
*
      END

ext: f

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
!==========================================================================
elemental function gsw_t_freezing_exact (sa, p, saturation_fraction)
!==========================================================================
!
!  Calculates the in-situ temperature at which seawater freezes. The
!  in-situ temperature freezing point is calculated from the exact
!  in-situ freezing temperature which is found by a modified Newton-Raphson
!  iteration (McDougall and Wotherspoon, 2013) of the equality of the
!  chemical potentials of water in seawater and in ice.
!
!  An alternative GSW function, gsw_t_freezing_poly, it is based on a
!  computationally-efficient polynomial, and is accurate to within -5e-4 K
!  and 6e-4 K, when compared with this function.
!
!  SA  =  Absolute Salinity                                        [ g/kg ]
!  p   =  sea pressure                                             [ dbar ]
!         ( i.e. absolute pressure - 10.1325 dbar )
!  saturation_fraction = the saturation fraction of dissolved air in
!                        seawater
!  (i.e., saturation_fraction must be between 0 and 1, and the default
!    is 1, completely saturated)
!
!  t_freezing = in-situ temperature at which seawater freezes.    [ deg C ]
!--------------------------------------------------------------------------

use gsw_mod_teos10_constants, only : gsw_sso

use gsw_mod_toolbox, only : gsw_gibbs_ice, gsw_chem_potential_water_t_exact
use gsw_mod_toolbox, only : gsw_t_deriv_chem_potential_water_t_exact
use gsw_mod_toolbox, only : gsw_t_freezing_poly

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: sa, p, saturation_fraction

real (r8) :: gsw_t_freezing_exact

real (r8) :: df_dt, p_r, sa_r, tf, tfm, tf_old, x, f

! The initial value of t_freezing_exact (for air-free seawater)
tf = gsw_t_freezing_poly(sa,p,polynomial=.true.)

df_dt = 1e3_r8*gsw_t_deriv_chem_potential_water_t_exact(sa,tf,p) - &
		gsw_gibbs_ice(1,0,tf,p)
! df_dt here is the initial value of the derivative of the function f whose
! zero (f = 0) we are finding (see Eqn. (3.33.2) of IOC et al (2010)).

tf_old = tf
f = 1e3_r8*gsw_chem_potential_water_t_exact(sa,tf_old,p) - &
		gsw_gibbs_ice(0,0,tf_old,p)
tf = tf_old - f/df_dt
tfm = 0.5_r8*(tf + tf_old)
df_dt = 1e3_r8*gsw_t_deriv_chem_potential_water_t_exact(sa,tfm,p) - &
		gsw_gibbs_ice(1,0,tfm,p)
tf = tf_old - f/df_dt

tf_old = tf
f = 1e3_r8*gsw_chem_potential_water_t_exact(sa,tf_old,p) - &
		gsw_gibbs_ice(0,0,tf_old,p)
tf = tf_old - f/df_dt

! Adjust for the effects of dissolved air
gsw_t_freezing_exact = tf - &
                saturation_fraction*(1e-3_r8)*(2.4_r8 - sa/(2.0_r8*gsw_sso))

return
end function

!--------------------------------------------------------------------------

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine verlet_bias: Changed version of verlet.f90 (dynamic), here
!        with incorporation of an additional bias potential for umbrella
!        samplings within rpmd.f90
!
!     part of EVB
!
subroutine verlet (istep,dt,derivs,epot,ekin,afm_force,analyze)
use general
use evb_mod
use qmdff
!  ambigious reference of z-coord to qmdff z(94) array!!
implicit none
integer::i,j,k,l,istep   ! loop indices and actual step number
real(kind=8)::dt,dt_2
real(kind=8)::etot,epot,ekin,epot1
real(kind=8)::eksum
logical::analyze  ! if temperature etc shall be calculated
real(kind=8)::temp,pres
real(kind=8),dimension(natoms)::charge  !dummy array for charges
!   for umbrella sampling
real(kind=8),dimension(3,natoms)::derivs_1d,q_1b
real(kind=8),dimension(3,natoms,nbeads)::derivs
real(kind=8),dimension(nat6)::act_int,int_ideal
real(kind=8)::poly(4,nbeads)
real(kind=8)::beta_n    ! the inverse temperature per bead
real(kind=8)::twown   ! inverse inverse temperature
real(kind=8)::pi_n    ! half circle section for each bead
real(kind=8)::wk,wt,wm  ! circle section scaled parameters
real(kind=8)::cos_wt,sin_wt   ! sinus and cosinus value of the circle section
!real(kind=8)::pi   ! the pi
real(kind=8)::infinity,dbl_prec_var  ! test if one of the coordinates is infinity
real(kind=8)::p_new ! the actual new momentum/coordinate component
real(kind=8),allocatable::q_old(:,:,:) ! for hard boxes and RPMD: store actual coord
real(kind=8)::centroid(3,natoms)  !the centroid with all the com's of all beads
integer::round  ! the actual window to sample (for calculation of avg/variance)
integer::constrain  ! if the trajectory shall be constrained to the actual xi value
                    ! (only for recrossing calculations)
                    !  0 : usual umbrella sampling
                    !  1 : constrain to dividing surface
                    !  2 : child trajectory: no restrain at all
                    !  3 : for pre sampling with tri/tetramolecular reactions
integer::num,modul
!   for AFM simulation
real(kind=8)::move_act(3),move_shall(3)
real(kind=8)::ekin1,ekin2
real(kind=8)::afm_bias(3),afm_force,dist
!   for kinetic energy
real(kind=8)::act_temp
!   for Nose-Hoover thermostat
real(kind=8)::massvec(3,natoms,nbeads)  ! vector with atom masses per dimension
real(kind=8)::ndoft ! number degrees of freedom
real(kind=8)::nose_zeta2
integer::tries  ! number of periodic corrections before throwing an error
integer::j_lower,j_upper
!   for Nose-Hoover barostat
real(kind=8)::e2,e4,e6,e8
real(kind=8)::eterm2,term,term2,resize,expterm
real(kind=8)::stress_ten(3,3),ekin_ten(3,3),factor
real(kind=8)::act_vol,act_dens
!   for Berendsen barostat
real(kind=8)::beren_scale
real(kind=8)::press_avg  ! the average pressure for printout
!parameter(pi=3.1415926535897932384626433832795029d0)

infinity = HUGE(dbl_prec_var)   ! set the larges possible real value

!
!     Calculate new box dimensions, if the NPT ensemble is used
!
if (periodic) then
   volbox=boxlen_x*boxlen_y*boxlen_z
end if
!
!     Reset the virial tensor if a pressure is applied
!
if (npt) then
   calc_vir=.true.
   vir_ten=0.d0
end if
!
!     For the first timestep: set the pressure to the desired value
!
if (istep .eq. 1) press_act=pressure

call get_centroid(centroid)

do k=1,nbeads
   do i=1,natoms
      do j=1,3
         massvec(j,i,k)=mass(i)
      end do
   end do
end do
!
!     For the Nose-Hoover chain thermostat: apply the full chain on the momenta
!     on the half-time step
!
if (thermostat .eq. 2) then
!
!     The NPT ensemble: correct momenta and pressure
!
   if (npt) then
      if (barostat .eq. 2) then
         call nhc_npt(dt,centroid,volbox)
      end if
!      call berendsen(dt,centroid,volbox)
!
!     The NVT ensemble: correct only momenta
!
   else
      call nhc(dt,centroid)
   end if
end if
!
!     update the momentum (half time step)
!     ---> (Nose-Hoover step D: second part of half time momentum update)
!
!if (thermostat .eq. 1 .or. thermostat .eq. 2) then

p_i=p_i-0.5d0*dt*derivs
!end if
!

!
!     calculate averaged kinetic energy for subset of the system
!
if (calc_ekin) then
!
!     first, calculate the kinetic energy by applying the centroid-momenta approximation
!
   ekin1=0d0
   do i=1,ekin_num
      k=ekin_atoms(i)
      do j=1,nbeads
       !  write(677,*) i,j,p_i(:,i,j),mass(i)
         ekin1=ekin1+dot_product(p_i(:,k,j),p_i(:,k,j))/(2d0*mass(k))
      end do
   end do
   ekin_avg=ekin_avg+ekin1

!
!     second, calculate the kinetic energy by applying the virial theorem
!
   ekin2=0d0
   beta_n = beta / nbeads
   do i=1,ekin_num
      k=ekin_atoms(i)
      do j=1,nbeads

         ekin2=ekin2+dot_product(q_i(:,k,j)-centroid(:,k),derivs(:,k,j))!+&
                    !  & mass(i)/beta_n**2/nbeads*((2*q_i(:,i,j)-q_i(:,i,j_upper)-q_i(:,i,j_lower))))

      end do
   end do
   ekin2=-ekin2/(2.d0*nbeads)
   ekin2=ekin2+ekin_num*3d0/(2d0*beta)
   ekin2_avg=ekin2_avg+ekin2
   write(156,*) ekin1/(nbeads**2*ekin_num),ekin2/ekin_num
end if
!
!     for AFM simulatons: fix only the anchor atom!
!
if (afm_run) then
   p_i(:,afm_fix_at,:)=0.d0
end if
!
!     If a barostat is used: prepare the volume update
!
if (npt) then
   if (barostat .eq. 2) then
      term = vbar*0.5d0*dt
      term2 = term*term
      expterm=exp(term)
      eterm2 = expterm * expterm
      e2 = 1.0d0 / 6.0d0
      e4 = e2 / 20.0d0
      e6 = e4 / 42.0d0
      e8 = e6 / 72.0d0
      resize = 1.0d0 + term2*(e2+term2*(e4+term2*(e6+term2*e8)))
      resize = expterm * resize * dt
   end if
end if
!
!     update the positions: for one bead, do the usual verlet procedure
!
!     ---> Nose-Hoover step E
if (nbeads .eq. 1) then
!
!     If a system with hard box walls is simulated: check, if one of the atoms
!     would move outside the box this timestep; if this is the case, revert
!     the momentum in the respective axis
!
   if (box_walls) then
      do i=1,natoms
         do j=1,3
            if (q_i(1,i,1)+p_i(1,i,1)*dt/massvec(1,i,1) .ge. walldim_x) then
               p_i(1,i,1)=-p_i(1,i,1)
            else if (q_i(2,i,1)+p_i(2,i,1)*dt/massvec(2,i,1) .ge. walldim_y) then
               p_i(2,i,1)=-p_i(2,i,1)
            else if (q_i(3,i,1)+p_i(3,i,1)*dt/massvec(3,i,1) .ge. walldim_z) then
               p_i(3,i,1)=-p_i(3,i,1)
            else if (q_i(j,i,1)+p_i(j,i,1)*dt/massvec(j,i,1) .le. 0.d0) then
               p_i(j,i,1)=-p_i(j,i,1)
            end if
         end do
      end do
   end if
!
!    For NPT (Nose-Hoover) or NVT/NVE ensembles
!

   if (npt) then
      if (barostat .eq. 2) then
         q_i=q_i*eterm2+p_i/massvec*resize
      else
         q_i=q_i+p_i*dt/massvec
      end if
   else
      q_i=q_i+p_i*dt/massvec
   end if
else
!
!     For the ring polymer: calculate the harmonic free ring polymer
!     interactions between the beats: do it in normal mode space
!     tranform with Fast Fourier transformation to the normal mode
!     space and back thereafter to the cartesian space
!
!     Transform to normal mode space
!     --> What is done there, exactly??
!
!     For hard boxes: store the actual coordinates, if momenta shall be reversed
   if (box_walls) then
      allocate(q_old(3,natoms,nbeads))
      q_old=q_i
   end if

   do i = 1, 3
      do j = 1, Natoms
!
!     For periodic systems: asure that all replicas of one atoms
!     fulfil the minimum image convention! No single replicas moved
!     on the other side of the box, generating a giant distance..
!     --> possible violations of the box size will be corrected
!      automatically in the next section!
!

!         if (periodic) then
!            tries=0
!            do k=2,nbeads
!               do while (abs(q_i(i,j,1)-q_i(i,j,k)) .gt. box_len2)
!                  q_i(i,j,k)=q_i(i,j,k)+sign(box_len,q_i(i,j,1)-q_i(i,j,k))
!                  tries=tries+1
!                  if (tries .gt. 20) then
!                     write(*,*) "Too many correction steps needed for periodic dynamics!"
!                     write(*,*) "The system seems to be exploded! Check your settings!"
!                     call fatal
!                  end if
!               end do
!            end do
!         end if
         call rfft(p_i(i,j,:), Nbeads)
         call rfft(q_i(i,j,:), Nbeads)
      end do
   end do
   do j = 1, Natoms

      poly(1,1) = 1.0d0
      poly(2,1) = 0.0d0
      poly(3,1) = dt / mass(j)
      poly(4,1) = 1.0d0

      if (Nbeads .gt. 1) then
         beta_n = beta / nbeads
         twown = 2.0d0 / beta_n
         pi_n = pi / Nbeads
         do k = 1, nbeads / 2
            wk = twown * dsin(k * pi_n)
            wt = wk * dt
            wm = wk * mass(j)
            cos_wt = dcos(wt)
            sin_wt = dsin(wt)
            poly(1,k+1) = cos_wt
            poly(2,k+1) = -wm*sin_wt
            poly(3,k+1) = sin_wt/wm
            poly(4,k+1) = cos_wt
         end do
         do k = 1, (Nbeads - 1) / 2
            poly(1,Nbeads-k+1) = poly(1,k+1)
            poly(2,Nbeads-k+1) = poly(2,k+1)
            poly(3,Nbeads-k+1) = poly(3,k+1)
            poly(4,Nbeads-k+1) = poly(4,k+1)
         end do
      end if

      do k = 1, Nbeads
         do i = 1, 3
            p_new = p_i(i,j,k) * poly(1,k) + q_i(i,j,k) * poly(2,k)
!
!    For NPT ensemble: resize the box dimensions!
!
            if (npt) then
               if (barostat .eq. 2) then
                  q_i(i,j,k) = p_i(i,j,k) * poly(3,k)*resize + q_i(i,j,k) *  &
                           & poly(4,k) * eterm2
               else
                  q_i(i,j,k) = p_i(i,j,k) * poly(3,k) + q_i(i,j,k) * poly(4,k)
               end if
            else
               q_i(i,j,k) = p_i(i,j,k) * poly(3,k) + q_i(i,j,k) * poly(4,k)
            end if
            p_i(i,j,k) = p_new
         end do
      end do
   end do
!
!     Transform back to Cartesian space
!
   do i = 1, 3
      do j = 1, Natoms
         call irfft(p_i(i,j,:), Nbeads)
         call irfft(q_i(i,j,:), Nbeads)
      end do
   end do
!
!     If a system with hard box walls is simulated: check, if one of the atoms
!     has moved outside the box this timestep; if this is the case, revert
!     the momentum in the respective axis
!
   if (box_walls) then
      do k=1,nbeads
         do i=1,natoms
            do j=1,3
               if (q_i(1,i,1) .ge. walldim_x) then
                  q_i(1,i,k)=2.d0*q_old(1,i,k)-q_i(1,i,k)
                  p_i(1,i,k)=-p_i(1,i,k)
               else if (q_i(2,i,1) .ge. walldim_y) then
                  q_i(2,i,k)=2.d0*q_old(2,i,k)-q_i(2,i,k)
                  p_i(2,i,k)=-p_i(2,i,k)
               else if (q_i(3,i,1) .ge. walldim_z) then
                  q_i(3,i,k)=2.d0*q_old(3,i,k)-q_i(3,i,k)
                  p_i(3,i,k)=-p_i(3,i,k)
               else if (q_i(j,i,1)  .le. 0.d0) then
                  q_i(j,i,k)=q_old(j,i,k)-2.d0*q_i(j,i,k)
                  p_i(j,i,k)=-p_i(j,i,k)
               end if
            end do
         end do
      end do
      deallocate(q_old)
   end if

end if
!
!     For Nose-Hoover barostat in NPT ensemble, align the dimensions of the box
!
if (npt) then
   if (barostat .eq. 2) then
      boxlen_x=boxlen_x*eterm2
      boxlen_y=boxlen_y*eterm2
      boxlen_z=boxlen_z*eterm2
      boxlen_x2=boxlen_x*0.5d0
      boxlen_y2=boxlen_y*0.5d0
      boxlen_z2=boxlen_z*0.5d0
      volbox=boxlen_x*boxlen_y*boxlen_z
   end if
end if

!
!     For periodic systems: shift all atoms that were moved outside the
!     box on the other side!  If several RPMD beads are involved, change the
!     coordinates of all of them in the same way
!
if (periodic) then
   do i=1,nbeads
      do j=1,natoms
         do k=1,3
            tries=0
            do while (q_i(k,j,i) .lt. 0)
               if (k .eq. 1) then
                  q_i(k,j,:)=q_i(k,j,:)+boxlen_x
               else if (k .eq. 2) then
                  q_i(k,j,:)=q_i(k,j,:)+boxlen_y
               else if (k .eq. 3) then
                  q_i(k,j,:)=q_i(k,j,:)+boxlen_z
               end if
               tries=tries+1
               if (tries .gt. 20) then
                  write(*,*) "Too many correction steps needed for periodic dynamics! (lower)"
                  write(*,*) "The system seems to be exploded! Check your settings!"
                  call fatal
               end if
            end do
            if (k .eq. 1) then
               do while (q_i(1,j,i) .gt. boxlen_x)
                  q_i(1,j,:)=q_i(1,j,:)-boxlen_x
                  tries=tries+1
                  if (tries .gt. 20) then
                     write(*,*) "Too many correction steps needed for periodic dynamics! (x, upper)"
                     write(*,*) "The system seems to be exploded! Check your settings!"
                     call fatal
                  end if
               end do
            else if (k .eq. 2) then
               do while (q_i(2,j,i) .gt. boxlen_y)
                  q_i(2,j,:)=q_i(2,j,:)-boxlen_y
                  tries=tries+1
                  if (tries .gt. 20) then
                     write(*,*) "Too many correction steps needed for periodic dynamics! (y,upper)"
                     write(*,*) "The system seems to be exploded! Check your settings!"
                     call fatal
                  end if
               end do
            else if (k .eq. 3) then
               do while (q_i(3,j,i) .gt. boxlen_z)
                  q_i(3,j,:)=q_i(3,j,:)-boxlen_z
                  tries=tries+1
                  if (tries .gt. 20) then
                     write(*,*) "Too many correction steps needed for periodic dynamics! (z,upper)"
                     write(*,*) "The system seems to be exploded! Check your settings!"
                     call fatal
                  end if
               end do
            end if
         end do
      end do
   end do
end if
!
!     get the potential energy and atomic forces
!     for each bead at once: define its current structure and collect
!     all results in a global derivs array thereafter
!
epot=0.d0
if (energysplit) then
   e_cov_split=0.d0
   e_noncov_split=0.d0
end if
do i=1,nbeads
   q_1b=q_i(:,:,i)
   call gradient (q_1b,epot1,derivs_1d,i)
   derivs(:,:,i)=derivs_1d
   epot=epot+epot1
end do
!
!   write out the gradients to file (verbose)
!
if (verbose) then
   if (mod(istep,iwrite) .eq. 0) then
      write(29,*) natoms*nbeads
      write(29,*)
      do k=1,nbeads
         do i=1,natoms
            write(29,*) i,k,derivs(:,i,k)
         end do
      end do
   end if

!
!    write out the velocities to file (verbose)
!

   if (mod(istep,iwrite) .eq. 0) then
      write(30,*) natoms*nbeads
      write(30,*)
      do k=1,nbeads
         do i=1,natoms
            write(30,*) i,k,p_i(:,i,k)/mass(i)
         end do
      end do
   end if
end if
!
!     For Mechanochemistry calculations: Add the additional forces here
!     to the gradient of the reference method!
!

if (add_force) then
   do i=1,nbeads
      derivs(:,force1_at,i)=derivs(:,force1_at,i)+force1_v(:)*force1_k
      derivs(:,force2_at,i)=derivs(:,force2_at,i)+force2_v(:)*force2_k
   end do
end if

!
!     For an AFM simulation run: determine the current position of the
!     bias force and its vector to be applied to the system (and in order
!     to calculate the effective AFM force)
!

if (afm_run) then
   call get_centroid(centroid)
   move_act=centroid(:,afm_move_at)
   move_shall=afm_move_first(:)+afm_move_v(:)*afm_move_dist*(real(istep)/real(afm_steps))

   afm_bias=(move_act-move_shall)*afm_k

   afm_force=sqrt(afm_bias(1)**2+afm_bias(2)**2+afm_bias(3)**2)

   do i=1,nbeads
      derivs(:,afm_move_at,i)=derivs(:,afm_move_at,i)+afm_bias
   end do
   afm_force=afm_force/newton2au
end if
!
!     update the momentum (full time step)
!
!     ---> (Nose-Hoover step F: first part of full time momentum update)
!
p_i=p_i-0.5d0*dt*derivs

!
!     For the Nose-Hoover chain thermostat: apply the full chain on the momenta
!     on the full-time step
!
call get_centroid(centroid)
if (thermostat .eq. 2) then
!
!     The NPT ensemble: correct momenta and pressure
!
   if (npt) then
      if (barostat .eq. 2) then
         call nhc_npt(dt,centroid,volbox)
      end if
!
!     The NVT ensemble: correct only momenta
!
   else
      call nhc(dt,centroid)
   end if
end if
!
!     Calculate the kinetic energy part of the stress tensor
!
if (npt) then
   ekin_ten=0.d0
   do i=1,nbeads
      do j=1,natoms
         do k=1,3
            do l=1,3
               ekin_ten(k,l)=ekin_ten(k,l)*p_i(k,j,i)*p_i(l,j,i)/mass(j)
            end do
         end do
      end do
   end do
!
!     Now calculate the actual pressure from the virial tensor / stress tensor
!
   factor = 1.d0 / volbox
   do i=1,3
      do j=1,3
         stress_ten(j,i)=factor*(2d0*ekin_ten(j,i)-vir_ten(j,i))
      end do
   end do
   press_act=(stress_ten(1,1)+stress_ten(2,2)+stress_ten(3,3)) / 3.d0
end if
!
!     If the Berendsen barostat is chosen, rescale the box and the coordinates
!
if (npt) then
  if (barostat .eq. 1) then
     beren_scale=(1.d0+(dt*compress/taupres)*(press_act-pressure))**(1.d0/3.d0)
  !   write(*,*) "scale",beren_scale,(dt*compress/taupres),(press_act-pressure)
     boxlen_x=boxlen_x*beren_scale
     boxlen_y=boxlen_y*beren_scale
     boxlen_z=boxlen_z*beren_scale
     boxlen_x2=boxlen_x*0.5d0
     boxlen_y2=boxlen_y*0.5d0
     boxlen_z2=boxlen_z*0.5d0
     volbox=boxlen_x*boxlen_y*boxlen_z
     q_i=q_i*beren_scale

  end if
end if
!
!     Apply andersen thermostat to apply an random hit and
!     change the momentum (only every andersen_step steps!)
!     --> Replace momenta with a fresh sampling from a Gaussian
!     distribution at the temperature of interest
!
if (thermostat .eq. 0) then
   if (mod(istep,andersen_step) .eq. 0) then
      call andersen
   end if
end if
!
!     calculate temperature and other dynamical parameters for each dump step
!
press_avg=press_avg+press_act
if (analyze) then
   ekin1=0.d0
   ekin=0.d0
   do i=1,natoms
      do j=1,nbeads
       !  write(677,*) i,j,p_i(:,i,j),mass(i)
         ekin1=ekin1+dot_product(p_i(:,i,j),p_i(:,i,j))/(2d0*mass(i))
      end do
   end do
   ekin=ekin1
   act_temp=2d0*ekin1/3d0/0.316679D-5/natoms/nbeads/nbeads
   if (npt) then
      act_vol=volbox*bohr**3
      act_dens=(mass_tot*1.6605402E-24*emass)/(act_vol*1E-24)

      write(*,'(i8,f12.5,a,f10.4,a,f12.4,f14.5,a,f12.6)') istep,epot,"     ",act_temp,"  ",act_vol, &
               & press_avg*prescon/iwrite,"    ",act_dens
   else
      write(*,'(a,i8,a,f12.5,a,f10.4,a)') " Step: ",istep,"  --  pot. energy: ", &
                & epot," Hartree  --  temperature: ",act_temp," K"
   end if
   write(236,*) istep,act_temp
   temp_test=temp_test+act_temp
   press_avg=0.d0
end if

!
!     If an error occured in the trajectory and one of the entries is either NaN
!     or Infinity, the error variable is set to 1
!
!do i=1,3
!   do j=1,natoms
!      do k=1,nbeads
!         if (q_i(i,j,k) .ne. q_i(i,j,k) .or. q_i(i,j,k) .gt. infinity) then
!            traj_error=1
!         end if
!      end do
!   end do
!end do
!write(*,*) "q_i",q_i
!write(*,*) "p_i",p_i
!write(*,*) "derivs",derivs

!
!     Remove total translation and rotation of the system during the dynamics
!      (especially needed for Nose-Hoover thermostat!)
!

call transrot(centroid)

return
end subroutine verlet

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


!> @brief This module contains the storage package methods
!!
!! This module contains the methods used to add the effects of storage
!! on the groundwater flow equation. The contribution of specific
!! storage and specific yield can be represented.
!!
!<
module GwfStoModule

  use KindModule, only: DP, I4B, LGP
  use ConstantsModule, only: DZERO, DEM6, DEM4, DHALF, DONE, DTWO,               &
                             LENBUDTXT, LINELENGTH
  use SimVariablesModule, only: errmsg
  use SimModule, only: store_error, count_errors
  use SmoothingModule, only: sQuadraticSaturation, &
                             sQuadraticSaturationDerivative, &
                             sQSaturation, sLinearSaturation
  use BaseDisModule, only: DisBaseType
  use NumericalPackageModule, only: NumericalPackageType
  use BlockParserModule, only: BlockParserType
  use GwfStorageUtilsModule, only: SsCapacity, SyCapacity, SsTerms, SyTerms
  use InputOutputModule, only: GetUnit, openfile
  use TvsModule, only: TvsType, tvs_cr

  implicit none
  public :: GwfStoType, sto_cr

  character(len=LENBUDTXT), dimension(2) :: budtxt = & !< text labels for budget terms
                                            ['          STO-SS', '          STO-SY']

  type, extends(NumericalPackageType) :: GwfStoType
    integer(I4B), pointer                            :: istor_coef => null()     !< indicates if ss is the storage coefficient
    integer(I4B), pointer                            :: iconf_ss => null()       !< indicates if ss is 0 below the top of a layer
    integer(I4B), pointer                            :: iorig_ss => null()       !< indicates if the original storage specific storage formulation should be used
    integer(I4B), pointer                            :: iss => null()            !< steady state flag: 1 = steady, 0 = transient
    integer(I4B), pointer                            :: iusesy => null()         !< flag set if any cell is convertible (0, 1)
    integer(I4B), dimension(:), pointer, contiguous  :: iconvert => null()       !< confined (0) or convertible (1)
    real(DP), dimension(:), pointer, contiguous      :: ss => null()             !< specfic storage or storage coefficient
    real(DP), dimension(:), pointer, contiguous      :: sy => null()             !< specific yield
    real(DP), dimension(:), pointer, contiguous      :: strgss => null()         !< vector of specific storage rates
    real(DP), dimension(:), pointer, contiguous      :: strgsy => null()         !< vector of specific yield rates
    integer(I4B), dimension(:), pointer, contiguous  :: ibound => null()         !< pointer to model ibound
    real(DP), pointer                                :: satomega => null()       !< newton-raphson saturation omega
    integer(I4B), pointer                         :: integratechanges => null()  !< indicates if mid-simulation ss and sy changes should be integrated via an additional matrix formulation term
    integer(I4B), pointer                            :: intvs => null()          !< TVS (time-varying storage) unit number (0 if unused)
    type(TvsType), pointer                           :: tvs => null()            !< TVS object
    real(DP), dimension(:), pointer, contiguous, private :: oldss => null()      !< previous time step specific storage
    real(DP), dimension(:), pointer, contiguous, private :: oldsy => null()      !< previous time step specific yield
  contains
    procedure :: sto_ar
    procedure :: sto_rp
    procedure :: sto_ad
    procedure :: sto_fc
    procedure :: sto_fn
    procedure :: sto_cq
    procedure :: sto_bd
    procedure :: sto_save_model_flows
    procedure :: sto_da
    procedure          :: allocate_scalars
    procedure, private :: allocate_arrays
    !procedure, private :: register_handlers
    procedure, private :: read_options
    procedure, private :: read_data
    procedure, private :: save_old_ss_sy
  end type

contains

  !> @ brief Create a new package object
  !!
  !!  Create a new storage (STO) object
  !!
  !<
  subroutine sto_cr(stoobj, name_model, inunit, iout)
    ! -- dummy variables
    type(GwfStoType), pointer :: stoobj         !< GwfStoType object
    character(len=*), intent(in) :: name_model  !< name of model
    integer(I4B), intent(in) :: inunit          !< package input file unit
    integer(I4B), intent(in) :: iout            !< model listing file unit
    !
    ! -- Create the object
    allocate (stoobj)
    !
    ! -- create name and memory path
    call stoobj%set_names(1, name_model, 'STO', 'STO')
    !
    ! -- Allocate scalars
    call stoobj%allocate_scalars()
    !
    ! -- Set variables
    stoobj%inunit = inunit
    stoobj%iout = iout
    !
    ! -- Initialize block parser
    call stoobj%parser%Initialize(stoobj%inunit, stoobj%iout)
    !
    ! -- return
    return
  end subroutine sto_cr

  !> @ brief Allocate and read method for package
  !!
  !!  Method to allocate and read static data for the STO package.
  !!
  !<
  subroutine sto_ar(this, dis, ibound)
    ! -- modules
    use MemoryManagerModule, only: mem_setptr
    use MemoryHelperModule, only: create_mem_path
    ! -- dummy variables
    class(GwfStoType)                       :: this            !< GwfStoType object
    class(DisBaseType), pointer, intent(in) :: dis             !< model discretization object
    integer(I4B), dimension(:), pointer, contiguous :: ibound  !< model ibound array
    ! -- local variables
    ! -- formats
    character(len=*), parameter :: fmtsto = &
      "(1x,/1x,'STO -- STORAGE PACKAGE, VERSION 1, 5/19/2014',                 &
      &' INPUT READ FROM UNIT ', i0, //)"
    !
    ! --print a message identifying the storage package.
    write (this%iout, fmtsto) this%inunit
    !
    ! -- store pointers to arguments that were passed in
    this%dis => dis
    this%ibound => ibound
    !
    ! -- set pointer to gwf iss
    call mem_setptr(this%iss, 'ISS', create_mem_path(this%name_model))
    !
    ! -- Allocate arrays
    call this%allocate_arrays(dis%nodes)
    !!
    !! -- Register side effect handlers
    !call this%register_handlers()
    !
    ! -- Read storage options
    call this%read_options()
    !
    ! -- read the data block
    call this%read_data()
    !
    ! -- TVS
    if (this%intvs /= 0) then
      call this%tvs%ar(this%dis)
    end if
    !
    ! -- return
    return
  end subroutine sto_ar

  !> @ brief Read and prepare method for package
  !!
  !!  Method to read and prepare stress period data for the STO package.
  !!
  !<
  subroutine sto_rp(this)
    ! -- modules
    use TdisModule, only: kper, nper
    implicit none
    ! -- dummy variables
    class(GwfStoType) :: this  !< GwfStoType object
    ! -- local variables
    integer(I4B) :: ierr
    logical :: isfound, readss, readsy, endOfBlock
    character(len=16) :: css(0:1)
    character(len=LINELENGTH) :: line, keyword
    ! -- formats
    character(len=*), parameter :: fmtlsp = &
                                   "(1X,/1X,'REUSING ',A,' FROM LAST STRESS PERIOD')"
    character(len=*), parameter :: fmtblkerr = &
                                   "('Error.  Looking for BEGIN PERIOD iper.  Found ', a, ' instead.')"
    ! -- data
    data css(0)/'       TRANSIENT'/
    data css(1)/'    STEADY-STATE'/
! ------------------------------------------------------------------------------
    !
    ! -- Store ss and sy values from end of last stress period if needed
    if(this%integratechanges /= 0) then
      call this%save_old_ss_sy()
    end if
    !
    ! -- get stress period data
    if (this%ionper < kper) then
      !
      ! -- get period block
      call this%parser%GetBlock('PERIOD', isfound, ierr, &
                                supportOpenClose=.true.)
      if (isfound) then
        !
        ! -- read ionper and check for increasing period numbers
        call this%read_check_ionper()
      else
        !
        ! -- PERIOD block not found
        if (ierr < 0) then
          ! -- End of file found; data applies for remainder of simulation.
          this%ionper = nper + 1
        else
          ! -- Found invalid block
          call this%parser%GetCurrentLine(line)
          write (errmsg, fmtblkerr) adjustl(trim(line))
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
      end if
    end if
    !
    ! -- read data if ionper == kper
    ! are these here to anticipate reading ss,sy per stress period?
    readss = .false.
    readsy = .false.

    !stotxt = aname(2)
    if (this%ionper == kper) then
      write (this%iout, '(//,1x,a)') 'PROCESSING STORAGE PERIOD DATA'
      do
        call this%parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('STEADY-STATE')
          this%iss = 1
        case ('TRANSIENT')
          this%iss = 0
        case default
          write (errmsg, '(4x,a,a)') 'Unknown STORAGE data tag: ', &
            trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END PROCESSING STORAGE PERIOD DATA'
      !else
      !  write(this%iout,fmtlsp) 'STORAGE VALUES'
    end if

    write (this%iout, '(//1X,A,I0,A,A,/)') &
      'STRESS PERIOD ', kper, ' IS ', trim(adjustl(css(this%iss)))
    !
    ! -- TVS
    if (this%intvs /= 0) then
      call this%tvs%rp()
    end if
    !
    ! -- return
    return
  end subroutine sto_rp

  !> @ brief Advance the package
  !!
  !!  Advance data in the STO package.
  !!
  !<
  subroutine sto_ad(this)
    ! -- modules
    use TdisModule, only: kstp
    ! -- dummy variables
    class(GwfStoType) :: this  !< GwfStoType object
    !
    ! -- Store ss and sy values from end of last time step if needed
    if(this%integratechanges /= 0 .and. kstp > 1) then
      call this%save_old_ss_sy()
    end if
    !
    ! -- TVS
    if(this%intvs /= 0) then
      call this%tvs%ad()
    endif
    !
    ! -- return
    return
  end subroutine sto_ad

  !> @ brief Fill A and right-hand side for the package
  !!
  !!  Fill the coefficient matrix and right-hand side with the STO package terms.
  !!
  !<
  subroutine sto_fc(this, kiter, hold, hnew, njasln, amat, idxglo, rhs)
    ! -- modules
    use TdisModule, only: delt
    ! -- dummy variables
    class(GwfStoType) :: this                           !< GwfStoType object
    integer(I4B), intent(in) :: kiter                   !< outer iteration number
    real(DP), intent(in), dimension(:) :: hold          !< previous heads
    real(DP), intent(in), dimension(:) :: hnew          !< current heads
    integer(I4B), intent(in) :: njasln                  !< size of the A matrix for the solution
    real(DP), dimension(njasln), intent(inout) :: amat  !< A matrix
    integer(I4B), intent(in), dimension(:) :: idxglo    !< global index model to solution
    real(DP), intent(inout), dimension(:) :: rhs        !< right-hand side
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: tled
    real(DP) :: sc1
    real(DP) :: sc2
    real(DP) :: rho1
    real(DP) :: rho2
    real(DP) :: sc1old
    real(DP) :: sc2old
    real(DP) :: rho1old
    real(DP) :: rho2old
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: snold
    real(DP) :: snnew
    real(DP) :: aterm
    real(DP) :: rhsterm
    ! -- formats
    character(len=*), parameter :: fmtsperror =                                &
      &"('DETECTED TIME STEP LENGTH OF ZERO.  GWF STORAGE PACKAGE CANNOT BE ', &
      &'USED UNLESS DELT IS NON-ZERO.')"
    !
    ! -- test if steady-state stress period
    if (this%iss /= 0) return
    !
    ! -- Ensure time step length is not zero
    if (delt == DZERO) then
      write (errmsg, fmtsperror)
      call store_error(errmsg, terminate=.TRUE.)
    end if
    !
    ! -- set variables
    tled = DONE/delt
    !
    ! -- loop through and calculate storage contribution to hcof and rhs
    do n = 1, this%dis%nodes
      idiag = this%dis%con%ia(n)
      if (this%ibound(n) < 1) cycle
      !
      ! -- aquifer elevations and thickness
      tp = this%dis%top(n)
      bt = this%dis%bot(n)
      !
      ! -- aquifer saturation
      if (this%iconvert(n) == 0) then
        snold = DONE
        snnew = DONE
      else
        snold = sQuadraticSaturation(tp, bt, hold(n), this%satomega)
        snnew = sQuadraticSaturation(tp, bt, hnew(n), this%satomega)
      end if
      !
      ! -- storage coefficients
      sc1 = SsCapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
      rho1 = sc1 * tled
      !
      if(this%integratechanges /= 0) then
        ! -- Integration of storage changes (e.g. when using TVS):
        !    separate the old (start of time step) and new (end of time step)
        !    primary storage capacities
        sc1old = SsCapacity(this%istor_coef, tp, bt, this%dis%area(n), &
                            this%oldss(n))
        rho1old = sc1old * tled
      else
        ! -- No integration of storage changes: old and new values are
        !    identical => normal MF6 storage formulation
        rho1old = rho1
      end if
      !
      ! -- calculate specific storage terms
      call SsTerms(this%iconvert(n), this%iorig_ss, this%iconf_ss, tp, bt, &
                   rho1, rho1old, snnew, snold, hnew(n), hold(n), &
                   aterm, rhsterm)
      !
      ! -- add specific storage terms to amat and rhs
      amat(idxglo(idiag)) = amat(idxglo(idiag)) + aterm
      rhs(n) = rhs(n) + rhsterm
      !
      ! -- specific yield
      if (this%iconvert(n) /= 0) then
        rhsterm = DZERO
        !
        ! -- secondary storage coefficient
        sc2 = SyCapacity(this%dis%area(n), this%sy(n))
        rho2 = sc2 * tled
        !
        if(this%integratechanges /= 0) then
        ! -- Integration of storage changes (e.g. when using TVS):
        !    separate the old (start of time step) and new (end of time step)
        !    secondary storage capacities
          sc2old = SyCapacity(this%dis%area(n), this%oldsy(n))
          rho2old = sc2old * tled
        else
        ! -- No integration of storage changes: old and new values are
        !    identical => normal MF6 storage formulation
          rho2old = rho2
        end if
        !
        ! -- calculate specific storage terms
        call SyTerms(tp, bt, rho2, rho2old, snnew, snold, &
                     aterm, rhsterm)
!
        ! -- add specific yield terms to amat and rhs
        amat(idxglo(idiag)) = amat(idxglo(idiag)) + aterm
        rhs(n) = rhs(n) + rhsterm
      end if
    end do
    !
    ! -- return
    return
  end subroutine sto_fc

  !> @ brief Fill Newton-Raphson terms in A and right-hand side for the package
  !!
  !!  Fill the coefficient matrix and right-hand side with STO package
  !!  with Newton-Raphson terms.
  !!
  !<
  subroutine sto_fn(this, kiter, hold, hnew, njasln, amat, idxglo, rhs)
    ! -- modules
    use TdisModule, only: delt
    ! -- dummy variables
    class(GwfStoType) :: this                            !< GwfStoType object
    integer(I4B), intent(in) :: kiter                    !< outer iteration number
    real(DP), intent(in), dimension(:) :: hold           !< previous heads
    real(DP), intent(in), dimension(:) :: hnew           !< current heads
    integer(I4B), intent(in) :: njasln                   !< size of the A matrix for the solution
    real(DP), dimension(njasln), intent(inout) :: amat   !< A matrix
    integer(I4B), intent(in), dimension(:) :: idxglo     !< global index model to solution
    real(DP), intent(inout), dimension(:) :: rhs         !< right-hand side
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: tled
    real(DP) :: sc1
    real(DP) :: sc2
    real(DP) :: rho1
    real(DP) :: rho2
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: tthk
    real(DP) :: h
    real(DP) :: snnew
    real(DP) :: derv
    real(DP) :: rterm
    real(DP) :: drterm
    !
    ! -- test if steady-state stress period
    if (this%iss /= 0) return
    !
    ! -- set variables
    tled = DONE/delt
    !
    ! -- loop through and calculate storage contribution to hcof and rhs
    do n = 1, this%dis%nodes
      idiag = this%dis%con%ia(n)
      if (this%ibound(n) <= 0) cycle
      !
      ! -- aquifer elevations and thickness
      tp = this%dis%top(n)
      bt = this%dis%bot(n)
      tthk = tp - bt
      h = hnew(n)
      !
      ! -- aquifer saturation
      snnew = sQuadraticSaturation(tp, bt, h)
      !
      ! -- storage coefficients
      sc1 = SsCapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
      sc2 = SyCapacity(this%dis%area(n), this%sy(n))
      rho1 = sc1*tled
      rho2 = sc2*tled
      !
      ! -- calculate newton terms for specific storage
      !    and specific yield
      if (this%iconvert(n) /= 0) then
        !
        ! -- calculate saturation derivative
        derv = sQuadraticSaturationDerivative(tp, bt, h)
        !
        ! -- newton terms for specific storage
        if (this%iconf_ss == 0) then
          if (this%iorig_ss == 0) then
            drterm = -rho1 * derv * (h - bt) + rho1 * tthk * snnew * derv
          else
            drterm = -(rho1 * derv * h)
          end if
          amat(idxglo(idiag)) = amat(idxglo(idiag)) + drterm
          rhs(n) = rhs(n) + drterm*h
        end if
        !
        ! -- newton terms for specific yield
        !    only calculated if the current saturation
        !    is less than one
        if (snnew < DONE) then
          ! -- calculate newton terms for specific yield
          if (snnew > DZERO) then
            rterm = -rho2 * tthk * snnew
            drterm = -rho2 * tthk * derv
            amat(idxglo(idiag)) = amat(idxglo(idiag)) + drterm + rho2
            rhs(n) = rhs(n) - rterm + drterm * h + rho2 * bt
          end if
        end if
      end if
    end do
    !
    ! -- return
    return
  end subroutine sto_fn

  !> @ brief Calculate flows for package
  !!
  !!  Flow calculation for the STO package components. Components include
  !!  specific storage and specific yield storage.
  !!
  !<
  subroutine sto_cq(this, flowja, hnew, hold)
    ! -- modules
    use TdisModule, only: delt
    ! -- dummy variables
    class(GwfStoType) :: this                                    !< GwfStoType object
    real(DP), dimension(:), contiguous, intent(inout) :: flowja  !< connection flows
    real(DP), dimension(:), contiguous, intent(in) :: hnew       !< current head
    real(DP), dimension(:), contiguous, intent(in) :: hold       !< previous head
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: rate
    real(DP) :: tled
    real(DP) :: sc1
    real(DP) :: sc2
    real(DP) :: rho1
    real(DP) :: rho2
    real(DP) :: sc1old
    real(DP) :: sc2old
    real(DP) :: rho1old
    real(DP) :: rho2old
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: snold
    real(DP) :: snnew
    real(DP) :: aterm
    real(DP) :: rhsterm
    !
    ! -- initialize strg arrays
    do n = 1, this%dis%nodes
      this%strgss(n) = DZERO
      this%strgsy(n) = DZERO
    end do
    !
    ! -- Set strt to zero or calculate terms if not steady-state stress period
    if (this%iss == 0) then
      !
      ! -- set variables
      tled = DONE/delt
      !
      ! -- Calculate storage change
      do n = 1, this%dis%nodes
        if (this%ibound(n) <= 0) cycle
        ! -- aquifer elevations and thickness
        tp = this%dis%top(n)
        bt = this%dis%bot(n)
        !
        ! -- aquifer saturation
        if (this%iconvert(n) == 0) then
          snold = DONE
          snnew = DONE
        else
          snold = sQuadraticSaturation(tp, bt, hold(n), this%satomega)
          snnew = sQuadraticSaturation(tp, bt, hnew(n), this%satomega)
        end if
        !
        ! -- primary storage coefficient
        sc1 = SsCapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
        rho1 = sc1*tled
        !
        if(this%integratechanges /= 0) then
          ! -- Integration of storage changes (e.g. when using TVS):
          !    separate the old (start of time step) and new (end of time step)
          !    primary storage capacities
          sc1old = SsCapacity(this%istor_coef, tp, bt, this%dis%area(n), &
                              this%oldss(n))
          rho1old = sc1old * tled
        else
          ! -- No integration of storage changes: old and new values are
          !    identical => normal MF6 storage formulation
          rho1old = rho1
        end if
        !
        ! -- calculate specific storage terms and rate
        call SsTerms(this%iconvert(n), this%iorig_ss, this%iconf_ss, tp, bt, &
                      rho1, rho1old, snnew, snold, hnew(n), hold(n), &
                      aterm, rhsterm, rate)
        !
        ! -- save rate
        this%strgss(n) = rate
        !
        ! -- add storage term to flowja
        idiag = this%dis%con%ia(n)
        flowja(idiag) = flowja(idiag) + rate
        !
        ! -- specific yield
        rate = DZERO
        if (this%iconvert(n) /= 0) then
          !
          ! -- secondary storage coefficient
          sc2 = SyCapacity(this%dis%area(n), this%sy(n))
          rho2 = sc2 * tled
          !
          if(this%integratechanges /= 0) then
            ! -- Integration of storage changes (e.g. when using TVS):
            !    separate the old (start of time step) and new (end of time
            !    step) secondary storage capacities
            sc2old = SyCapacity(this%dis%area(n), this%oldsy(n))
            rho2old = sc2old * tled
          else
            ! -- No integration of storage changes: old and new values are
            !    identical => normal MF6 storage formulation
            rho2old = rho2
          end if
          !
          ! -- calculate specific yield storage terms and rate
          call SyTerms(tp, bt, rho2, rho2old, snnew, snold, &
                       aterm, rhsterm, rate)

        end if
        this%strgsy(n) = rate
        !
        ! -- add storage term to flowja
        idiag = this%dis%con%ia(n)
        flowja(idiag) = flowja(idiag) + rate
      end do
    end if
    !
    ! -- return
    return
  end subroutine sto_cq

  !> @ brief Model budget calculation for package
  !!
  !!  Budget calculation for the STO package components. Components include
  !!  specific storage and specific yield storage.
  !!
  !<
  subroutine sto_bd(this, isuppress_output, model_budget)
    ! -- modules
    use TdisModule, only: delt
    use BudgetModule, only: BudgetType, rate_accumulator
    ! -- dummy variables
    class(GwfStoType) :: this                          !< GwfStoType object
    integer(I4B), intent(in) :: isuppress_output       !< flag to suppress model output
    type(BudgetType), intent(inout) :: model_budget    !< model budget object
    ! -- local variables
    real(DP) :: rin
    real(DP) :: rout
    !
    ! -- Add confined storage rates to model budget
    call rate_accumulator(this%strgss, rin, rout)
    call model_budget%addentry(rin, rout, delt, budtxt(1), &
                               isuppress_output, '         STORAGE')
    !
    ! -- Add unconfined storage rates to model budget
    if (this%iusesy == 1) then
      call rate_accumulator(this%strgsy, rin, rout)
      call model_budget%addentry(rin, rout, delt, budtxt(2), &
                                 isuppress_output, '         STORAGE')
    end if
    !
    ! -- return
    return
  end subroutine sto_bd

  !> @ brief Save model flows for package
  !!
  !!  Save cell-by-cell budget terms for the STO package.
  !!
  !<
  subroutine sto_save_model_flows(this, icbcfl, icbcun)
    ! -- dummy variables
    class(GwfStoType) :: this           !< GwfStoType object
    integer(I4B), intent(in) :: icbcfl  !< flag to output budget data
    integer(I4B), intent(in) :: icbcun  !< cell-by-cell file unit number
    ! -- local variables
    integer(I4B) :: ibinun
    integer(I4B) :: iprint, nvaluesp, nwidthp
    character(len=1) :: cdatafmp = ' ', editdesc = ' '
    real(DP) :: dinact
    !
    ! -- Set unit number for binary output
    if (this%ipakcb < 0) then
      ibinun = icbcun
    elseif (this%ipakcb == 0) then
      ibinun = 0
    else
      ibinun = this%ipakcb
    end if
    if (icbcfl == 0) ibinun = 0
    !
    ! -- Record the storage rates if requested
    if (ibinun /= 0) then
      iprint = 0
      dinact = DZERO
      !
      ! -- storage(ss)
      call this%dis%record_array(this%strgss, this%iout, iprint, -ibinun, &
                                 budtxt(1), cdatafmp, nvaluesp, &
                                 nwidthp, editdesc, dinact)
      !
      ! -- storage(sy)
      if (this%iusesy == 1) then
        call this%dis%record_array(this%strgsy, this%iout, iprint, -ibinun, &
                                   budtxt(2), cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      end if
    end if
    !
    ! -- return
    return
  end subroutine sto_save_model_flows

  !> @ brief Deallocate package memory
  !!
  !!  Deallocate STO package scalars and arrays.
  !!
  !<
  subroutine sto_da(this)
    ! -- modules
    use MemoryManagerModule, only: mem_deallocate
    ! -- dummy variables
    class(GwfStoType) :: this  !< GwfStoType object
    !
    ! -- TVS
    if (this%intvs /= 0) then
      call this%tvs%da()
      deallocate(this%tvs)
    end if
    !
    ! -- Deallocate arrays if package is active
    if (this%inunit > 0) then
      call mem_deallocate(this%iconvert)
      call mem_deallocate(this%ss)
      call mem_deallocate(this%sy)
      call mem_deallocate(this%strgss)
      call mem_deallocate(this%strgsy)
      !
      ! -- deallocate TVS arrays
      if(associated(this%oldss)) then
        call mem_deallocate(this%oldss)
      end if
      if(associated(this%oldsy)) then
        call mem_deallocate(this%oldsy)
      end if
    end if
    !
    ! -- Deallocate scalars
    call mem_deallocate(this%istor_coef)
    call mem_deallocate(this%iconf_ss)
    call mem_deallocate(this%iorig_ss)
    call mem_deallocate(this%iusesy)
    call mem_deallocate(this%satomega)
    call mem_deallocate(this%integratechanges)
    call mem_deallocate(this%intvs)
    !
    ! -- deallocate parent
    call this%NumericalPackageType%da()
    !
    ! -- return
    return
  end subroutine sto_da

  !> @ brief Allocate scalars
  !!
  !! Allocate and initialize scalars for the STO package. The base numerical
  !! package allocate scalars method is also called.
  !!
  !<
  subroutine allocate_scalars(this)
    ! -- modules
    use MemoryManagerModule, only: mem_allocate, mem_setptr
    ! -- dummy variables
    class(GwfStoType) :: this  !< GwfStoType object
    !
    ! -- allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- allocate scalars
    call mem_allocate(this%istor_coef, 'ISTOR_COEF', this%memoryPath)
    call mem_allocate(this%iconf_ss, 'ICONF_SS', this%memoryPath)
    call mem_allocate(this%iorig_ss, 'IORIG_SS', this%memoryPath)
    call mem_allocate(this%iusesy, 'IUSESY', this%memoryPath)
    call mem_allocate(this%satomega, 'SATOMEGA', this%memoryPath)
    call mem_allocate(this%integratechanges, 'INTEGRATECHANGES', &
                      this%memoryPath)
    call mem_allocate(this%intvs, 'INTVS', this%memoryPath)
    !
    ! -- initialize scalars
    this%istor_coef = 0
    this%iconf_ss = 0
    this%iorig_ss = 0
    this%iusesy = 0
    this%satomega = DZERO
    this%integratechanges = 0
    this%intvs = 0
    !
    ! -- return
    return
  end subroutine allocate_scalars

  !> @ brief Allocate package arrays
  !!
  !!  Allocate and initialize STO package arrays.
  !!
  !<
  subroutine allocate_arrays(this, nodes)
    ! -- modules
    use MemoryManagerModule, only: mem_allocate
    ! -- dummy variables
    class(GwfStoType), target :: this !< GwfStoType object
    integer(I4B), intent(in) :: nodes !< active model nodes
    ! -- local variables
    integer(I4B) :: n
    !
    ! -- Allocate arrays
    call mem_allocate(this%iconvert, nodes, 'ICONVERT', this%memoryPath)
    call mem_allocate(this%ss, nodes, 'SS', this%memoryPath)
    call mem_allocate(this%sy, nodes, 'SY', this%memoryPath)
    call mem_allocate(this%strgss, nodes, 'STRGSS', this%memoryPath)
    call mem_allocate(this%strgsy, nodes, 'STRGSY', this%memoryPath)
    !
    ! -- Initialize arrays
    this%iss = 0
    do n = 1, nodes
      this%iconvert(n) = 1
      this%ss(n) = DZERO
      this%sy(n) = DZERO
      this%strgss(n) = DZERO
      this%strgsy(n) = DZERO
      if (this%integratechanges /= 0) then
        this%oldss(n) = DZERO
        if (this%iusesy /= 0) then
          this%oldsy(n) = DZERO
        end if
      end if
    end do
    !
    ! -- return
    return
  end subroutine allocate_arrays

  !> @ brief Read options for package
  !!
  !!  Read options block for STO package.
  !!
  !<
  subroutine read_options(this)
    ! -- modules
    ! -- dummy variables
    class(GwfStoType) :: this  !< GwfStoType object
    ! -- local variables
    character(len=LINELENGTH) :: keyword, fname
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    ! -- formats
    character(len=*), parameter :: fmtisvflow = &
                                   "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE "// &
                                   "WHENEVER ICBCFL IS NOT ZERO.')"
    character(len=*), parameter :: fmtflow = &
                                   "(4x, 'FLOWS WILL BE SAVED TO FILE: ', a, /4x, 'OPENED ON UNIT: ', I7)"
    character(len=*), parameter :: fmtorigss = &
      "(4X,'ORIGINAL_SPECIFIC_STORAGE OPTION:',/,                              &
      &1X,'The original specific storage formulation will be used')"
    character(len=*), parameter :: fmtstoc = &
      "(4X,'STORAGECOEFFICIENT OPTION:',/,                                     &
      &1X,'Read storage coefficient rather than specific storage')"
    character(len=*), parameter :: fmtconfss = &
      "(4X,'SS_CONFINED_ONLY OPTION:',/,                                       &
      &1X,'Specific storage changes only occur under confined conditions')"
    !
    ! -- get options block
    call this%parser%GetBlock('OPTIONS', isfound, ierr, &
                              supportOpenClose=.true., blockRequired=.false.)
    !
    ! -- parse options block if detected
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING STORAGE OPTIONS'
      do
        call this%parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('SAVE_FLOWS')
          this%ipakcb = -1
          write (this%iout, fmtisvflow)
        case ('STORAGECOEFFICIENT')
          this%istor_coef = 1
          write (this%iout, fmtstoc)
        case ('SS_CONFINED_ONLY')
          this%iconf_ss = 1
          this%iorig_ss = 0
          write (this%iout, fmtconfss)
        case ('TVS6')
          if (this%intvs /= 0) then
            errmsg = 'Multiple TVS6 keywords detected in OPTIONS block.' // &
                     ' Only one TVS6 entry allowed.'
            call store_error(errmsg, terminate=.TRUE.)
          end if
          call this%parser%GetStringCaps(keyword)
          if(trim(adjustl(keyword)) /= 'FILEIN') then
            errmsg = 'TVS6 keyword must be followed by "FILEIN" ' // &
                     'then by filename.'
            call store_error(errmsg, terminate=.TRUE.)
          endif
          call this%parser%GetString(fname)
          this%intvs = GetUnit()
          call openfile(this%intvs, this%iout, fname, 'TVS')
          call tvs_cr(this%tvs, this%name_model, this%intvs, this%iout)
        !
        ! -- right now these are options that are only available in the
        !    development version and are not included in the documentation.
        !    These options are only available when IDEVELOPMODE in
        !    constants module is set to 1
        case ('DEV_ORIGINAL_SPECIFIC_STORAGE')
          this%iorig_ss = 1
          write (this%iout, fmtorigss)
        case ('DEV_OLDSTORAGEFORMULATION')
          call this%parser%DevOpt()
          this%iconf_ss = 1
          this%iorig_ss = 0
          write (this%iout, fmtconfss)
        case default
          write (errmsg, '(4x,a,a)') 'Unknown STO option: ', &
            trim(keyword)
          call store_error(errmsg, terminate=.TRUE.)
        end select
      end do
      write (this%iout, '(1x,a)') 'END OF STORAGE OPTIONS'
    end if
    !
    ! -- set omega value used for saturation calculations
    if (this%inewton > 0) then
      this%satomega = DEM6
    end if
    !
    ! -- return
    return
  end subroutine read_options

  !> @ brief Read data for package
  !!
  !!  Read griddata block for STO package.
  !!
  !<
  subroutine read_data(this)
    ! -- modules
    ! -- dummy variables
    class(GwfStotype) :: this  !< GwfStoType object
    ! -- local variables
    character(len=LINELENGTH) :: keyword
    character(len=:), allocatable :: line
    character(len=LINELENGTH) :: cellstr
    integer(I4B) :: istart, istop, lloc, ierr
    logical :: isfound, endOfBlock
    logical :: readiconv
    logical :: readss
    logical :: readsy
    logical :: isconv
    character(len=24), dimension(4) :: aname
    integer(I4B) :: n
    ! -- formats
    !data
    data aname(1)/'                ICONVERT'/
    data aname(2)/'        SPECIFIC STORAGE'/
    data aname(3)/'          SPECIFIC YIELD'/
    data aname(4)/'     STORAGE COEFFICIENT'/
    !
    ! -- initialize
    isfound = .false.
    readiconv = .false.
    readss = .false.
    readsy = .false.
    isconv = .false.
    !
    ! -- get stodata block
    call this%parser%GetBlock('GRIDDATA', isfound, ierr)
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING GRIDDATA'
      do
        call this%parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call this%parser%GetStringCaps(keyword)
        call this%parser%GetRemainingLine(line)
        lloc = 1
        select case (keyword)
        case ('ICONVERT')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%iconvert, &
                                        aname(1))
          readiconv = .true.
        case ('SS')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%ss, &
                                        aname(2))
          readss = .true.
        case ('SY')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%sy, &
                                        aname(3))
          readsy = .true.
        case default
          write (errmsg, '(4x,a,a)') 'Unknown GRIDDATA tag: ', &
            trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END PROCESSING GRIDDATA'
    else
      write (errmsg, '(1x,a)') 'Required GRIDDATA block not found.'
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- Check for ICONVERT
    if (.not. readiconv) then
      write (errmsg, '(a, a, a)') 'Error in GRIDDATA block: ', &
        trim(adjustl(aname(1))), ' not found.'
      call store_error(errmsg)
    else
      isconv = .false.
      do n = 1, this%dis%nodes
        if (this%iconvert(n) /= 0) then
          isconv = .true.
          this%iusesy = 1
          exit
        end if
      end do
    end if
    !
    ! -- Check for SS
    if (.not. readss) then
      write (errmsg, '(a, a, a)') 'Error in GRIDDATA block: ', &
        trim(adjustl(aname(2))), ' not found.'
      call store_error(errmsg)
    end if
    !
    ! -- Check for SY
    if (.not. readsy .and. isconv) then
      write (errmsg, '(a, a, a)') 'Error in GRIDDATA block: ', &
        trim(adjustl(aname(3))), ' not found.'
      call store_error(errmsg)
    end if
    !
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- Check SS and SY for negative values
    do n = 1, this%dis%nodes
      if (this%ss(n) < DZERO) then
        call this%dis%noder_to_string(n, cellstr)
        write (errmsg, '(a,2(1x,a),1x,g0,1x,a)') &
          'Error in SS DATA: SS value in cell', trim(adjustl(cellstr)), &
          'is less than zero (', this%ss(n), ').'
        call store_error(errmsg)
      end if
      if (readsy) then
        if (this%sy(n) < DZERO) then
          call this%dis%noder_to_string(n, cellstr)
          write (errmsg, '(a,2(1x,a),1x,g0,1x,a)') &
            'Error in SY DATA: SY value in cell', trim(adjustl(cellstr)), &
            'is less than zero (', this%sy(n), ').'
          call store_error(errmsg)
        end if
      end if
    end do
    !
    ! -- return
    return
  end subroutine read_data

  !> @ brief Save old storage property values
  !!
  !!  Save ss and sy values from the previous time step for use with storage
  !!  integration when integratechanges is non-zero.
  !!
  !<
  subroutine save_old_ss_sy(this)
    ! -- modules
    use MemoryManagerModule, only: mem_allocate
    ! -- dummy variables
    class(GwfStoType) :: this  !< GwfStoType object
    ! -- local variables
    integer(I4B) :: n
    !
    ! -- Allocate TVS arrays if needed
    if(.not. associated(this%oldss)) then
      call mem_allocate(this%oldss, this%dis%nodes, 'OLDSS', this%memoryPath)
    end if
    if(this%iusesy == 1 .and. .not. associated(this%oldsy)) then
      call mem_allocate(this%oldsy, this%dis%nodes, 'OLDSY', this%memoryPath)
    end if
    !
    ! -- Save current specific storage
    do n = 1, this%dis%nodes
      this%oldss(n) = this%ss(n)
    end do
    !
    ! -- Save current specific yield, if used
    if(this%iusesy == 1) then
      do n = 1, this%dis%nodes
        this%oldsy(n) = this%sy(n)
      end do
    end if
    !
    ! -- Return
    return
  end subroutine save_old_ss_sy

end module GwfStoModule


text:
Subroutine color_palette(number_color,Color_matrix,color)
	use dflib
	integer::number_color,Color_matrix(4,8),color(1:256)
	integer(4),dimension(:),allocatable::Primary_color
	real(8)::red_step,green_step,blue_step,color_step
	allocate(primary_color(number_color))
	do i=1,number_color
		color(color_matrix(4,i))=rgbtointeger(color_matrix(1,i),color_matrix(2,i),color_matrix(3,i))
	end do
	i=1
	do while (i/=number_color)
		color_step=color_matrix(4,i+1)-color_matrix(4,i)
		red_step=(color_matrix(1,i+1)-color_matrix(1,i))/color_step
		green_step=(color_matrix(2,i+1)-color_matrix(2,i))/color_step
		blue_step=(color_matrix(3,i+1)-color_matrix(3,i))/color_step
		do j=color_matrix(4,i)+1,color_matrix(4,i+1)-1
			color(j)=rgbtointeger(int(color_matrix(1,i)+((j-color_matrix(4,i)+1)*red_step)),int(color_matrix(2,i)+((j-color_matrix(4,i)+1)*green_step)),int(color_matrix(3,i)+((j-color_matrix(4,i)+1)*blue_step)))
		end do
		i=i+1
	end do
end subroutine color_palette

max_stars_repo_path: Src/Color Palette.f90

max_stars_repo_name: Armanimani/FractalGenerator

max_issues_repo_path: Src/Color Palette.f90

max_issues_repo_name: Armanimani/FractalGenerator

max_forks_repo_path: Src/Color Palette.f90

max_forks_repo_name: Armanimani/FractalGenerator


text:
! THIS VERSION: GALAHAD 2.1 - 22/03/2007 AT 09:00 GMT.
! Updated 09/01/2004: format number corrected

!-*-*-*-*-*-*-*-  L A N C E L O T  -B-  FRNTL  M O D U L E  *-*-*-*-*-*-*-*-*

!  Nick Gould, for GALAHAD productions
!  Copyright reserved
!  February 2nd 1995

   MODULE LANCELOT_FRNTL_double

     USE GALAHAD_SMT_double
     USE GALAHAD_SILS_double
     USE GALAHAD_SCU_double, ONLY : SCU_matrix_type, SCU_data_type,            &
       SCU_info_type, SCU_restart_m_eq_0, SCU_solve, SCU_append
     USE LANCELOT_ASMBL_double
     USE LANCELOT_MDCHL_double

     IMPLICIT NONE

     PRIVATE
     PUBLIC :: FRNTL_get_search_direction

!  Set precision

     INTEGER, PARAMETER :: wp = KIND( 1.0D+0 )

!  Set other parameters

     INTEGER, PARAMETER :: liwmin = 1, lwmin = 1
     REAL ( KIND = wp ), PARAMETER :: zero = 0.0_wp
     REAL ( KIND = wp ), PARAMETER :: one = 1.0_wp
     REAL ( KIND = wp ), PARAMETER :: four = 4.0_wp
     REAL ( KIND = wp ), PARAMETER :: half = 0.5_wp
     REAL ( KIND = wp ), PARAMETER :: ten = 10.0_wp

   CONTAINS

!-*-*  L A N C E L O T  -B-   FRNTL_get_search_direction   M O D U L E  *-*-*

     SUBROUTINE FRNTL_get_search_direction(                                    &
        n , ng, nel   , ntotel, nnza  , maxsel, nvargp, buffer, INTVAR,        &
        IELVAR, nvrels, INTREP, IELING, ISTADG, ISTAEV, A     , ICNA  ,        &
        ISTADA, GUVALS, lnguvl, HUVALS, lnhuvl, ISTADH, GXEQX , GVALS2,        &
        GVALS3, IFREE , nfree , GMODEL, P , XT, BND   , fmodel, GSCALE,        &
        ESCALE, X0    , boundx, nobnds, dxsqr , radius, gstop , number,        &
        next  , modchl, RANGE , nsemib, ratio , iprint, error , out   ,        &
        status, alloc_status, bad_alloc , ITYPEE, DIAG  , OFFDIA, IVUSE ,      &
        RHS, RHS2, P2, ISTAGV, ISVGRP, lirnh, ljcnh, lh, LINK_col, POS_in_H,   &
        llink, lpos, IW_asmbl, W_ws, W_el, W_in, H_el, H_in,                   &
        matrix, SILS_data, SILS_cntl, SILS_infoa, SILS_infof, SILS_infos,      &
        SCU_matrix, SCU_data, SCU_info, SA, skipg, KNDOFG )

!  Use a factorization of H, the assembled Hessian matrix of a group
!  partially separable function, restricted to the IFREE variables,
!  to obtain the following search directions p:

!  If H is positive definite,
!     solve H * p = GMODEL,

!  If H is indefinite,
!     find a direction p such that

!      p( transpose ) * H * p < 0  AND  p( transpose ) * GMODEL > 0

!  If H is singular but H * p = GMODEL is consistent
!     solve H * p = GMODEL,

!  If H is singular but H * p = GMODEL is inconsistent
!     find a descent direction P such that

!     H * p = 0 and p( transpose ) * GMODEL > 0

!    Fortran 77 version:  M. Lescrenier

!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------

     INTEGER :: n , ng, maxsel, nvrels, ntotel, lnguvl, nvargp, status,        &
                lnhuvl, nfree , number, iprint, error, out, nnza, nel, buffer
     INTEGER, INTENT( IN ) :: nsemib
     INTEGER, INTENT( OUT ) :: alloc_status
     REAL ( KIND = wp ) fmodel, gstop, dxsqr, radius
     REAL ( KIND = wp ), INTENT( OUT ) :: ratio
     LOGICAL next, modchl, boundx, nobnds, skipg
     CHARACTER ( LEN = 24 ), INTENT( OUT ) :: bad_alloc
     INTEGER, DIMENSION( n       ) :: IFREE
     INTEGER, DIMENSION( ng  + 1 ) :: ISTADA, ISTADG
     INTEGER, DIMENSION( nel + 1 ) :: INTVAR, ISTAEV, ISTADH
     INTEGER, DIMENSION( nvrels  ) :: IELVAR
     INTEGER, DIMENSION( ntotel  ) :: IELING
     INTEGER, DIMENSION( nnza    ) :: ICNA
     INTEGER, INTENT( IN ), DIMENSION( : ) :: ITYPEE
     REAL ( KIND = wp ), DIMENSION( n ) :: GMODEL, P, X0, XT
     REAL ( KIND = wp ), DIMENSION( ng ) :: GVALS2, GVALS3, GSCALE
     REAL ( KIND = wp ), DIMENSION( : , : ) :: BND
     REAL ( KIND = wp ), DIMENSION( ntotel ) :: ESCALE
     REAL ( KIND = wp ), DIMENSION( nnza ) :: A
     REAL ( KIND = wp ), DIMENSION( lnguvl ) :: GUVALS
     REAL ( KIND = wp ), DIMENSION( lnhuvl ) :: HUVALS
     LOGICAL, DIMENSION( nel ) :: INTREP
     LOGICAL, DIMENSION( ng  ) :: GXEQX
     TYPE ( SCU_matrix_type ), INTENT( INOUT ) :: SCU_matrix
     TYPE ( SCU_data_type ), INTENT( INOUT ) :: SCU_data
     TYPE ( SCU_info_type ), INTENT( INOUT ) :: SCU_info
     TYPE ( SMT_type ), INTENT( INOUT ) :: matrix
     TYPE ( SILS_factors ), INTENT( INOUT ) :: SILS_data
     TYPE ( SILS_control ), INTENT( INOUT ) :: SILS_cntl
     TYPE ( SILS_ainfo ), INTENT( INOUT ) :: SILS_infoa
     TYPE ( SILS_finfo ), INTENT( INOUT ) :: SILS_infof
     TYPE ( SILS_sinfo ), INTENT( INOUT ) :: SILS_infos
     TYPE ( ASMBL_save_type ), INTENT( INOUT ) :: SA

     INTEGER, INTENT( IN ), OPTIONAL, DIMENSION( ng ) :: KNDOFG

!---------------------------------------------------------------
!   D u m m y   A r g u m e n t s   f o r   W o r k s p a c e
!--------------------------------------------------------------

     INTEGER, ALLOCATABLE, DIMENSION( : ) :: IVUSE
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: RHS, RHS2, P2, DIAG
     REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : , : ) :: OFFDIA

     INTEGER, INTENT( IN ), DIMENSION( : ) :: ISVGRP
     INTEGER, INTENT( IN ), DIMENSION( : ) :: ISTAGV

     INTEGER, INTENT( INOUT ) :: lirnh, ljcnh, lh, llink, lpos
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: LINK_col
     INTEGER, ALLOCATABLE, DIMENSION( : ) :: POS_in_H

     INTEGER, INTENT( OUT ), DIMENSION( : ) :: IW_asmbl
     REAL ( KIND = wp ), INTENT( OUT ), DIMENSION( : ) :: W_ws
     REAL ( KIND = wp ), INTENT( OUT ), DIMENSION( : ) :: W_el
     REAL ( KIND = wp ), INTENT( OUT ), DIMENSION( : ) :: W_in
     REAL ( KIND = wp ), INTENT( OUT ), DIMENSION( : ) :: H_el
     REAL ( KIND = wp ), INTENT( OUT ), DIMENSION( : ) :: H_in

!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------

     INTERFACE
       SUBROUTINE RANGE( ielemn, transp, W1, W2, nelvar, ninvar, ieltyp,       &
                         lw1, lw2 )
       INTEGER, INTENT( IN ) :: ielemn, nelvar, ninvar, ieltyp, lw1, lw2
       LOGICAL, INTENT( IN ) :: transp
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( IN ), DIMENSION ( lw1 ) :: W1
       REAL ( KIND = KIND( 1.0D+0 ) ), INTENT( OUT ), DIMENSION ( lw2 ) :: W2
       END SUBROUTINE RANGE
     END INTERFACE

!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------

     INTEGER :: i, j, nnzh, max_sc
     INTEGER :: nfixed, scu_status, nsemiw
     REAL ( KIND = wp ) :: dotprd, pertur, alpha , sgrad , sths
     REAL ( KIND = wp ) :: pi, hp, ptp   , dxtp  , steptr, big
     REAL ( KIND = wp ) :: onepep, amodel, atemp , gnrmsq, epsmch
     LOGICAL :: pronel, prnter, consis, negate, reallocate

!-----------------------------------------------
!   A l l o c a t a b l e   A r r a y s
!-----------------------------------------------

     big = HUGE( one ) ;  epsmch = EPSILON( one )
     pronel = iprint == 2 ; prnter = iprint >= 5
     IF ( pronel ) WRITE( out, 2120 ) gstop
     onepep = one + ten * epsmch
     nfixed = 0
     max_sc = SCU_matrix%m_max

!  Assemble the Hessian restricted to the IFREE variables

     nsemiw = nsemib
     CALL ASMBL_assemble_hessian(                                              &
                n , ng, nel   , ntotel, nvrels, nnza  , maxsel,                &
                nvargp, nfree , IFREE , ISTADH, ICNA  ,                        &
                ISTADA, INTVAR, IELVAR, IELING, ISTADG, ISTAEV,                &
                ISTAGV, ISVGRP, A     , GUVALS, lnguvl, HUVALS,                &
                lnhuvl, GVALS2, GVALS3, GSCALE, ESCALE, GXEQX ,                &
                ITYPEE, INTREP, RANGE , iprint, error, out, buffer, .FALSE.,   &
                .TRUE., .FALSE., nsemiw, status, alloc_status, bad_alloc, SA,  &
                lirnh, ljcnh, lh, matrix%row, matrix%col, matrix%val,          &
                LINK_col, POS_in_H,                                            &
                llink, lpos, IW_asmbl, W_ws, W_el, W_in, H_el, H_in,           &
                skipg, nnzh = nnzh, KNDOFG = KNDOFG )

     IF ( status /= 0 ) RETURN

!  Choose initial values for the control parameters

     IF ( iprint >= 1000 ) THEN
       SILS_cntl%lp = 6 ; SILS_cntl%mp = 6
       SILS_cntl%wp = 6 ; SILS_cntl%sp = 6
       SILS_cntl%ldiag = 1
     ELSE
       SILS_cntl%lp = 0 ; SILS_cntl%mp = 0
       SILS_cntl%wp = 0 ; SILS_cntl%sp = 0
       SILS_cntl%ldiag = 0
     END IF

!  Choose pivots for Gaussian elimination

     matrix%n = nfree
     matrix%ne = nnzh
     CALL SILS_analyse( matrix, SILS_data, SILS_cntl, SILS_infoa )

!  Allocate arrays for MA27B or MCFA, MA27C, MISC_negcur or MISC_dscdir

     IF ( modchl ) THEN

!  Allocate further workspace

       reallocate = .TRUE.
       IF ( ALLOCATED( DIAG ) ) THEN
          IF ( SIZE( DIAG ) < nfree ) THEN ; DEALLOCATE( DIAG )
          ELSE ; reallocate = .FALSE.
          END IF
       END IF
       IF ( reallocate ) THEN
          ALLOCATE( DIAG( nfree ), STAT = alloc_status )
          IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'DIAG' ; GO TO 980
          END IF
       END IF
     END IF

!  Factorize the restricted Hessian

     IF ( modchl ) SILS_cntl%pivoting = 4
     CALL SILS_factorize( matrix, SILS_data, SILS_cntl, SILS_infof )
     pertur = SILS_infof%maxchange

!  Test that the factorization succeeded

     IF ( SILS_infof%flag < 0 ) RETURN

!  Record the relative fill-in

     IF ( nnzh > 0 ) THEN
       ratio = DBLE( FLOAT( SILS_infof%nrlbdu ) ) / DBLE( FLOAT( nnzh ) )
     ELSE
       ratio = one
     END IF

!  Allocate further workspace

     reallocate = .TRUE.
     IF ( ALLOCATED( RHS ) ) THEN
       IF ( SIZE( RHS ) < nfree ) THEN
         DEALLOCATE( RHS ) ; ELSE ; reallocate = .FALSE. ; END IF
     END IF
     IF ( reallocate ) THEN
       ALLOCATE( RHS( nfree ), STAT = alloc_status )
       IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'RHS' ; GO TO 980 ; END IF
     END IF

     SCU_matrix%n = nfree
     SCU_matrix%m = nfixed
     CALL SCU_restart_m_eq_0( SCU_data, SCU_info )

     reallocate = .TRUE.
     IF ( ALLOCATED( RHS2 ) ) THEN
        IF ( SIZE( RHS2 ) < nfree + max_sc ) THEN ; DEALLOCATE( RHS2 )
        ELSE ; reallocate = .FALSE.
        END IF
     END IF
     IF ( reallocate ) THEN
        ALLOCATE( RHS2( nfree + max_sc ), STAT = alloc_status )
        IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'RHS2' ; GO TO 980
        END IF
     END IF

     reallocate = .TRUE.
     IF ( ALLOCATED( P2 ) ) THEN
        IF ( SIZE( P2 ) <  nfree + max_sc ) THEN ; DEALLOCATE( P2 )
        ELSE ; reallocate = .FALSE.
        END IF
     END IF
     IF ( reallocate ) THEN
        ALLOCATE( P2(  nfree + max_sc ), STAT = alloc_status )
        IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'P2' ; GO TO 980
        END IF
     END IF

     IF ( .NOT. modchl ) THEN
       reallocate = .TRUE.
       IF ( ALLOCATED( OFFDIA ) ) THEN
         IF ( SIZE( OFFDIA, 1 ) /= 2 .OR.                                &
              SIZE( OFFDIA, 2 ) < nfree ) THEN ; DEALLOCATE( OFFDIA )
         ELSE ; reallocate = .FALSE. ; END IF
       END IF
       IF ( reallocate ) THEN
         ALLOCATE( OFFDIA( 2, nfree ), STAT = alloc_status )
         IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'OFFDIA'; GO TO 980
         END IF
       END IF

       reallocate = .TRUE.
       IF ( ALLOCATED( IVUSE ) ) THEN
          IF ( SIZE( IVUSE ) < nfree ) THEN
            DEALLOCATE( IVUSE ) ; ELSE ; reallocate = .FALSE.
          END IF
       END IF
       IF ( reallocate ) THEN
         ALLOCATE( IVUSE( nfree ), STAT = alloc_status )
         IF ( alloc_status /= 0 ) THEN ; bad_alloc = 'IVUSE' ; GO TO 980 ; END IF
       END IF

     END IF

     IF ( modchl ) THEN
       status = 1
     ELSE
       status = MDCHL_block_type( nfree, SILS_infof%rank, SILS_data,           &
                                  IVUSE( : nfree ), OFFDIA )
     END IF

!  Define the model of the gradient restricted to the free variables in RHS

     RHS( : nfree ) = - GMODEL( IFREE( : nfree ) )
     IF ( status == 1 ) THEN

!  The restricted Hessian is positive definite -
!  solve the linear system H * s = - GMODEL. s is output in GMODEL.

       CALL SILS_solve( matrix, SILS_data, RHS, SILS_cntl, SILS_infos )
       negate = .FALSE.
     END IF
     IF ( status == 2 ) THEN

!  The restricted Hessian is not positive definite. Compute a direction of
!  negative curvature

       CALL MDCHL_get_donc( nfree, SILS_infof%rank, SILS_data,                 &
                            IVUSE( : nfree ), OFFDIA, number, next, RHS, sths, &
                            SILS_cntl, alloc_status )
       IF ( alloc_status /= 0 ) THEN
         bad_alloc = 'SILS_data%w'; GO TO 980 ; END IF

!  Check if the dot product of this direction times GMODEL is > 0

!      dotprd = - DOT_PRODUCT( GMODEL( IFREE( : nfree ) ),                     &
!                              RHS( : nfree ) )
       dotprd = zero
       DO i = 1, nfree
          dotprd = dotprd - GMODEL( IFREE( i ) ) * RHS( i )
       END DO
       negate = dotprd <= zero
     END IF
     IF ( status == 3 ) THEN
       CALL MDCHL_get_singular_direction(                                      &
                 nfree, SILS_infof%rank, SILS_data, IVUSE( : nfree ), OFFDIA,  &
                 RHS( : nfree ), consis, SILS_cntl, alloc_status )
       IF ( alloc_status /= 0 ) THEN
         bad_alloc = 'SILS_data%w'; GO TO 980 ; END IF
       negate = .FALSE.
       IF ( .NOT. consis ) THEN
         status = 4
         sths = zero
       END IF
     END IF

!  Construct the solution P from RHS

     IF ( negate ) THEN
       DO j = 1, nfree
         P( IFREE( j ) ) = - RHS( j )
       END DO
     ELSE
       DO j = 1, nfree
         P( IFREE( j ) ) = RHS( j )
       END DO
     END IF

!  Calculate the slope of the quadratic model from the Cauchy point

!    sgrad = DOT_PRODUCT( P     ( IFREE( : nfree ) ),                      &
!                         GMODEL( IFREE( : nfree ) ) )
     sgrad = zero
     DO i = 1, nfree
       sgrad = sgrad + P( IFREE( i ) ) * GMODEL( IFREE( i ) )
     END DO

!  Obtain the curvature

     IF ( status == 1 .OR. status == 3 ) sths = - sgrad

!  Modify the curvature if the Hessian has been perturbed

     IF ( pertur > zero ) THEN
        status = 5
!       sths = sths - DOT_PRODUCT( DIAG( : nfree ),                            &
!                                  RHS ( : nfree ) ** 2 )
        DO i = 1, nfree
           sths = sths - DIAG( i ) * RHS ( i ) ** 2
        END DO
     END IF
 100 CONTINUE

!  Compute the steplength to the minimizer of the model along the
!  current search direction

     IF ( status == 1 .OR. status == 3 .OR. status == 5 ) THEN
       IF ( status == 5 ) THEN
         IF ( sths > zero ) THEN
           alpha = - sgrad / sths
         ELSE
           alpha = big
         END IF
       ELSE
         alpha = one
       END IF
     ELSE
       alpha = big
     END IF

!  If required, print details of the current step

     IF ( iprint >= 20 ) THEN
       WRITE( out, 2060 )
       DO j = 1, nfree
         i = IFREE( j )
         IF ( i > 0 ) WRITE( out, 2070 )                                      &
           i, XT( i ), GMODEL( i ), P( i ), BND( i, 1 ), BND( i, 2 )
       END DO
     END IF

!  Find the largest feasible step in the direction P from XT

     amodel = alpha
     IF ( .NOT. nobnds ) THEN
       DO j = 1, nfree
         i = IFREE( j )
         IF ( i > 0 ) THEN
           IF ( ABS( P( i ) ) > epsmch ) THEN
             IF ( P( i ) > zero ) THEN
               alpha = MIN( alpha, ( BND( i, 2 ) - XT( i ) ) / P( i ) )
             ELSE
               alpha = MIN( alpha, ( BND( i, 1 ) - XT( i ) ) / P( i ) )
             END IF
           END IF
         END IF
       END DO
     END IF

     IF ( boundx ) THEN
       dxtp = zero ; ptp = zero
!DIR$ IVDEP
       DO j = 1, nfree
         i = IFREE( j )
         IF ( i > 0 ) THEN
           pi = P( i )
           dxtp = dxtp + pi * ( XT( i ) - X0( i ) )
           ptp = ptp + pi * pi
         END IF
       END DO

!  Now compute the distance to the spherical boundary, steptr, and
!  find the smaller of this and the distance to the boundary of the
!  feasible box, STEPMX

       steptr = ( SQRT( dxtp * dxtp - ptp * ( dxsqr - radius ** 2 ) ) -  &
                  dxtp ) / ptp
       alpha = MIN( steptr, alpha )
     END IF

!  Update the model function value

     fmodel = fmodel + alpha * ( sgrad + alpha * half * sths )

!  If required, print the model value, slope and the step taken

     IF ( prnter ) WRITE( out, 2040 ) amodel, alpha, fmodel, sgrad, sths

!  Check to see if the boundary is encountered

     gnrmsq = zero
     IF ( alpha < amodel ) THEN

!  A free variable has encountered a bound. Make a second pass to determine
!  which variable and compute the size of the model gradient at the new point

       DO j = 1, nfree
         i = IFREE( j )
         IF ( i > 0 ) THEN
           IF ( ABS( P( i ) ) > epsmch ) THEN
             atemp = big
             IF ( .NOT. nobnds ) THEN
               IF ( P( i ) > zero ) THEN
                 atemp = ( BND( i, 2 ) - XT( i ) ) / P( i )
               ELSE
                 atemp = ( BND( i, 1 ) - XT( i ) ) / P( i )
               END IF
             END IF

!  Variable I encounters its bound

             IF ( atemp <= alpha * onepep ) THEN
               IF ( prnter ) WRITE( out, 2090 ) i
               IFREE( j ) = - i
             ELSE

!  Update the gradient

                IF ( status == 1 .OR. status >= 4 ) hp = - GMODEL( i )
                IF ( status == 4 ) hp = zero
                IF ( status == 5 ) hp = hp - DIAG( j ) * P( i )
                GMODEL( i ) = GMODEL( i ) + alpha * hp
                gnrmsq = gnrmsq + GMODEL( i ) ** 2
              END IF
             XT( i ) = XT( i ) + alpha * P( i )
           END IF
         END IF
       END DO
     ELSE

!  Step to the new point

       DO j = 1, nfree
         i = IFREE( j )
         IF ( i > 0 ) THEN
           IF ( ABS( P( i ) ) > epsmch ) XT( i ) = XT( i ) + alpha * P( i )
         END IF
       END DO
     END IF
     IF ( boundx ) THEN
       IF ( alpha == steptr ) THEN
         IF ( prnter ) WRITE( out, 2050 ) gnrmsq, gstop ; GO TO 500 ; END IF
     END IF

!  If the model gradient is sufficiently small, exit

     IF ( prnter ) WRITE( out, 2050 ) gnrmsq, gstop
     IF ( pronel ) WRITE( out, 2130 ) nfixed, fmodel, gnrmsq
     IF ( gnrmsq <= gstop ) GO TO 500

!  Continue the minimization in the restricted subspace

     IF ( .NOT. modchl ) GO TO 500

!  Exit if there are no bounds

     IF ( alpha < amodel .AND. nobnds ) GO TO 500

!  Set up further workspace addresses to partition the unused portions of W
!  and IW. Calculate the largest number of variables which can be fixed

!  Determine which variables are to be fixed

     DO j = 1, nfree
       i = IFREE( j )
       IF ( i < 0 ) THEN

!  If more than max_sc variables have been fixed, return

         IF ( nfixed >= max_sc ) GO TO 500

!  Update the factorization of the Schur complement to allow for the removal of
!  the J-th row and column of the original Hessian - this removal is effected
!  by appending the J-th row and column of the identity matrix to the Hessian

         SCU_matrix%BD_val( nfixed + 1 ) = one
         SCU_matrix%BD_row( nfixed + 1 ) = j
         SCU_matrix%BD_col_start( nfixed + 2 ) = nfixed + 2
         scu_status = 1
 210     CONTINUE

!  Call SCU_append to update the Schur-complement

         CALL SCU_append( SCU_matrix, SCU_data, RHS, scu_status, SCU_info )

!  SCU_append requires additional information. Compute the solution to
!  the equation H * S = RHS, returning the solution S in RHS

         IF ( scu_status > 0 ) THEN
            CALL SILS_solve( matrix, SILS_data, RHS, SILS_cntl, SILS_infos )
            GO TO 210
         END IF
         IF ( scu_status < 0 ) THEN
            IF ( prnter ) WRITE( out, 2100 ) scu_status
            GO TO 500
         END IF
         nfixed = nfixed + 1
       END IF
     END DO

!  Define the model of the gradient restricted to the free variables in RHS

     DO j = 1, nfree
       i = IFREE( j )
       IF ( i > 0 ) THEN
         RHS2( j ) = - GMODEL( i )
       ELSE
         RHS2( j ) = zero
       END IF
     END DO
     RHS2( nfree + 1 : nfree + nfixed ) = zero

!  Solve the new linear system H * P2 = RHS2

     scu_status = 1
 310 CONTINUE

!  Call SCU_solve to solve the system

     CALL SCU_solve( SCU_matrix, SCU_data, RHS2, P2, RHS, scu_status )
     IF ( scu_status > 0 ) THEN

!  SCU_solve requires additional information. Compute the solution to
!  the equation H * S = RHS, returning the solution S in RHS

       CALL SILS_solve( matrix, SILS_data, RHS, SILS_cntl, SILS_infos )
       GO TO 310
     END IF

!  Store the corrections to the free variables

     DO j = 1, nfree
       i = IFREE( j )
       IF ( i < 0 ) IFREE( j ) = 0
       IF ( i > 0 ) P( i ) = P2( j )
     END DO

!  Calculate the slope of the quadratic model from the Cauchy point

     sgrad = zero
     DO j = 1, nfree
       i = IFREE( j )
       IF ( i > 0 ) sgrad = sgrad + P( i ) * GMODEL( i )
     END DO

!  Obtain the curvature

     IF ( status == 1 .OR. status == 3 .OR. status == 5 ) sths = - sgrad

!  Modify the curvature if the Hessian has been perturbed

     IF ( pertur > zero ) THEN
       DO j = 1, nfree
         IF ( IFREE( j ) > 0 ) sths = sths - DIAG( j ) * P2( j ) ** 2
       END DO
     END IF
     GO TO 100

!  Successful return

 500 CONTINUE
     RETURN

!  Unsuccessful returns

 980 CONTINUE
     WRITE( error, 2990 ) bad_alloc, alloc_status
     RETURN

!  Non-executable statements

 2040  FORMAT( /,' Model step and actual step = ', 2ES12.4, /,                 &
                 ' fmodel, sgrad, sths = ', 3ES12.4 )
 2050  FORMAT( /,' Model gradient ** 2 ', ES12.4, ' gstop ', ES12.4 )
 2060  FORMAT( /,'    i      XT          G           P           BL',          &
                '          BU' )
 2070  FORMAT( I5, 5ES12.4 )
 2090  FORMAT( ' Variable ',I6,' has encountered a bound in FRNTL ' )
 2100  FORMAT( ' ** Message from -FRNTL-', /,                                  &
               '    value of status after SCU_append = ', I3 )
 2120  FORMAT( /,'    ** FRNTL entered ** nfixed     MODEL   ',                &
                 '   Gradient gstop = ', ES12.4 )
 2130  FORMAT( 24X, I7, 3ES12.4 )
 2990  FORMAT( ' ** Message from -FRNTL_get_search_directio-', /,              &
               ' Allocation error (status = ', I6, ') for ', A24 )

!  End of subroutine FRNTL_get_search_direction

     END SUBROUTINE FRNTL_get_search_direction

!  End of module LANCELOT_FRNTL

   END MODULE LANCELOT_FRNTL_double

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
! { dg-do compile }
! PR42999  bogus error: Parameter 'i' at (1) has not been declared
! or is a variable, which does not reduce to a constant expression
 TYPE DD
  INTEGER :: I
 END TYPE DD
 TYPE(DD) :: X(2)=(/(DD(I),I=1,2)/)
 END

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      PROGRAM xmglin
C     driver for routine mglin
      INTEGER JMAX,NSTEP
      PARAMETER(JMAX=33,NSTEP=4)
      INTEGER i,j,midl
      DOUBLE PRECISION f(JMAX,JMAX),u(JMAX,JMAX)
      do 12 i=1,JMAX
        do 11 j=1,JMAX
          u(i,j)=0.d0
11      continue
12    continue
      midl=JMAX/2+1
      u(midl,midl)=2.d0
      call mglin(u,JMAX,2)
      write(*,'(1x,a)') 'MGLIN Solution:'
      do 13 i=1,JMAX,NSTEP
        write(*,'(1x,9f8.4)') (u(i,j),j=1,JMAX,NSTEP)
13    continue
      write(*,'(/1x,a)') 'Test that solution satisfies Difference Eqns:'
      do 15 i=NSTEP+1,JMAX-1,NSTEP
        do 14 j=NSTEP+1,JMAX-1,NSTEP
          f(i,j)=u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-4.d0*u(i,j)
14      continue
        write(*,'(7x,7f8.4)')
     *       (f(i,j)*(JMAX-1)*(JMAX-1),j=NSTEP+1,JMAX-1,NSTEP)
15    continue
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xmglin.for


text:
MODULE SUB_Swap
!!#### PURPOSE
!! Swaps the two values of the two arguments.

!!#### USAGE
!! The code snippet
!
!    x = a
!    y = b
!    CALL Swap( x , y )
!
!! yields <x=b> and <y=a>.

!!#### EXTERNAL KINDS
USE KND_IntrinsicTypes,ONLY: KIND_Rsp,KIND_Rdp,& !!((01-A-KND_IntrinsicTypes.f90))
              KIND_Csp,KIND_Cdp,&
              KIND_I1,KIND_I2,KIND_I4,KIND_I8,&
              KIND_L1,KIND_L2,KIND_L4,&
              KIND_S

!!#### FORTRAN STANDARDS
USE ISO_varying_string                           !!((03-A-ISO_varying_string.f90))

!!#### DEFAULT IMPLICIT
IMPLICIT NONE

!!#### DEFAULT ACCESS
PRIVATE

!!#### SUBROUTINE OVERLOADING
INTERFACE Swap
 MODULE PROCEDURE Swap_L1
 MODULE PROCEDURE Swap_L2
 MODULE PROCEDURE Swap_L4
 MODULE PROCEDURE Swap_I1
 MODULE PROCEDURE Swap_I2
 MODULE PROCEDURE Swap_I4
 MODULE PROCEDURE Swap_I8
 MODULE PROCEDURE Swap_Rsp
 MODULE PROCEDURE Swap_Rdp
 MODULE PROCEDURE Swap_Csp
 MODULE PROCEDURE Swap_Cdp
 MODULE PROCEDURE Swap_VS
 MODULE PROCEDURE Swap_S
END INTERFACE

!!#### PUBLIC ACCESS LIST
PUBLIC :: Swap

CONTAINS

SUBROUTINE Swap_L1( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L1     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_L.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_L2( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L2     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_L.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_L4( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L4     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_L.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_I1( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I1     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_I.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_I2( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I2     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_I.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_I4( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I4     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_I.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_I8( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I8     !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_I.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_Rsp( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rsp    !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_R.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_Rdp( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rdp    !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_R.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_Csp( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_C=>KIND_Csp    !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_C.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_Cdp( x , y )
USE KND_IntrinsicTypes,ONLY: KIND_C=>KIND_Cdp    !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "04-A-SUB_Swap_C.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_S( x , y )
INCLUDE "04-A-SUB_Swap_S.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

SUBROUTINE Swap_VS( x , y )
INCLUDE "04-A-SUB_Swap_VS.f90.hdr"
!!--begin--
INCLUDE "04-A-SUB_Swap.f90.bdy"
!!--end--
END SUBROUTINE

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
c ccplib.f
C     ==================================
      REAL FUNCTION STHLSQ(IH, IK, IL)
C     ==================================
C
      IMPLICIT NONE
C
      EXTERNAL STHLSQ1
C
C     PARAMETERS
      INTEGER IH, IK, IL
C
C     reso RETURN value
      REAL RESO
C
      CALL STHLSQ1(RESO, IH, IK, IL)
      STHLSQ=RESO
      RETURN
      END


text:
!
!	Get item scores.  Use numcatobs to find dimension and defaults.
!	numcatobs is the number of observed categories per item
function getscores(numcatobs)
implicit none
integer,intent(in)::numcatobs(:)
real(kind=8)::getscores(sum(numcatobs))
!	Error flag is io and category counter is cat.  Item count is item.  Position is placement of
!	entry for item score for a specific item and category.
interface
    subroutine readscores(scores)
        implicit none
        real(kind=8),intent(inout)::scores(:)
    end subroutine readscores
end interface
integer::cat,io,item,position
!	copy of getscores.
real(kind=8),allocatable::scores(:)

!	Default setting.
position=1
allocate(scores(sum(numcatobs)),stat=io)
if(io/=0) stop 'Allocation failed for item scores.'
do item=1,size(numcatobs)
	do cat=1,numcatobs(item)
		scores(position)=cat-1
		position=position+1
	end do
end do
call readscores(scores)

getscores=scores
return
end function getscores
subroutine readscores(scores)
implicit none
integer::io
real(kind=8),intent(inout)::scores(:)
namelist/itemscores/scores
read(*,nml=itemscores,iostat=io)
if(io/=0) stop 'Item scores not read successfully.'
return
end subroutine readscores

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


text:
      SUBROUTINE CREST

C   create a structure network

      REAL        XA(2), YA(2)
      DATA        XA/0.,0./, YA/0.5,0.5/

      CALL POPST (100)
      CALL PPL (2, XA, YA)
      CALL PEXST (101)
      CALL PCLST
      CALL POPST (101)
      CALL PSCHH (0.1)
      CALL PTX (0.55, 0.55, 'TEXT')
      CALL PCLST

      END