text:
! =============================================================================
!               This module contains all ellipse operations.
! =============================================================================
module parcel_ellipse
    use constants, only : pi   &
                        , fpi  &
                        , zero &
                        , two  &
                        , f12  &
                        , f14
    implicit none

    contains

        ! Obtain the largest eigenvalue (i.e. the semi-major axis squared [a**2])
        ! @param[in] B11 is the B matrix (1, 1) element
        ! @param[in] B12 is the B matrix (1, 2) element
        ! @param[in] B22 is the B matrix (2, 2) element
        ! @returns the largest eigenvalue
        function get_eigenvalue(B11, B12, B22) result(a2)
            double precision, intent(in) :: B11
            double precision, intent(in) :: B12
            double precision, intent(in) :: B22
            double precision             :: a2

            a2 = f12 * (B11 + B22) + dsqrt(f14 * (B11 - B22) ** 2 + B12 ** 2)
        end function get_eigenvalue

        ! Obtain the eigenvector of the largest eigenvalue
        ! @param[in] a2 is the largest eigenvalue
        ! @param[in] B11 is the B matrix (1, 1) element
        ! @param[in] B12 is the B matrix (1, 2) element
        ! @param[in] B22 is the B matrix (2, 2) element
        ! @returns the eigenvector
        function get_eigenvector(a2, B11, B12, B22) result(evec)
            double precision, intent(in) :: a2
            double precision, intent(in) :: B11
            double precision, intent(in) :: B12
            double precision, intent(in) :: B22
            double precision             :: evec(2)

            evec(1) = a2 - B22
            evec(2) = B12

            if (dabs(evec(1)) + dabs(evec(2)) == zero) then
                if (B11 > B22) then
                    evec(1) = evec(1) + epsilon(evec(1))
                else
                    evec(2) = evec(2) + epsilon(evec(2))
                endif
            endif

            evec = evec / norm2(evec)

        end function get_eigenvector

        ! used in unit tests only
        function get_angle(B11, B12, volume) result(angle)
            double precision, intent(in) :: B11
            double precision, intent(in) :: B12
            double precision, intent(in) :: volume
            double precision             :: B22
            double precision             :: a2
            double precision             :: evec(2)
            double precision             :: angle

            B22 = get_B22(B11, B12, volume)

            a2 = get_eigenvalue(B11, B12, B22)

            evec = get_eigenvector(a2, B11, B12, B22)

            angle = datan2(evec(2), evec(1))

        end function get_angle

        ! Obtain the B22 matrix element
        ! @param[in] B11 is the B matrix (1, 1) element
        ! @param[in] B12 is the B matrix (1, 2) element
        ! @param[in] volume of the parcel
        ! @returns B22
        elemental function get_B22(B11, B12, volume) result(B22)
            double precision, intent(in) :: B11
            double precision, intent(in) :: B12
            double precision, intent(in) :: volume
            double precision             :: B22
            double precision             :: ab

            ab = get_ab(volume)

            B22 = (ab ** 2 + B12 ** 2) / B11

        end function get_B22

        ! Obtain the product of the semi-minor and semi-major axis.
        ! @param[in] volume of the parcel
        ! @returns ab = volume / pi
        elemental function get_ab(volume) result(ab)
            double precision, intent(in) :: volume
            double precision             :: ab

            ab = volume * fpi
        end function

        ! Obtain the aspect ratio a/b of the parcel(s).
        ! @param[in] a2 is the largest eigenvalue
        ! @param[in] volume of the parcel(s)
        ! @returns lam = a/b
        elemental function get_aspect_ratio(a2, volume) result(lam)
            double precision, intent(in) :: a2
            double precision, intent(in) :: volume
            double precision             :: lam

            lam = (a2 / volume) * pi
        end function

        ! Obtain the ellipse support points for par2grid and grid2par
        ! @param[in] position vector of the parcel
        ! @param[in] volume of the parcel
        ! @param[in] B matrix elements of the parcel
        ! @returns the parcel support points
        function get_ellipse_points(position, volume, B) result(points)
            double precision, intent(in) :: position(2)
            double precision, intent(in) :: volume
            double precision, intent(in) :: B(2)        ! B11, B12
            double precision             :: B22
            double precision             :: c
            double precision             :: a2
            double precision             :: evec(2)
            double precision             :: h(2)
            double precision             :: points(2, 2)

            B22 = get_B22(B(1), B(2), volume)

            a2 = get_eigenvalue(B(1), B(2), B22)

            c = dsqrt(dabs(two * a2 - B(1) - B22))

            evec = get_eigenvector(a2, B(1), B(2), B22)

            h = f12 * c * evec

            points(:, 1) = position - h
            points(:, 2) = position + h

        end function get_ellipse_points
end module parcel_ellipse

max_stars_repo_name: matt-frey/epic

max_issues_repo_name: matt-frey/epic

max_forks_repo_name: matt-frey/epic


text:
      subroutine cappr(iz,ndummy)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C**********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 Calculate capillary pressure.
CD1
C**********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 05-20-92     G. Zyvoloski   00022   Initial implementation.
CD2                                        However, previous non-YMP
CD2                                        versions of FEHM exist, and
CD2                                        the current version differs
CD2                                        from these in minor ways.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/cappr.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:42:22   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 08:55:16   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:04:58   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:22:00   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 11:56:04   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:39:06 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.6   Wed Feb 14 11:35:26 1996   zvd
CD2 Updated purpose to reflect input now done with rlp
CD2
CD2    Rev 1.5   Mon Jan 29 13:16:58 1996   hend
CD2 Updated Requirements Traceability
CD2
CD2    Rev 1.4   12/13/95 08:30:04   gaz
CD2 cap pressure model inputted in macro rlp
CD2
CD2    Rev 1.3   01/23/95 10:58:26   robinson
CD2 BAR: Fixed bug in computation of linear cap. pressure model
CD2
CD2    Rev 1.2   09/02/94 12:14:28   robinson
CD2 Fixed bug in reading of parameters - BAR
CD2
CD2    Rev 1.1   03/18/94 16:11:24   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   01/20/94 10:21:32   pvcs
CD2 original version in process of being certified
CD2
C**********************************************************************
CD3
CD3 INTERFACES
CD3
CD3 Formal Calling Parameters
CD3
CD3 Identifier              Type     Use  Description
CD3
CD3 iz                      int       I   Control parameter used to
CD3                                          determine the reason for
CD3                                          this call
CD3 ndummy                  int       I   Parameter to set the correct
CD3                                          node number for dual
CD3                                          porosity or dpdp
CD3
CD3 Interface Tables
CD3
CD3 None
CD3
CD3 Files
CD3
CD3 Name           Description
CD3
CD3 input file     Contains all input data from FEHMN macros in ASCII
CD3                form
CD3
C**********************************************************************
CD4
CD4 GLOBAL OBJECTS
CD4
CD4 Global Constants
CD4
CD4 Identifier  Type     Description
CD4
CD4
CD4 Global Types
CD4
CD4 NONE
CD4
CD4 Global Variables
CD4
CD4 icapp, inpt, wdd1, icapt, cp1f, cp2f, cp3f, narrays, pointer, itype,
CD4 default, igroup, ireturn,
CD4 macroread
CD4
CD4                      COMMON
CD4 Identifier   Type    Block   Description
CD4
CD4
CD4
CD4 Global Subprograms
CD4
CD4 Name      Type     Description
CD4
CD4 initdata  N/A      Reads and sets input values defined at each node
CD4
C**********************************************************************
CD5
CD5 LOCAL IDENTIFIERS
CD5
CD5 Local Constants
CD5
CD5 None
CD5
CD5 Local Types
CD5
CD5 None
CD5
CD5 Local Variables
CD5
CD5 Identifier   Type        Description
CD5
CD5 i            int         Do loop parameter
CD5 mid          int         Do loop index
CD5 mi           int         Current node number
CD5 it           int         Current region number of capillary
CD5                             pressure model
CD5 itp          int         Current capillary pressure model number
CD5 itperm       int         Current model number for relative
CD5                             permeability model
CD5 itpperm      int         Current relative permeability model number
CD5 cp1          real*8      Capillary pressure model parameter
CD5 cp3          real*8      Capillary pressure model parameter
CD5 im           int         Do loop parameter
CD5
CD5 Local Subprograms
CD5
CD5 None
CD5
C**********************************************************************
CD6
CD6 FUNCTIONAL DESCRIPTION
CD6
CD6
C**********************************************************************
CD7
CD7 ASSUMPTIONS AND LIMITATIONS
CD7
CD7 None
CD7
C**********************************************************************
CD8
CD8 SPECIAL COMMENTS
CD8
CD8  Requirements from SDN: 10086-RD-2.20-00
CD8    SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD8    FEHM Application Version 2.20
CD8
C**********************************************************************
CD9
CD9 REQUIREMENTS TRACEABILITY
CD9
CD9 2.4.4 Relative-permeability and capillary-pressure functions
CD9
C**********************************************************************
CDA
CDA REFERENCES
CDA
CDA See FEHMN SRS, MMS, and SDD
CDA
C**********************************************************************
CPS
CPS PSEUDOCODE
CPS
CPS BEGIN cappr
CPS
CPS IF capillary pressures are used in this simulation
CPS
CPS   IF the data are to be read in on this pass
CPS
CPS     LOOP to read in first line of data
CPS
CPS       null - Check for null line or zero
CPS       IF there is more data to read
CPS         Read in capillary pressure parameters
CPS         IF the capillary pressure parameters were set correctly
CPS           Compute slope of linear capillary pressure model
CPS         ELSE the parameters are incorrect
CPS           Set slope to 0 to zero out capillary pressure at that node
CPS         ENDIF
CPS       ELSE
CPS         IF no data was present
CPS           Set parameter to neglect capillary pressures
CPS         ENDIF
CPS     EXITIF no more data in this group
CPS       ENDIF
CPS
CPS     ENDLOOP
CPS
CPS     Set parameters for reading capillary pressure flag
CPS
CPS     initdata - read flag and set at each node
CPS
CPS     Set flag to denote that the cap macro has been called
CPS
CPS   ELSE parameter values are being set
CPS
CPS     FOR each node
CPS       Set integer parameters
CPS       IF a van Ganuchten relative permeability model is used
CPS         Set parameter to skip capillary pressure calculation
CPS       ELSE
CPS         Set capillary pressure model flag
CPS       ENDIF
CPS       IF a linear model is assumed
CPS         Compute capillary pressure and derivative
CPS         IF the computed value is less than 0
CPS           Set value and derivative to 0
CPS         ENDIF
CPS       ENDIF
CPS     ENDFOR
CPS
CPS     FOR each node
CPS       IF this is a compressed liquid
CPS         Set capillary pressure and derivative to 0
CPS       ENDIF
CPS     ENDFOR
CPS
CPS   ENDIF the data are to be read in on this pass
CPS
CPS ENDIF capillary pressures are used in this simulation
CPS
CPS END cappr
CPS
C**********************************************************************
c
c subroutine to calculate the capillary pressure
c
      use combi
      use comci
      use comdi
      use comei
      use comfi
      use comgi
      use comii
      use comdti
      use comai
      use comki
      use comwt

      implicit none

      integer i,iz,ndummy,mid,it,mi,itp,itperm,itpperm,im, iz_temp
	real*8 hp, sl, dhp,rp1,rp2,rp3,rp4,sucut,slcut
      real*8 cp1,cp3,cap_max
      integer jj, ireg
	real*8 s_corr,ds_corrp,s_vg,ds_vgp
	real*8 ph_dum,hmin,hmax,pl,del_h,del_s,dr_nr
	real*8 ac3, ac4, bc3, bc4, hlcut, hucut, dslh
	real*8 smcut,smcutm,scutm,alamda,alpi,hcut,alpha,beta
c PJJohnson code changes - variable definitions
      real*8 pjni, pjSri, pjCpmaxi, pjn, pjlmax, pjsat, pjm, pjb
      real*8 pjSr, pjCpmax, pjconstant
c end PJJohnson variable changes
	parameter (scutm = 1.d-03)

c read in data
      if(icapp.eq.0) return
         if(iz.gt.0) then
c
c     load nodal capillary pressures
c
            do mid=1,neq
               mi=mid+ndummy
               it=icap(mi)
               itperm = irlp(mi)
               itpperm = irlpt(itperm)
              if(it.eq.0.or.(itpperm .ge. 3 .and. itpperm .le. 20)) then
                  itp = itpperm
              else
                  itp=icapt(it)
              end if
               if(itp.eq.1) then
c
c     linear forsythe(1988) model
c
c * * * * * * * * * * PJJOHNSON CODE MODIFICATIONS, SUMMER 2017 * * * * * * *
c added four flags for por-dependent capillary functions
c these are tied to the linear rlp model - user assigns as normal except Slmax is given a flag, -111, -333, -666, or -999
c which activates the f(sat, n) one
c pjn, pjni, pjsr, pjSri, pjcpmax, pjcpmaxi are variables used for that
c ALL PJJOHNSON variables used for calcs throughout FEHM are noted by prefix pj
c

                   if(irlpt(it).eq.-666) then
c set residual saturation and max capillary pressure
c                      open(666,file='Cp_out.txt')

                       pjSri= rp1f(it)
                       pjCpmaxi=cp1f(it)

c read porosity from global variables
                       pjn=ps(mi)
                      if(pjn.ge. 0.9 .or. pjn.le.1e-6) then
                          pcp(mi) = 0
                          dpcef(mi) = 0
                      else

c flag for initial porosity: if 0, then read it from global
c otherwise, user can specify value to use (e.g. for same material with different properties)
                      if(rp5f(it) .gt. 0) then
                          pjni = rp7f(it)
                      else

                       pjni = psini(mi)
                      endif

c flag for changing saturation above which Pc = 0; if 0, change it, otherwise stays same
                      if(cp3f(it).ne.0) then
                          pjlmax = 1
                      else
                          pjlmax = 1-pjn
                      endif
c linear residual saturation calculation
                       pjSr=(pjSri/(1-pjni))*(1-pjn)

c extrapolate max capillary pressure - note that this is at sat = Sr, not sat = 0 (different from rlp 1)
c                       pjCpmax=(pjCpmaxi/pjsri)*(pjsr)
                       pjCpmax=pjCpmaxi*(1-pjn)/(1-pjni)
c read saturation
                       pjsat=s(mi)
c                         pjb = pjCpmax + (pjCpmax/(pjlmax-pjSr)*pjSr)
c calculate Pc
                       if(pjsat.ge.pjlmax) then
                           pcp(mi)=0
                       else if(pjsat.le.pjSr) then
                          pcp(mi)=pjCpmax
                       else
                         pjm = -1*(pjCpmax / (pjlmax - pjSr))
c define derivative
                         pcp(mi) = pjCpmax*(pjlmax-pjsat)/(pjlmax-pjSr)
c
c
                          dpcef(mi) = pjm
c
                       end if
                       end if

                  else if(irlpt(it).eq. -333) then
c Leverett function

c if user enters a number for porosity entry in rlperm.f, use that
c if they enter 0, then use the global variable
c this allows them to have multiple different units be compared to same standard
                      if(rp7f(it) .gt. 0) then
                          pjni = rp7f(it)
                      else

                       pjni = psini(mi)

                      endif
c take initial maximum capillary pressure and residual saturation from input
                       pjCpmaxi = cp1f(it)

                       pjSri= rp1f(it)
c
c read current porosity from global variables
                       pjn=ps(mi)
c Permeability comes from salt macro, but is sometimes fed 0 (e.g. timestep 1)
c so address that if need be
                       if(pjk(mi) .eq. 0) then
                           pjk(mi) = pjki(it)
                       endif

c Calculate leverett capillary pressure
                       pjCpmax = pjCpmaxi*(sqrt(pjki(it)/pjni) /
     &                 (sqrt(pjk(mi)/pjn)))

                          pjlmax = 1

c handle residual saturation based on user preference
c if they use 0, hold Sr constant; otherwise, vary with porosity
                      if(cp3f(it).ne.0) then
                       pjSr=(pjSri/(1-pjni))*(1-pjn)

                      else
                          pjSr = pjSri
                      endif
c truncate Sr in case it tries to go >1 or to .99999 with function divergence issues
                       if(pjSr .ge. 0.99) then
                           pjSr = 0.99
                       endif
c read saturation
                       pjsat=s(mi)
c calculate Pc
c max sat is currently 1 but if modified, this line handles that
                       if(pjsat.ge.pjlmax) then
                           pcp(mi)=0
                           dpcef(mi)=0
c if Sl < Sr, use maximum value
                       else if(pjsat.le.pjSr) then
                          pcp(mi)=pjCpmax
                          dpcef(mi)=0
c otherwise fit saturation function
                       else
                         pjm = -1*(pjCpmax / (pjlmax - pjSr))
                         pcp(mi) = pjCpmax*(pjlmax-pjsat)/(pjlmax-pjSr)

c I added a truncation here to hold Pc to no more than an order of magnitude higher than the starting value
c to avoid crashes
                         if(pcp(mi) .gt. (10*pjCpmaxi)) then
                             pcp(mi) = 10*pjCpmaxi
c failsafe in case of model crash
                         else if (pcp(mi).lt.0) then
                             pcp(mi) = 0
                         else
                             pcp(mi) = pcp(mi)
                         end if

                         pjm = -1*(pjCpmax / (pjlmax - pjSr))
                         dpcef(mi) = pjm
                       end if
c following comment lines are output files that can be activated to check calcs
c                      write(333,*) pjni, pjn
c                      write(666,*) pjki(it), pjk(mi)
c                      write(999,*) pjCpmaxi, pjCpmax
c                      write(111,*) pjSri, pjSr
c                      write(777,*) pjsat, pcp(mi)
c gaz debug 011219
c                       end if

                  else
                  pjn = 0
                  pjsat = 0
c *********END PJJOHNSON CHANGES **********
c
c gaz debug 011219
                  endif
                  cp1=cp1f(it)
                  cp3=cp3f(it)
                  if(ieos(mi).eq.2) then
                     pcp(mi)=cp1*(cp3-s(mi))
                     dpcef(mi)=-cp1
                     if(pcp(mi).lt.0.0) then
                        pcp(mi)=0.0
                        dpcef(mi)=0.0
                     endif
                   else if(ieos(mi).eq.3) then
                     pcp(mi)=cp1*cp3
                     dpcef(mi)=0.0
                   else if(ieos(mi).eq.1) then
                     pcp(mi)=0.0
                     dpcef(mi)=0.0
                   endif

               end if
               if(itp.eq.3) then
c
c     van Genucten (1980) model (calculated in rlperm)
c
               endif
             continue
            enddo

         else if(iz.lt.0) then
c
c inverse models (031707 gaz)
c
c     load nodal capillary pressures
c

         do mid=1,neq
              mi=mid+ndummy
c inverse calculation valid only for wtsi nodes
              iz_temp = izone_free_nodes(mi)
	        if(iz_temp .gt. 1) then
               it=icap(mi)
               itperm = irlp(mi)
               itpperm = irlpt(itperm)
               if(itpperm .ge. 3 .and. itpperm .le. 20) then
                  itp = itpperm
               else
                  itp=icapt(it)
               end if
	   if(itp.eq.1) then
c
c     linear forsyth(1988) model
c     inverse solved explicitly
c
                  cp1=cp1f(it)
                  cp3=cp3f(it)
c  "+" removes negative cap pressure
                  hmin = head12(mi,1)+cp1*cp3
		        hmax = head12(mi,2)
	            del_h = hmax-hmin
	            pl = phi(mi)
	             if(irich.eq.0) then
c richard's equation  with wtsi correction
                    if(iz_temp.eq.2) then
                      rlxyf(mi) = (pl-(hmin-cp1*cp3))/
     &				 (del_h+cp1) + rlptol
			        drlxyf(mi) = 1.d0/(del_h+cp1)
				    s(mi) = rlxyf(mi) - rlptol
				  else
				   	rlxyf(mi) = 0.0
					s(mi) = 0.0
				    drlxyf(mi) = 1.d0/(del_h+cp1)
				  endif
	            else
c richard's equation only
                   call vgcap_inv_calc(2,mi)
	            endif

          else if(itp.eq.3) then
c
c     van Genucten (1980) model
c
c s = cap_inv  explicit calc
c hmin and hmax are "pressures"
c
c richard's equation  with wtsi correction
	      if(irich.eq.0) then
c get richards equation terms
	        call vgcap_inv_calc(2,mi)
c wtsi terms  (still needs work gaz-112007)
             	  hmin = head12(mi,1)
		      hmax = head12(mi,2)
	          del_h = hmax-hmin
	          pl = phi(mi)
                s_corr = (pl-hmin)/del_h
                ds_corrp = 1.d0/del_h
	          if(s_corr.lt.0.0d0) then
	            s_corr = 0.0d0
	            ds_corrp = 0.0d0
	          else if(s_corr.gt.1.0d0) then
	           s_corr = 1.0d0
	           ds_corrp = 0.0d0
 	          endif
c
c  gaz dslh is not defined
	           s_vg = s(mi)
                 ds_vgp = drlxyf(mi)
c
c Now the combination
c
			   	rlxyf(mi) = s_corr*(1.d0-s_vg) + s_vg + rlptol
                  s(mi) = rlxyf(mi) - rlptol
                  drlxyf(mi) = ds_corrp*(1.d0-s_vg)
     &				+ (1.d0-s_corr)*ds_vgp

	       else
c
c richard's equation only (test solution)
c
                 call vgcap_inv_calc(2,mi)
             endif
        endif
c end itp block
        endif
c end wtsi block
       enddo
      endif

      return

      end

ext: f

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   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.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
!!!!   work (education/research/industry/development/...) by citing the ParaMonte
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> \brief
!> This module contains the classes and procedures for setting up the `burninAdaptationMeasure` attribute of samplers of class [ParaDRAM_type](@ref paradram_mod::paradram_type).
!> For more information, see the description of this attribute in the body of the module.
!> \author Amir Shahmoradi

module SpecDRAM_BurninAdaptationMeasure_mod

    use Constants_mod, only: RK
    implicit none

    character(*), parameter         :: MODULE_NAME = "@SpecDRAM_BurninAdaptationMeasure_mod"

    real(RK)                        :: burninAdaptationMeasure ! namelist input

    type                            :: BurninAdaptationMeasure_type
        real(RK)                    :: val
        real(RK)                    :: def
        real(RK)                    :: null
        character(:), allocatable   :: desc
    contains
        procedure, pass             :: set => setBurninAdaptationMeasure, checkForSanity, nullifyNameListVar
    end type BurninAdaptationMeasure_type

    interface BurninAdaptationMeasure_type
        module procedure            :: constructBurninAdaptationMeasure
    end interface BurninAdaptationMeasure_type

    private :: constructBurninAdaptationMeasure, setBurninAdaptationMeasure, nullifyNameListVar

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function constructBurninAdaptationMeasure(methodName) result(BurninAdaptationMeasureObj)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: constructBurninAdaptationMeasure
#endif
        use Constants_mod, only: NULL_RK
        use String_mod, only: num2str
        implicit none
        character(*), intent(in)            :: methodName
        type(BurninAdaptationMeasure_type)  :: BurninAdaptationMeasureObj
        BurninAdaptationMeasureObj%def = 1._RK
        BurninAdaptationMeasureObj%null = NULL_RK
        BurninAdaptationMeasureObj%desc = &
        "burninAdaptationMeasure is a 64-bit real number between 0 and 1, representing the adaptation measure threshold below &
        &which the simulated Markov chain will be used to generate the output " //methodName// " sample. In other words, any &
        &point in the output Markov Chain that has been sampled during significant adaptation of the proposal distribution (as &
        &determined by burninAdaptationMeasure) will not be included in the construction of the final " //methodName// " output &
        &sample. This is to ensure that the generation of the output sample will be based on the part of the simulated chain that &
        &is practically guaranteed to be Markovian and ergodic. If this variable is set to 0, then the output sample will be &
        &generated from the part of the chain where no proposal adaptation has occurred. This non-adaptive or &
        &minimally-adaptive part of the chain may not even exist if the total adaptation period of the simulation &
        &(as determined by adaptiveUpdateCount and adaptiveUpdatePeriod input variables) is longer than the total length &
        &of the output MCMC chain. In such cases, the resulting output sample may have a zero size. In general, when &
        &good mixing occurs (e.g., when the input variable chainSize is very large) any specific value of burninAdaptationMeasure &
        &becomes practically irrelevant. &
        &The default value for burninAdaptationMeasure is "//num2str(BurninAdaptationMeasureObj%def)//", implying that the entire &
        &chain (with the exclusion of an initial automatically-determined burnin period) will be used to generate the final output sample."
    end function constructBurninAdaptationMeasure

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine nullifyNameListVar(BurninAdaptationMeasureObj)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: nullifyNameListVar
#endif
        implicit none
        class(BurninAdaptationMeasure_type), intent(in)  :: BurninAdaptationMeasureObj
        burninAdaptationMeasure = BurninAdaptationMeasureObj%null
    end subroutine nullifyNameListVar

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine setBurninAdaptationMeasure(BurninAdaptationMeasureObj,burninAdaptationMeasure)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: setBurninAdaptationMeasure
#endif
        use Constants_mod, only: RK
        implicit none
        class(BurninAdaptationMeasure_type), intent(inout)  :: BurninAdaptationMeasureObj
        real(RK), intent(in)                                :: burninAdaptationMeasure
        BurninAdaptationMeasureObj%val = burninAdaptationMeasure
        if (BurninAdaptationMeasureObj%val==BurninAdaptationMeasureObj%null) BurninAdaptationMeasureObj%val = BurninAdaptationMeasureObj%def
    end subroutine setBurninAdaptationMeasure

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine checkForSanity(BurninAdaptationMeasureObj,Err,methodName)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: checkForSanity
#endif
        use Constants_mod, only: RK
        use Err_mod, only: Err_type
        use String_mod, only: num2str
        implicit none
        class(BurninAdaptationMeasure_type), intent(in) :: BurninAdaptationMeasureObj
        type(Err_type), intent(inout)       :: Err
        character(*), intent(in)            :: methodName
        character(*), parameter             :: PROCEDURE_NAME = "@checkForSanity()"
        if ( BurninAdaptationMeasureObj%val<0._RK ) then
            Err%occurred = .true.
            Err%msg =   Err%msg // &
                        MODULE_NAME // PROCEDURE_NAME // ": Error occurred. &
                        &The input variable burninAdaptationMeasure (" // num2str(BurninAdaptationMeasureObj%val) // &
                        ") cannot be less than 0. If you are not sure of the appropriate value for burninAdaptationMeasure, drop it &
                        &from the input list. " // methodName // " will automatically assign an appropriate value to it.\n\n"
        end if
        if ( BurninAdaptationMeasureObj%val>1._RK ) then
            Err%occurred = .true.
            Err%msg =   Err%msg // &
                        MODULE_NAME // PROCEDURE_NAME // ": Error occurred. &
                        &The input variable burninAdaptationMeasure (" // num2str(BurninAdaptationMeasureObj%val) // &
                        ") cannot be larger than 1. If you are not sure of the appropriate value for burninAdaptationMeasure, drop it &
                        &from the input list. " // methodName // " will automatically assign an appropriate value to it.\n\n"
        end if
    end subroutine checkForSanity

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end module SpecDRAM_BurninAdaptationMeasure_mod ! LCOV_EXCL_LINE

max_stars_repo_name: cdslaborg/paramonte

max_issues_repo_name: cdslaborg/paramonte

max_forks_repo_name: cdslaborg/paramonte


text:
	SUBROUTINE MT_GRPT  ( bultin, lenb, bultyp, ibpnt,
     +			      rpttyp, report, lenr, corflg, iret )
C************************************************************************
C* MT_GRPT							        *
C*							 	        *
C* This subroutine gets the next report from a METAR/SPECI bulletin.    *
C* Reports must begin with a four-character station id, and should end  *
C* with '='.  Reports are checked to see if they conform to the 	*
C* SAO format.  The report type is saved to common.			*
C*								        *
C* MT_GRPT  ( BULTIN, LENB, BULTYP, IBPNT, RPTTYP, REPORT, LENR,        *
C* 	      CORFLG, IRET ) 						*
C*								        *
C* Input parameters:						        *
C*	BULTIN		CHAR*		Bulletin		        *
C*	LENB		INTEGER		Bulletin length		        *
C*	BULTYP		CHAR*		Type of bulletin                *
C*								        *
C* Input and output parameters:					        *
C*	IBPNT		INTEGER		Pointer in bulletin	        *
C*								        *
C* Output parameters:						        *
C*	RIVALS(IRTHRP)	REAL		Report type flag                *
C*	RPTTYP		CHAR*		Type of report                  *
C*	REPORT		CHAR*		Report			        *
C*	LENR		INTEGER		Length of report	        *
C*	CORFLG		LOGICAL		Correction flag                 *
C*	IRET		INTEGER		Return code		        *
C*					  0 = normal return	        *
C*					 -1 = no station id found       *
C*					 -2 = no more reports	        *
C*					 -3 = NIL report		*
C**								        *
C* Log:								        *
C* D. Kidwell/NCEP     	10/95      				        *
C* D. Kidwell/NCEP 	12/95  	Check for id at start of US reports	*
C* D. Kidwell/NCEP 	 3/96	Simplified station id search          	*
C* D. Kidwell/NCEP 	 8/96	Removed references to TESTM, TESTS    	*
C* K. Tyle/GSC		11/96	Check for SAO's			 	*
C* K. Tyle/GSC		 1/97	Reorganize header and comments		*
C* K. Tyle/GSC		 3/97	Check length of station id; increase	*
C*				dimension of stid to 8			*
C* K. Tyle/GSC		 4/97	Check for NIL reports and reports that	*
C*				begin with WMO bulletin ID		*
C* D. Kidwell/NCEP 	 6/97	Removed call to ST_LSTR for stid; added *
C*				explicit character length specification *
C* D. Kidwell/NCEP 	 6/97	Cleaned up code                         *
C* D. Kidwell/NCEP 	 4/98	New interface; removed iasos flag       *
C* D. Kidwell/NCEP 	 9/05	CSC to add corflg                       *
C* S. Chiswell/Unidata	10/07	Added LWIS check			*
C************************************************************************
	INCLUDE		'mtcmn.cmn'
C*
	CHARACTER*(*)	bultin, report, bultyp, rpttyp
	LOGICAL		corflg
C*
	LOGICAL 	found, saoflg, wmo
	CHARACTER	type*8, satype*3, stid*8
C-----------------------------------------------------------------------
        iret   = 0
	lenr   = 0
	saoflg = .false.
	wmo    = .true.
	corflg = .false.
C
C*	Check for end of bulletin.
C
        IF ( ibpnt + 7 .ge. lenb ) THEN
	    report = ' '
            iret   = -2
            RETURN
        END IF
C
C*	Get report type.
C
	rpttyp = ' '
	ibeg   = ibpnt
	type   = bultin ( ibeg:ibeg + 7 )
	indxmt = INDEX ( type, 'METAR' )
	indxsp = INDEX ( type, 'SPECI' )
	indxlw = INDEX ( type, 'LWIS' )
C
C*	Check for SAO format--found in pos. 5-6 in Canadian
C*	and 6-7 in Alaskan rpts.
C
	satype = bultin ( ibeg + 4: ibeg + 6 )
	indxsa = INDEX ( satype, 'SA' )
	indxap = INDEX ( satype, 'SP' )
	IF ( ( indxmt .ne. 0 ) .or.
     +	     ( indxlw .ne. 0 ) ) THEN
	    rpttyp = 'MT'
	  ELSE IF ( indxsp .ne. 0 ) THEN
	    rpttyp = 'SP'
	  ELSE IF ( indxsa .ne. 0 ) THEN
	    rpttyp = 'SA'
	    saoflg = .true.
	  ELSE IF ( indxap .ne. 0 ) THEN
	    rpttyp = 'SX'
	    saoflg = .true.
	END IF
C
	IF ( rpttyp .eq. ' ' ) THEN
	    rpttyp = bultyp
	    IF ( bultyp .eq. ' ' ) rpttyp = 'MT'
	  ELSE IF ( indxap .eq. 0 .and. indxsa .eq. 0 ) THEN
	    IF ( indxlw .eq. 0 ) THEN
	      indx = MAX ( indxmt, indxsp )
	      ibeg = ibeg + indx + 5
	    ELSE
	      indx = indxlw
	      ibeg = ibeg + indx + 4
	    END IF
	END IF
C
C*      Find the next station id.
C
	DO WHILE ( wmo )
	    IF ( rpttyp .ne. ' ' ) THEN
C
C*		Check length of station id.
C
		iblk = INDEX ( bultin ( ibeg:ibeg+8 ), ' ' )
		IF ( iblk .eq. 0 ) THEN
		    stid = ' '
		  ELSE
		    stid = bultin ( ibeg:ibeg + iblk - 1 )
		END IF
		iens = ibeg + iblk
		lens = iblk - 1
		IF ( lens .eq. 4 ) THEN
		    found = .true.
		    wmo   = .false.
		  ELSE IF ( lens .eq. 3 .and. saoflg ) THEN
C
C*		    Canadian SAO report.
C
		    found = .true.
		    wmo   = .false.
		  ELSE
		    found = .false.
C
C*		    Check if report is mistakenly started with the WMO
C*		    bulletin ID, or for COR preceding the station id.
C
		    IF ( lens .gt. 0 ) THEN
		        IF ( stid .eq. wmohdr ( :lens ) ) THEN
			    ibeg = ibeg + lens + 1
			  ELSE IF ( stid .eq. 'COR' ) THEN
C
C*			    COR precedes the station id (WMO 11/05).
C
			    corflg = .true.
			    ibeg   = ibeg + lens + 1
		          ELSE
			    wmo  = .false.
			END IF
		      ELSE IF ( lens .eq. 0 ) THEN
C
C*			Leading character is a blank, so keep looking.
C
			ibeg = ibeg + 1
		      ELSE
			wmo  = .false.
		    END IF
		END IF
	      ELSE
		found = .false.
		wmo   = .false.
	    END IF
	END DO
C
C*      Look for '=' as report terminator.
C
	IF ( ibeg .gt. lenb ) THEN
	    iret = -2
	    RETURN
	  ELSE
	    iend1 = INDEX ( bultin ( ibeg:lenb ), '=' )
	    IF ( iend1 .ne. 0 ) THEN
		iend = iend1 + ibeg - 2
	      ELSE
		iend = lenb
	    END IF
	END IF
C
C*	Reset pointer in bulletin.
C
	ibpnt  = iend + 3
C
C*	Construct report and pointers.
C
	IF ( ibeg .gt. iend ) THEN
	    iend = ibpnt
	    IF ( ibeg .gt. ibpnt ) THEN
		iret = -2
		RETURN
	    END IF
	END IF
	report = bultin ( ibeg : iend )
	CALL ST_LSTR ( report ( :iend - ibeg + 1 ), lenr, ier )
C
C*	Don't decode a report containing "NIL" after station ID.
C
	lenmin = MIN ( lenb, iens + 10 )
	inil   = INDEX ( bultin ( iens:lenmin ), 'NIL' )
	IF ( inil .ne. 0 ) THEN
	    iret = -3
	    CALL DC_WLOG ( 4, 'DCMETR', 6, report(:lenr), ier)
	    RETURN
	END IF
C
C*	Was a station id found?
C
	IF ( .not. found ) THEN
C
C*	    No station id found where expected.
C
	    iret = -1
	  ELSE
C
C*	    Have a four character alpha string - assume station id.
C*	    Initialize common variables to missing.
C
	    CALL MT_INIT ( mret )
C
C*	    Store report type.
C
	    IF ( rpttyp .eq. 'MT' ) THEN
		rivals ( irthrp ) = 0.
	      ELSE IF ( rpttyp .eq. 'SP' ) THEN
		rivals ( irthrp ) = 1.
	    END IF
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text: California Networking appears to have taken over the location once occupied by D Z Solutions, a company that was equally mysterious.  Listings elsewhere on the interwebs seem to indicate that theyre called California Networking (you probably already figured that one out), that theyre located in Davis (if youre good, you got that one, too!), and that they focus on billing services (GOTCHA!).  Not sure how strong a connection there is between networking and billing services, but there you have it.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module file_reader
    use element
    implicit none
contains
    function read_file(path)
        character(*), intent(in) :: path
        integer, parameter :: number = 10
        type(string_element), allocatable, dimension(:) :: read_file
        character(128), allocatable :: line
        character(:), allocatable :: trimmed
        integer :: count, i
        line = ""
        trimmed = ""
        count = 0
        i = 1

        open(number, file=path, status='old')
        do
            read(number, *, end=100)
            count = count + 1
        end do
100     close(number)

        allocate(read_file(count))
        open(number, file=path, status='old')
        do
            read(number, "(A)", end=200)  line
            trimmed = trim(line)
            read_file(i)%string = trimmed
            i = i + 1
        end do
200     close(number)
    end function
end module


text:
The Fall 2012 ASUCD Election was held November 13th16th online at http://elections.ucdavis.edu. There were 14 candidates competing for 6 ASUCD Senate seats. Two slates, SMART and NOW, were responsible for 10 of the candidates. The additional 4 candidates ran independent.  NOW received 40% of the vote, SMART received 36% and the Independents received 24%.  Unlike the Winter 2012 ASUCD Election the six top vote getters and the six elected senators were the same.

Voter turnout was http://www.jtleathers.com/asucd/voterturnout.html just over 20%, the highest turnout for a Fall ASUCD Election since 2002.

In order, the winners were:
1. Armando Figueroa
2. Alyson Sagala
3. Tal Topf
4. Felicia Ong
5. Liam Burke
6. Maxwell Kappes

 Beergate

## The name of this scandal would have been Steinsgate but I didnt think anyone would get wiki:wikipedia:Steins;Gate the reference.
After the election, the dark truth about the SMART slate finally came out: some of them drink beer at parties.  But the worst part about this was the fact that not a single one of them cited their expenses for the booze on their campaign expenditure forms.  So Senator Bradley Bottoms, respectable citizen that he is, filed a complaint against the slate with the Elections Committee.  The SMARTies claimed that the campaign event was BYOB, and thus the hopsladen beverages were personal expenses only.

Would it be that this sad tale ended here, I would be most happy.  But woe to me, dear reader, as I am compelled to write the horrible tale of the unfortunate events that followed.  In fact, http://www.theaggie.org/2012/11/28/complaintsfiledagainstasucdsenateslates/ according to the Aggie, there were no less than eight electoral complaints filed in this election.  Their rivals in NOW also stood accused of the heinous crime of not reporting their expenses on the brewski line item.

Furthermore, candidate Felicia Ong accused a campaigner of her rival, Armando Figueroa, of calling her a racist in the supposed safety of a Unitrans bus.  This cruel remark was clearly possibly maybe a violation of the ASUCD Bill of Rights.  To add insult to insult, no one in ASUCD actually knows how the Bill of Rights figures into anything  it seems to have been adopted without adding it to any existing document, so its legal force is questionable at best.  As the Elections Committee only enforces Elections Codes, this sadly left Ong with no redress except winning a Senate seat.

Complaints were flying:  Did NOW repeat the tragedy of the Campaigning in Dorms Controversy?  Did SMARTs campaigning near the bus terminal count as an illegal endorsement from Unitrans?  (Okay, thats ones a definite no.)  As Jared CrisologoSmith put it so wisely, This is an outrage, and represents an wiki:tvtropes:egregious offense to the principles of democracy that ASUCD was founded upon.  It appeared that not one candidate stayed clean in this sad affair.

And then the other shoe finally fell  and I warn you dear reader, this is not a pretty sight  Alyson Noele Sagala was disqualified by the Elections Committee.  Shunned as a scapegoat, sacrificed to pay for all of the evils of the election  unreported Natty Light (Correction, Natty Light was not provided! We have standards), statements of opinion, and shaking hands at the bus stop  she was condemned to be banished from the ASUCD Senate.

Seeing all of the injustice, IAC Chair Sergio Cano took the case to the highest levels (well, the only other level), and filed an appeal with the ASUCD Court.  He claimed that even the Elections Committee had ignored the bylaws.  In so doing so, it meant that no new Senator could be seated, for all would grind to a halt.  And so it was that a summit was arranged; Cano convened a council with Justin Goss,  Aaron Hsu, the Elections Committee, and even the ASUCD Business Manager in hopes of settling the tragic tale once and for all.  Serg pleaded his case: most of the complaints had been filed against slates, and not against candidates.  The bylaws specify that only candidates can be charged or censured  not slates  in order to prevent far wiki:berkeley:ASUC greater travesties.

Canos argument won the day, and all of the election complaints were dropped.  Everyone was forgiven, and allowed to take office!  It was like the whole affair had never happened, except that these people are http://www.theaggie.org/2012/12/04/editorialsenate/ now supposed to work with each other every day.  Huzzah!

 The California Aggie Aggie Endorsements
1. Olivia Brown
2. Liam Burke
3. Armando Figueroa
4. Felicia Ong
5. Jonathan Yip
6. Alyson Sagala

 Alternate Universe
For the following sections bold means they were one of the winning candidates. These are not the actual election results, but fun we can have with the data. The usual disclaimers apply, such as voter behavior being different in a system without choice voting, etc.

Condorcet Winner
The http://en.wikipedia.org/wiki/Condorcet_method Condorcet winner (the least objectionable candidate) of this election was Alyson Sagala.

 Top 6
What if the top six votegetters were elected, counting peoples top 6 rankings?

 Most Hated
Who got the most lastranked votes?

Countback Winners
These Countback elections were run under the assumption that all eligible candidates were willing to serve if elected in the Countback election and that there had been no prior Countback election for the Fall 2012 pool of Senate candidates. Both of those assumptions can alter the outcome of an actual Countback election. This list shows the current Senators and who would replace them in a Countback election if both of the preceding assumptions were true.

  Armando Figueroa would be replaced by Olivia Brown
  Alyson Noele Sagala would be replaced by Olivia Brown
  Tal Topf would be replaced by Kevin Pelstring
  Felicia Ong would be replaced by Robert Helfend
  Liam Burke would be replaced by Robert Helfend (Helfend defeats Pelstring by less than 1 vote in the final round of the Countback election)
  Maxwell Kappes would be replaced by Robert Helfend

The ASUCD Constitution requires that the winner of a Countback election acquire at least 50% of the threshold for the Countback election to serve as a Senator. All six candidates met that requirement.

 Links
  http://www.jtleathers.com/asucd/2012_Fall/fall2012overview.html Fall 2012 Ballot Data and Reports


text:
      SUBROUTINE MB01RX( SIDE, UPLO, TRANS, M, N, ALPHA, BETA, R, LDR,
     $                   A, LDA, B, LDB, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute either the upper or lower triangular part of one of the
C     matrix formulas
C        _
C        R = alpha*R + beta*op( A )*B,                               (1)
C        _
C        R = alpha*R + beta*B*op( A ),                               (2)
C                                             _
C     where alpha and beta are scalars, R and R are m-by-m matrices,
C     op( A ) and B are m-by-n and n-by-m matrices for (1), or n-by-m
C     and m-by-n matrices for (2), respectively, and op( A ) is one of
C
C        op( A ) = A   or   op( A ) = A',  the transpose of A.
C
C     The result is overwritten on R.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     SIDE    CHARACTER*1
C             Specifies whether the matrix A appears on the left or
C             right in the matrix product as follows:
C                     _
C             = 'L':  R = alpha*R + beta*op( A )*B;
C                     _
C             = 'R':  R = alpha*R + beta*B*op( A ).
C
C     UPLO    CHARACTER*1                               _
C             Specifies which triangles of the matrices R and R are
C             computed and given, respectively, as follows:
C             = 'U':  the upper triangular part;
C             = 'L':  the lower triangular part.
C
C     TRANS   CHARACTER*1
C             Specifies the form of op( A ) to be used in the matrix
C             multiplication as follows:
C             = 'N':  op( A ) = A;
C             = 'T':  op( A ) = A';
C             = 'C':  op( A ) = A'.
C
C     Input/Output Parameters
C
C     M       (input) INTEGER           _
C             The order of the matrices R and R, the number of rows of
C             the matrix op( A ) and the number of columns of the
C             matrix B, for SIDE = 'L', or the number of rows of the
C             matrix B and the number of columns of the matrix op( A ),
C             for SIDE = 'R'.  M >= 0.
C
C     N       (input) INTEGER
C             The number of rows of the matrix B and the number of
C             columns of the matrix op( A ), for SIDE = 'L', or the
C             number of rows of the matrix op( A ) and the number of
C             columns of the matrix B, for SIDE = 'R'.  N >= 0.
C
C     ALPHA   (input) DOUBLE PRECISION
C             The scalar alpha. When alpha is zero then R need not be
C             set before entry.
C
C     BETA    (input) DOUBLE PRECISION
C             The scalar beta. When beta is zero then A and B are not
C             referenced.
C
C     R       (input/output) DOUBLE PRECISION array, dimension (LDR,M)
C             On entry with UPLO = 'U', the leading M-by-M upper
C             triangular part of this array must contain the upper
C             triangular part of the matrix R; the strictly lower
C             triangular part of the array is not referenced.
C             On entry with UPLO = 'L', the leading M-by-M lower
C             triangular part of this array must contain the lower
C             triangular part of the matrix R; the strictly upper
C             triangular part of the array is not referenced.
C             On exit, the leading M-by-M upper triangular part (if
C             UPLO = 'U'), or lower triangular part (if UPLO = 'L') of
C             this array contains the corresponding triangular part of
C                                 _
C             the computed matrix R.
C
C     LDR     INTEGER
C             The leading dimension of array R.  LDR >= MAX(1,M).
C
C     A       (input) DOUBLE PRECISION array, dimension (LDA,k), where
C             k = N  when  SIDE = 'L', and TRANS = 'N', or
C                          SIDE = 'R', and TRANS = 'T';
C             k = M  when  SIDE = 'R', and TRANS = 'N', or
C                          SIDE = 'L', and TRANS = 'T'.
C             On entry, if SIDE = 'L', and TRANS = 'N', or
C                          SIDE = 'R', and TRANS = 'T',
C             the leading M-by-N part of this array must contain the
C             matrix A.
C             On entry, if SIDE = 'R', and TRANS = 'N', or
C                          SIDE = 'L', and TRANS = 'T',
C             the leading N-by-M part of this array must contain the
C             matrix A.
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,l), where
C             l = M  when  SIDE = 'L', and TRANS = 'N', or
C                          SIDE = 'R', and TRANS = 'T';
C             l = N  when  SIDE = 'R', and TRANS = 'N', or
C                          SIDE = 'L', and TRANS = 'T'.
C
C     B       (input) DOUBLE PRECISION array, dimension (LDB,p), where
C             p = M  when  SIDE = 'L';
C             p = N  when  SIDE = 'R'.
C             On entry, the leading N-by-M part, if SIDE = 'L', or
C             M-by-N part, if SIDE = 'R', of this array must contain the
C             matrix B.
C
C     LDB     INTEGER
C             The leading dimension of array B.
C             LDB >= MAX(1,N), if SIDE = 'L';
C             LDB >= MAX(1,M), if SIDE = 'R'.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The matrix expression is evaluated taking the triangular
C     structure into account. BLAS 2 operations are used. A block
C     algorithm can be easily constructed; it can use BLAS 3 GEMM
C     operations for most computations, and calls of this BLAS 2
C     algorithm for computing the triangles.
C
C     FURTHER COMMENTS
C
C     The main application of this routine is when the result should
C     be a symmetric matrix, e.g., when B = X*op( A )', for (1), or
C     B = op( A )'*X, for (2), where B is already available and X = X'.
C
C     CONTRIBUTORS
C
C     V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1999.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Mar. 2004.
C
C     KEYWORDS
C
C     Elementary matrix operations, matrix algebra, matrix operations.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      CHARACTER         SIDE, TRANS, UPLO
      INTEGER           INFO, LDA, LDB, LDR, M, N
      DOUBLE PRECISION  ALPHA, BETA
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), R(LDR,*)
C     .. Local Scalars ..
      LOGICAL           LSIDE, LTRANS, LUPLO
      INTEGER           J
C     .. External Functions ..
      LOGICAL           LSAME
      EXTERNAL          LSAME
C     .. External Subroutines ..
      EXTERNAL          DGEMV, DLASCL, DLASET, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         MAX
C     .. Executable Statements ..
C
C     Test the input scalar arguments.
C
      INFO   = 0
      LSIDE  = LSAME( SIDE,  'L' )
      LUPLO  = LSAME( UPLO,  'U' )
      LTRANS = LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' )
C
      IF(      ( .NOT.LSIDE  ).AND.( .NOT.LSAME( SIDE,  'R' ) ) )THEN
         INFO = -1
      ELSE IF( ( .NOT.LUPLO  ).AND.( .NOT.LSAME( UPLO,  'L' ) ) )THEN
         INFO = -2
      ELSE IF( ( .NOT.LTRANS ).AND.( .NOT.LSAME( TRANS, 'N' ) ) )THEN
         INFO = -3
      ELSE IF( M.LT.0 ) THEN
         INFO = -4
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDR.LT.MAX( 1, M ) ) THEN
         INFO = -9
      ELSE IF( LDA.LT.1 .OR.
     $   ( ( (      LSIDE .AND. .NOT.LTRANS ) .OR.
     $       ( .NOT.LSIDE .AND.      LTRANS ) ) .AND. LDA.LT.M ) .OR.
     $   ( ( (      LSIDE .AND.      LTRANS ) .OR.
     $       ( .NOT.LSIDE .AND. .NOT.LTRANS ) ) .AND. LDA.LT.N ) ) THEN
         INFO = -11
      ELSE IF( LDB.LT.1 .OR.
     $       (      LSIDE .AND. LDB.LT.N ) .OR.
     $       ( .NOT.LSIDE .AND. LDB.LT.M ) ) THEN
         INFO = -13
      END IF
C
      IF ( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'MB01RX', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( M.EQ.0 )
     $   RETURN
C
      IF ( BETA.EQ.ZERO .OR. N.EQ.0 ) THEN
         IF ( ALPHA.EQ.ZERO ) THEN
C
C           Special case alpha = 0.
C
            CALL DLASET( UPLO, M, M, ZERO, ZERO, R, LDR )
         ELSE
C
C           Special case beta = 0 or N = 0.
C
            IF ( ALPHA.NE.ONE )
     $         CALL DLASCL( UPLO, 0, 0, ONE, ALPHA, M, M, R, LDR, INFO )
         END IF
         RETURN
      END IF
C
C     General case: beta <> 0.
C     Compute the required triangle of (1) or (2) using BLAS 2
C     operations.
C
      IF( LSIDE ) THEN
         IF( LUPLO ) THEN
            IF ( LTRANS ) THEN
               DO 10 J = 1, M
                  CALL DGEMV( TRANS, N, J, BETA, A, LDA, B(1,J), 1,
     $                        ALPHA, R(1,J), 1 )
   10          CONTINUE
            ELSE
               DO 20 J = 1, M
                  CALL DGEMV( TRANS, J, N, BETA, A, LDA, B(1,J), 1,
     $                        ALPHA, R(1,J), 1 )
   20          CONTINUE
            END IF
         ELSE
            IF ( LTRANS ) THEN
               DO 30 J = 1, M
                  CALL DGEMV( TRANS, N, M-J+1, BETA, A(1,J), LDA,
     $                        B(1,J), 1, ALPHA, R(J,J), 1 )
   30          CONTINUE
            ELSE
               DO 40 J = 1, M
                  CALL DGEMV( TRANS, M-J+1, N, BETA, A(J,1), LDA,
     $                        B(1,J), 1, ALPHA, R(J,J), 1 )
   40          CONTINUE
            END IF
         END IF
C
      ELSE
         IF( LUPLO ) THEN
            IF( LTRANS ) THEN
               DO 50 J = 1, M
                  CALL DGEMV( 'NoTranspose', J, N, BETA, B, LDB, A(J,1),
     $                        LDA, ALPHA, R(1,J), 1 )
   50          CONTINUE
            ELSE
               DO 60 J = 1, M
                  CALL DGEMV( 'NoTranspose', J, N, BETA, B, LDB, A(1,J),
     $                        1, ALPHA, R(1,J), 1 )
   60          CONTINUE
            END IF
         ELSE
            IF( LTRANS ) THEN
               DO 70 J = 1, M
                  CALL DGEMV( 'NoTranspose', M-J+1, N, BETA, B(J,1),
     $                        LDB, A(J,1), LDA, ALPHA, R(J,J), 1 )
   70           CONTINUE
            ELSE
               DO 80 J = 1, M
                  CALL DGEMV( 'NoTranspose', M-J+1, N, BETA, B(J,1),
     $                        LDB, A(1,J), 1, ALPHA, R(J,J), 1 )
   80          CONTINUE
            END IF
         END IF
      END IF
C
      RETURN
C *** Last line of MB01RX ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
       subroutine gensdp_switch
!***********************************************************************
! Copyright 2008 Los Alamos National Security, LLC  All rights reserved
! Unless otherwise indicated,  this information has been authored by an
! employee or employees of the Los Alamos National Security, LLC (LANS),
! operator of the  Los  Alamos National  Laboratory  under Contract  No.
! DE-AC52-06NA25396  with  the U. S. Department  of  Energy.  The  U. S.
! Government   has   rights  to  use,  reproduce,  and  distribute  this
! information.  The  public may copy  and  use this  information without
! charge, provided that this  Notice and any statement of authorship are
! reproduced on all copies.  Neither  the  Government nor LANS makes any
! warranty,   express   or   implied,   or   assumes  any  liability  or
! responsibility for the use of this information.
!***********************************************************************
!D1
!D1  PURPOSE
!D1
!D1  This subroutine solves the isothermal air-water equations with
!D1  full jacobian (unsymmetric, 2n by 2n). Richards Equation
!D1
!***********************************************************************
!D2 REVISION HISTORY
!D2
!D2 Initial implementation: Oct-2008, Programmer: G. Zyvoloski
!D2
!D2 $Log:   /pvcs.config/fehm90/src/gensdp_switch.f_a  $
!D2
!***********************************************************************
!D3
!D3  REQUIREMENTS TRACEABILITY
!D3
!D3  2.3.2 Heat- and mass-transfer equations
!D3  2.3.3 Noncondensible gas flow equations
!D3  2.5.2 Solve nonlinear equation set at each time step
!D3
!***********************************************************************
!D4
!D4  SPECIAL COMMENTS AND REFERENCES
!D4
!D4 Requirements from SDN: 10086-RD-2.20-00
!D4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
!D4   FEHM Application Version 2.20
!D4
!***********************************************************************

      use davidi
      use comhi
      use comgi
      use comfi
      use comei
      use comdi
      use comci
      use combi
      use comdti
      use comai
      implicit none

      integer nc(16),inorm
      integer nmatd,ndex(6)
      integer i, id, idl, neqp1, nrhs1, nrhs2, nsizea
      integer j, i1, i2
      real*8, allocatable :: dumn(:)
      real*8, allocatable :: sto1(:)
      real*8, allocatable :: sto5(:)
      real*8, allocatable :: dum(:)
      real*8, allocatable :: piv_gaz(:,:)
      real*8  facr, fdm, fdum2, tollr, tolls
      real*8 a_total
      parameter(fdm=20.0, inorm = 1)
c
c     call thermo for next level (i+neq)
c
      call dpdp(2)
c
c     adjust storage
c
      nstorepiv = neq*4
      neqp1=neq+1
      fdum2=0.
      nmatd=nelm(neqp1)-neqp1
      nsizea=8*nmatd
      do i=1,nsizea
         a(i)=0.0
      enddo
      do i=1,neq
         bp(i+nrhs(1))=0.0
         bp(i+nrhs(2))=0.0
         bp(i+nrhs(3))=0.0
         bp(i+nrhs(4))=0.0
      enddo
      nrhs1=nrhs(1)
      nrhs2=nrhs(2)
      if(.not.allocated(piv_gaz)) then
         allocate(piv_gaz(neq,4))
      endif
      do id=1,neq
         nmat(1)=0
         nmat(2)=nmatd
c         nmat(3)=nmatd*(5-1)
c         nmat(4)=nmatd*(6-1)
         nrhs(1)=nrhs1
         nrhs(2)=nrhs2
c
c     decide on equation type
c
         if(ico2.lt.0) then
            if(ieos(id).ne.4) then
               call geneq2_switch( id)
               call add_accumulation(id)
            else
               call geneq3( id)
               a(nelmdg(id)+nmat(1))=sx1(id)
               bp(id+nrhs(1))=0.0
            endif
         else if (ico2.eq.0) then
            if(ieos(id).ne.4) then
               call geneq1( id)
            else
               call geneq3( id)
               a(nelmdg(id)+nmat(1))=sx1(id)
               bp(id+nrhs(1))=0.0
            endif
         endif
c
c     decide on equation type(for node idl=id+neq)
c
         idl=id+neq
         nmat(1)=nmatd*(7-1)
         nmat(2)=nmatd*(8-1)
c         nmat(3)=nmatd*(15-1)
c         nmat(4)=nmatd*(16-1)
         nrhs(1)=nrhs(3)
         nrhs(2)=nrhs(4)
         if(ico2.lt.0) then
            if(ieos(idl).ne.4) then
               call geneq2_switch(idl)
               call add_accumulation(idl)
            else
               call geneq3(idl)
               a(nelmdg(id)+nmat(1))=sx1(id)
               bp(id+nrhs(1))=0.0
            endif
         else if (ico2.eq.0) then
            if(ieos(idl).ne.4) then
               call geneq1(idl)
            else
               call geneq3(idl)
               a(nelmdg(id)+nmat(1))=sx1(id)
               bp(id+nrhs(1))=0.0
            endif
         endif
         continue
      enddo
      do i=1,8
         nmat(i)=nmatd*(i-1)
      enddo
      nrhs(1)=nrhs1
      nrhs(2)=nrhs2
c
c     call md_modes to complete equations for multiply defined nodes
c
c     if(imdnode.ne.0) call md_nodes(3,0,0)
c
c     now get coupling for influence terms
c
      if(ico2.eq.0) then
         call dpdpfh
      else if(ico2.lt.0) then
         call dpdpfa
      endif

c
c     arrange arrays for 2 dof solve (ie variable switching)
c
      call air_combine(1)
c
c     normalize the equations
c
c     some changes to call sequence here
c
      if(islord.ne.0) then
         call switch(nmat,nmatb,islord,2,1,nmatd)
         call switchb(bp,nrhs,nrhsb,islord,2,1,neq)
      end if

      ndex(1) = 1
      if(inorm.ne.0) then
         allocate(dumn(36))
c     call nrmlz4(neq,a,nmat,bp,nrhs,nelm,fdum2,sto1)
         call normal_dof(neq,a,bp,nelm,nmat,nrhs,nelmdg
     &        ,ndex,2,dumn(1),dumn(17),dumn(33),0,fdum2)
         deallocate(dumn)
      else
	 fdum2 = 0.0
	 do i = 1,neq
            fdum2 = fdum2 + bp(i+nrhs(1))*bp(i+nrhs(1)) +
     &           bp(i+nrhs(2))*bp(i+nrhs(2))
	 enddo
      endif

c
c     return and stop if singular
c
      if(ndex(1).lt.0 ) then
         if (iout .ne. 0) then
            write(iout,*) '   '
            write(iout, 100)
            write(iout,*) '   '
c            write(iout,'(i8,3g12.3)')
c     &           abs(nmatb(1)),(cord(abs(nmatb(1)),i),i=1,3)
            write(iout,'(i8,3g12.3)')
     &           abs(ndex(1)),(cord(abs(ndex(1)),i),i=1,3)
         end if
         if(iptty.gt.0) then
            write(iptty,*) '   '
            write(iptty,100)
            write(iptty,*) '   '
c            write(iptty,'(i8,3g12.3)')
c     &           abs(nmatb(1)),(cord(abs(nmatb(1)),i),i=1,3)
            write(iptty,'(i8,3g12.3)')
     &           abs(ndex(1)),(cord(abs(ndex(1)),i),i=1,3)
         endif
         iad=maxit
         return
      endif
 100  format ('* singular matrix found during normalization *')

      fdum=sqrt(fdum2)
      if(fdum.eq.0.0) go to 999
      mink=n
      if(iad.eq.0) then
         f0=max(fdum*epe,tmch)
      endif
      if(fdum1.lt.0.0.and.iad.ne.0) then
         do i=1,neq
            if(abs(bp(i+nrhs(1))).gt.tmch) go to 99
            if(abs(bp(i+nrhs(2))).gt.tmch) go to 99
         enddo
         fdum=-1.0
         go to 999
 99      continue
         f0=-1.0
      endif
      if(fdum.le.f0.and.iad.ne.0) goto 999
      facr=1.0
      tolls=max(facr*f0,min(g1*fdum,g2*fdum2))
      tollr=tolls*g3
      if(g3*tmch.gt.tollr) tollr=g3*tmch
c
c     set maximum iterations in solve_new
c

      iter=maxsolve
      call storage_derivatives(1,1)
c     allocate(dum(neq*4*4))
      allocate(dum(neq*8))
      if(irdof.eq.0) then
c
c     full solution
c

         if (igauss .gt. 1) then
            call solve_new(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv_gaz
     *           ,h,c,ss,g,y,iter,iback,2,iptty,maxor,accm)
         else
            call solve_new(neq,a,b,bp,nmat,nmat,nrhs,nelm,nelm,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv_gaz
     *           ,h,c,ss,g,y,iter,iback,2,iptty,maxor,accm)
         end if
      else
         allocate(sto5(neq*4))
         if (igauss .gt. 1) then
            call  rdof_new(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv_gaz
     *           ,h,c,ss,g,y,iter,iback,2,iptty,maxor
     *           ,overf,irdof,icoupl,0,sto5,accm)
         else
            call  rdof_new(neq,a,b,bp,nmat,nmat,nrhs,nelm,nelm,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv_gaz
     *           ,h,c,ss,g,y,iter,iback,2,iptty,maxor
     *           ,overf,irdof,icoupl,0,sto5,accm)
         end if
         deallocate(sto5)
      endif

      deallocate(dum)
      call storage_derivatives(0,1)
      itert=itert+iter
      itotals=itotals+iter
      minkt=minkt+mink
 999  continue
      if(islord.ne.0) then
         call switch(nmat,nmatb,islord,2,2,nmatd)
         call switchb(bp,nrhs,nrhsb,islord,2,2,neq)
      end if

c     sort out variable updates
      call air_combine(2)

c     extract dpdp solution
      call dpdp(4)
c
c
c     times=times-second(0.0)
c     write(59,*) 'time to solve jacobian(full)', times

      return
      end


text:
!
!  Written by Leandro Martínez, 2009-2011.
!  Copyright (c) 2009-2011, Leandro Martínez, Jose Mario Martinez,
!  Ernesto G. Birgin.
!

!
! Subroutine that writes the last point obtained when
! a solution was not found
!

subroutine checkpoint(n,x)

  use sizes
  use compute_data
  use input
  use usegencan
  use ahestetic

  implicit none
  integer :: i, strlength
  integer :: n
  double precision :: x(n)
  double precision :: fx
  logical :: movebadprint

  ! All molecules are important

   do i = 1, ntfix
     comptype(i) = .true.
   end do

  ! Call the subroutine that computes de function value

  call computef(n,x,fx)

  write(*,dash3_line)
  write(*,"(&
            &' Packmol was not able to find a solution to your',/,&
            &' packing problem with the desired distance tolerance.',/,/,&
            &' First of all, be sure if the molecules fit in the',/,&
            &' regions specified and if the constraints were set',/,&
            &' correctly. ',/,/,&
            &' Secondly, try simply running it again with a different ',/,&
            &' seed for the random number generator of the initial ',/,&
            &' point. This is done by adding the keyword seed to the',/,&
            &' input file, as in: ',/,/,&
            &' seed 192911 ',/,/,&
            &' The best configuration found has a function value of',/,&
            &' f = ', e14.7,/,/,&
            &' IMPORTANT: ',/,&
            &' If the number of molecules and the restraints are',/,&
            &' correct, it is still very likely that the current point',/,&
            &' fits your needs if your purpose is to run a MD',/,&
            &' simulation.',/,&
            &' Therefore, we recommend to minimize the energy of the',/,&
            &' solution found, equilibrate it and run with it as well.',/&
            &)") fx
  write(*,dash3_line)

  call output(n,x)

  write(*,*) ' The solution with the best function value was '
  write(*,*) ' written to the output file: ', xyzout(1:strlength(xyzout))
  write(*,dash1_line)
  write(*,*) ' Forcing the solution to fit the constraints...'

  ! CALL GENCAN

  init1 = .true.
  do i = 1, nloop
    iprint1 = 0
    iprint2 = 0
    call pgencan(n,x,fx)
    movebadprint = .false.
    call movebad(n,x,fx,movebadprint)
  end do
  init1 = .false.

  write(*,*)
  write(*,dash1_line)
  xyzout = xyzout(1:strlength(xyzout))//'_FORCED'
  call output(n,x)

  write(*,*) ' The forced point was writen to the '
  write(*,*) ' output file: ', xyzout(1:strlength(xyzout)+7)
  write(*,*)
  write(*,*) ' If you want that the packing procedure continues'
  write(*,*) ' for a longer time, add the following keyword '
  write(*,*) ' to the input file: '
  write(*,*)
  write(*,*) ' nloop [integer]      (ex: nloop 200) '
  write(*,*)
  write(*,*) ' The default nloop value is 50 for each molecule.'
  write(*,*)

  write(*,hash1_line)
  write(*,*) ' ENDED WITHOUT PERFECT PACKING: '
  write(*,*) ' The output file:'
  write(*,*)
  write(*,*) '   ',xyzout(1:strlength(xyzout)-7)
  write(*,*)
  write(*,*) ' contains the best solution found. '
  write(*,*)
  write(*,*) ' Very likely, if the input data was correct, '
  write(*,*) ' it is a reasonable starting configuration.'
  write(*,*) ' Check commentaries above for more details. '
  write(*,hash1_line)

  return
end subroutine checkpoint

max_stars_repo_name: susilehtola/packmol


text: #redirect Users/PaulHarms


!//=========================================================================
!// Oslo CTM3 v2015.01
!//=========================================================================
!// Ole Amund Sovde, May 2015
!//=========================================================================
!// DUMMY.
!//=========================================================================
module soa_oslo
  !// ----------------------------------------------------------------------
  !// MODULE: soa_oslo
  !// DESCRIPTION: SOA DUMMY.
  !//
  !// Amund Sovde, July 2012
  !// ----------------------------------------------------------------------
  implicit none
  !// ----------------------------------------------------------------------
  integer,parameter :: ndep_soa = 1
  integer,dimension(ndep_soa),parameter   :: soa_deps = (/-1/)
  !// ----------------------------------------------------------------------
  private
  public soa_init, soa_setdrydep, SOA_v9_separate, soa_diag_drydep, &
       soa_diag_separate, soa_diag2file, soa_nopsdiag, &
       ndep_soa, soa_deps, soa_diag_lsscav, soa_diag_soagas
  !// ----------------------------------------------------------------------

contains

  !// ----------------------------------------------------------------------
  subroutine soa_init()
    !// --------------------------------------------------------------------
    !// Initialize SOA simulations.
    !//
    !// Amund Sovde, June 2012
    !// --------------------------------------------------------------------
    use cmn_size, only: LSOA
    use cmn_chem, only: TNAME
    use cmn_oslo, only: trsp_idx
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------

    !// Check for name of component id 12
    if (.not.LSOA .and. TNAME(trsp_idx(12)).ne.'C6HXR') then
       write(*,'(a)') '* Tracer 12 is wrongly named: '//trim(TNAME(trsp_idx(12)))
       write(*,'(a)') '  It must be named C6HXR when not including SOA'
       stop
    endif

    !// --------------------------------------------------------------------
  end subroutine soa_init
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine soa_setdrydep(VDEP,MP)
    !// --------------------------------------------------------------------
    !// DUMMY.
    !//
    !// Amund Sovde, July 2012
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use cmn_size, only: NPAR, IPAR, JPAR
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: MP
    real(r8), intent(in) :: VDEP(NPAR,IPAR,JPAR)
    !// --------------------------------------------------------------------

    return

    !// --------------------------------------------------------------------
  end subroutine soa_setdrydep
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine SOA_v9_separate( &
       ZC_MASS, & !// the mass of each species
       TEMP, &    !// temperature in column
       DV, &      !// volume in column
       ILOC, &    !// I-box
       JLOC, &    !// J-box
       LSTART, &  !// Start level for calculations
       LEND )     !// End level for calculations
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use cmn_size, only: LPAR, TRACER_ID_MAX
    !// --------------------------------------------------------------------
    !// DUMMY.
    !//
    !// Amund Sovde, July 2012
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: ILOC, JLOC, LSTART, LEND
    real(r8), intent(in) :: TEMP(LPAR), DV(LPAR)
    real(r8), intent(in) :: ZC_MASS(TRACER_ID_MAX,LPAR)
    !// --------------------------------------------------------------------

    return

    !// --------------------------------------------------------------------
  end subroutine SOA_v9_separate
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine soa_diag_drydep(DRYDEP,DV,I,J)
    !// --------------------------------------------------------------------
    !// DUMMY.
    !//
    !// Amund Sovde, July 2012
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use cmn_size, only: LPAR, TRACER_ID_MAX
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: I,J
    real(r8), intent(in) :: DRYDEP(TRACER_ID_MAX), DV(LPAR)
    !// --------------------------------------------------------------------

    return

    !// --------------------------------------------------------------------
  end subroutine soa_diag_drydep
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine soa_diag_separate(PROD,I,J,LSTART,LEND)
    !// --------------------------------------------------------------------
    !// DUMMY.
    !//
    !// Amund Sovde, July 2012
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use cmn_size, only: LPAR, TRACER_ID_MAX
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: I,J,LSTART,LEND
    real(r8), intent(in) :: PROD(TRACER_ID_MAX,LPAR)
    !// --------------------------------------------------------------------

    return

    !// --------------------------------------------------------------------
  end subroutine soa_diag_separate
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine soa_diag2file(NDAY,NDAYI)
    !// --------------------------------------------------------------------
    !// DUMMY.
    !//
    !// Amund Sovde, July 2012
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: NDAY, NDAYI
    !// --------------------------------------------------------------------

    return

    !// --------------------------------------------------------------------
  end subroutine soa_diag2file
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine soa_nopsdiag(NDAY,NMET,NOPS,JDAY,DTOPS)
    !// --------------------------------------------------------------------
    !// DUMMY.
    !//
    !// Amund Sovde, July 2012
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: NDAY, NMET, NOPS, JDAY
    real(r8), intent(in) :: DTOPS
    !// --------------------------------------------------------------------

    return

    !// --------------------------------------------------------------------
  end subroutine soa_nopsdiag
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine soa_diag_lsscav(BTT,BTTBCK,MP)
    !// --------------------------------------------------------------------
    !// DUMMY.
    !//
    !// Amund Sovde, August 2014
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use cmn_size, only: LPAR, NPAR, IDBLK, JDBLK
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: MP
    !// In
    real(r8), dimension(LPAR,NPAR,IDBLK,JDBLK), intent(in):: BTT
    real(r8), dimension(LPAR,NPAR,IDBLK,JDBLK), intent(in):: BTTBCK
    !// --------------------------------------------------------------------

    return

    !// --------------------------------------------------------------------
  end subroutine soa_diag_lsscav
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine soa_diag_soagas(PROD,I,J,LSTART,LEND)
    !// --------------------------------------------------------------------
    !// DUMMY.
    !//
    !// Amund Sovde, May 2016
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use cmn_size, only: TRACER_ID_MAX, LPAR
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: I,J,LSTART,LEND
    real(r8), intent(in) :: PROD(TRACER_ID_MAX,LPAR) !// [kg/gridbox]
    !// --------------------------------------------------------------------

    return

    !// --------------------------------------------------------------------
  end subroutine soa_diag_soagas
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
end module soa_oslo
!//=========================================================================


text:
c
c     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine madsum(lastkn,lastrs,ibegin,aij,
     >                  rslatmd,knlat_x,knlat_y,knlat_z,knlatsq,
     >                  eta,tau,madmat,
     >                  pi4,iprint,istop)
c     ================================================================
c
c     ****************************************************************
c     performs ewald summation for Madelung constant.
c     requires 'lastkn,lastrs,eta,knlat,aij'
c     returns 'madelung sum': madmat
c     ****************************************************************
      implicit   none
c
      character  istop*32
      character  sname*32
      parameter (sname='madsum')
c
      integer    lastkn
      integer    lastrs
      integer    ibegin
      integer    iprint
      integer    i
c
      real*8     knlat_x(lastkn)
      real*8     knlat_y(lastkn)
      real*8     knlat_z(lastkn)
      real*8     knlatsq(lastkn)
      real*8     rslatmd(lastrs)
      real*8     aij(3)
      real*8     eta
      real*8     tau
      real*8     madmat
      real*8     pi4
      real*8     fac
      real*8     zero
      real*8     four
      real*8     erfc
      real*8     erf
c
      parameter (zero=0.0d0)
      parameter (four=4.0d0)
c
c     ================================================================
c                                     2   2       -> ->
c              4*pi          1    -eta *Kq /4 - i*Kq*aij
c     term1 =  ---- * sum  ----- e
c              tau    q<>0    2
c                           Kq
c
c     note sum starts at 2 since kn=0.0 of 1/(rs-aij) is
c     canceled by kn=0.0 of the 1/rs sum.
c     ================================================================
      madmat = zero
      fac=eta*eta/four
      do i = 2,lastkn
         madmat=madmat+exp(-fac*knlatsq(i))/knlatsq(i)
     >                *cos( knlat_x(i)*aij(1)+knlat_y(i)*aij(2)
     >                     +knlat_z(i)*aij(3) )
c if(i.le.5)write(6,'(''i,a,k,f,m'',i5,10f10.4)')
c     >   i,aij(1),aij(2),aij(3),knlat_x(i),knlat_y(i)
c     >   ,knlat_z(i),exp(-fac*knlatsq(i))/knlatsq(i)*pi4/tau
c     >        ,madmat*pi4/tau*6./3.5449077d0
      enddo
      madmat = pi4/tau*madmat
c     ================================================================
c
c                           ->   ->
c                  1 - erf(|Rn + aij|/eta)
c     term2 = sum -------------------------
c              n         ->   ->
c                       |Rn + aij|
c
c     note for calculation of aij=0.0 term ibegin=2.
c     ================================================================
      do i=ibegin,lastrs
c     do i=ibegin,1
	 erf=erfc(rslatmd(i)/eta)/rslatmd(i)
         madmat = madmat + erf
c     if(i.lt.20.and.iprint.ge.0)write(6,'('' i,erf'',i5,3f15.7)')i,
c    >rslatmd(i),eta,erf
      enddo
c
      return
      end

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: hkershaw-brown/MuST