c***********************************************************************
c     dualfoil.f  (version 5.1)  April 7, 2008

c     Dual lithium ion insertion cell
c  Copyright Marc Doyle and John Newman 1998.
c  You may make a copy of this program which you may
c  personally and freely use in its unaltered form.
c  You may distribute this program subject to the
c  conditions that it be made freely available and
c  that any duplication of this program must be
c  essentially unaltered and must include this notice.
c
c  We make no warranties, express or implied, that
c  this program is free of errors or that it will
c  meet the requirements of your application.  The
c  author and publisher disclaim all liablility for
c  direct or consequential damages resulting from
c  use of this program.
c
c  Note: For lflag=0, the model works only for initially zero current.
c
c***********************************************************************
c Based on the original articles:
c Marc Doyle, Thomas F. Fuller, and John Newman,
c J. Electrochem. Soc., 140 (1993), 1526-1533.
c Thomas F. Fuller, Marc Doyle, and John Newman,
c J. Electrochem. Soc., 141 (1994), 1-10.
c Thomas F. Fuller, Marc Doyle, and John Newman,
c J. Electrochem. Soc., 141 (1994), 982-990.

c     This program was developed for lithium and lithium-ion
c     batteries, with the flexibility that a number of positive
c     and negative active materials of the insertion or intercalation
c     type can be selected or added, in addition to a lithium foil
c     electrode.  Several electrolytes are also provided, and
c     these can be added to by providing properties in subroutine
c     prop.  An aqueous system, such as a nickel/metal hydride
c     battery in a KOH solution, could also be added, since it
c     follows the pattern of being a dual-insertion system.  Primary
c     cells, as well as rechargeable cells, can be simulated with
c     this program.

c     This program should be very similar to Marc
c     Doyle's impedance program of 1995 to 1998.
c     See Doyle, Meyers, Darling, Newman, JES, 147 (2000),
c     99-110 and 147 (2000), 2930-2940.
c     The program has been modified to include a
c     zero time step in order to get the potentials
c     and currents correct and to have the impedance
c     independent of the initial values set in
c     subroutine guess.  It should run even with
c     a nonzero current density, but this is not
c     recommended and is not proper in the sense
c     of having a steady state as a basis for the
c     impedance simulation (one of the criteria for
c     the validity of the Kramers-Kronig relations).
c     The program has been cleaned up a lot by
c     rearranging equations ii2, i2div, and kp1 and
c     simplifying the impedance statements.
c
c     Includes FOIL and impedance option 7-1-95
c     Some modifications made to Bellcore stuff on 1-12-96
c     Error in EIS b.c.'s for eqns i2div and ip1 fixed 1-23-96
c     Some additional changes up to 7/96 incl. film resistance
c
c     Modified impedance section by Newman on 4/4/98
c     HALF-CELL version:  split off from impnew1.f on 4/19/98
c     Current newest version June 29, 1998
c     Has Jeremy Meyers modifications 6/98
c
c     Revised June, 1998, to include double-layer capacitance in
c     each electrode and to correct a factor of two in Ohm's law.
c     See J. Electrochem. Soc., 146 (1999), 4360-4365.
c
c     Revised by Karen Thomas, Feb. 12, 1999:
c     - if n1 = 0, then code treats the negative electrode as metal foil.
c     - To run, simply type "dualfoil", then enter
c       input and output file names when prompted.
c     - double-layer capacitance is not currently calculated at
c       a foil electrode
c     July 16, 2000:
c     -  Program modified to be compatible with g77 fortran compiler.
c     July 6, 2001:
c     - Prints total and irreversible heat generation; positive heat
c       is exothermic
c     August 2002:
c     - corrected double precision, added graphite, LiCoO2, and LiNi0.8Co0.2O2
c     May 2003:
c     - made compatible with Compaq visual fortran compiler, other minor changes
c     April, 2005:
c      -The version of Dualfoil dated (July 15, 2003) used a single cutoff
c      potential, which caused problems when a user attempted to charge the cell
c      following a discharge.  For example, if the global cutoff potential was
c      set to 3.0 V and the cell discharged beneath that value, upon charging the
c      cutoff potential would be triggered when the cell potential exceeded 3.0 V.
c      The revised version of dualfoil.f uses two global cutoff potentials,
c      one that serves as an upper bound and one that serves as a lower bound.
c      dualfoil.in has been modified to allow the user to specify two cutoff
c      potentials.
c     July, 2007:
c      -A series of major modifications have been made since the April, 2005 version.
c       See the guide posted on this page for more information.  Note that
c       all previous functionality is maintained by the current version.  Following is
c       an abbreviated list of the changes that have been made:
c        - Addition of specified load and power modes
c        - Ability to specify an upper and lower cutoff potential for each segment
c        - Addition of NiMH materials (MH anode, NiOOH cathode, KOH electrolyte)
c        - Addition of side reactions to model gas-phase recombination in NiMH system
c          or side reactions in another system (e.g., solvent reduction in a Li-ion cell).
c        - Ability to use a variable solid-phase diffusion coefficient in a
c          "pseudo-2d" arrangement
c        - Use of a new input file with activation energies for temperature dependences
c          of the transport properties (solid and electrolyte diffusion, electrolyte
c          conductivity), rate constants, and film resistances.  Where appropriate
c          there is a separate activation energy in each electrode.  The activation
c          energies are stored in the input file ebar.in
c        - Detailed profiles are now written as the program goes along (rather
c          than after the code has finished running) to a separate file, called
c          profiles.out, if the user selects il1=1
c        - A problem with the identity/capacity/density of the cathode in the
c          default input file has been corrected
c        - The unused subroutine util has been removed
c        - An additional equation has been added to improve convergence
c          when running constant potential, power, and load.
c        - The ability to have up to five particles per electrode has been added
c        - The variables and equations have been given names rather than numbers
c        - The energy balance is now done according to the method given
c          in Rao and Newman, J. Electrochem. Soc., 144 (1997), 2697-2704.
c          Cp(dT/dt)-Q=-Integral(sum over reactions (a*in,i*Uh,i))dx-IV
c      February, 2008:
c        The following changes have been made in going from version 5.0 to 5.1
c        - Subroutine guess has been removed
c        - Convergence has been improved and the dV/dt trigger has been removed
c        - fj(j) for cutting off part of the electrode has been removed
c        - The variables have been renamed
c        - A grid resistance, RG, has been added
c        - An error in the impedance calculation for anode and cathode has
c          been corrected; the original full-cell impedance was correct
c        - The equation structure for impedance mode has been made general,
c          facilitating future additions
c        - We have improved the output formatting
c        - We have improved the logic of the multiple-particle selection section

c***********************************************************************
      implicit real*8(a-h,o-z)
      character *30 filin, filout
      parameter(maxt=900)
      common /n/ tmmax,imp,ji,nx,nt,n1,n2,nj,n3,nconv,npa,iSk3,kj3
     &,ii2,ki2,kj,i2div,ip2,kp2,ip1,kp1,imb2,kS2,imb1,kS1,iSk1,kj1
     &,iSk2,kj2
      common/activ/ EbarD,Ebarkap,Ebarka,Ebarkc,Ebarr1,Ebarr3,
     &Ebars1,Ebars3,Ebarks1a,Ebarks1c,Ebarks2a,Ebarks2c,
     &Ebarks3a,Ebarks3c
      common /calc/ ai(maxt,5),ai2(maxt,5),ts(maxt),h,h1,h2,h3,hcn,
     1hcp,rr,rrmax,cuL,sumpa(5,221),Rad1pa(5),Rad3pa(5),area1pa(5)
     &,area3pa(5),mcL
      common/pindiv/pwfindiv(5,221),pwfindiv_d(5,221),vf1(5),vf3(5),izfl
      common/const/ fc,r,t,frt,cur,ep3,ep2,pi,ep1,epf3,epf1,
     &epp1,epp2,epp3,shape3,shape1,capp1,capp3,nneg,nprop,npos
      common/power/ ed,Vold,ranodesave,rcathdesave,heat,qlosstot
      common/ssblock/ xp0(16),xx0(16,221),term(221)
      common/var/ xp(16),xx(17,221),xt(16,221,maxt)
      common/imp/ vreal(110),vimag(110),phi(110),zmag(110),omega(510),
     1vreala(110),vrealc(110),vimaga(110),vimagc(110),
     1phia(110),phic(110),zmaga(110),zmagc(110)
      common/cprop/ sig3,area3,rka3,rka3save,ct3,dfs3,Rad3,cap1,cap3,
     1sig1,area1,rka1,rka1save,ct1,dfs1,Rad1,tw,dfs1save,dfs3save
      common/tprop/df(221),cd(221),tm(221),ddo2(221),ddh2(221),
     1ddf(221),dcd(221),dtm(221),dfu(221),d2fu(221),do2(221),dh2(221)
      common/temp/ thk,htc,dudt,Cp,dens,tam,g0,qq,qloss,residm,ncell,lht
      common/side/rksc1,c1init,c2init,rksa1,term_s1(221),vol,
     &rksa2,rksa3,rksc3,rksc2,UsO2,UsH2,cn2,term_s2(221),nside
      COMMON /vdc/ aa(1,1),bb(1,1),cc(1,150),dd(1,3),gg(1),xxx(1,1),
     1yy(1,1),hha,hhc,cssold(220,150),css(220,150),utz(5,221),
     &dsold(220,150),ds(220,150),time,nn,nnj,np,mvdc1,mvdc3,lims
      common/gas/ epg1,epg2,epg3
      common /maxpow/ lcount
      common/RG/ RG,RGn,RGp,RGext
      dimension terms(221),tt(200),cu(200),tot(200),mc(200),taper(200)
      dimension vcutlo(200),vcuthi(200)
      dimension term_s1s(221)
      dimension rad1paf(5),rad3paf(5)
   44 format(/'            mass  = ',f7.4,' kg/m2')
   45 format(' specific energy whole run  = ',f8.2,' W-h/kg')
   46 format(' specific power whole run  = ',f8.2,' W/kg')
   47 format(' total heat = ', f15.4,' W-h/m2')
   48 format(' specific energy segment = ', f8.2,' W-h/kg')
   49 format(' specific power segment  = ', f8.2,' W/kg')
   57 format(e15.6,', ',e15.6,', ',e15.6)
   89 format(f8.1,',',f8.1,',',f8.1,',',i3,',',g10.5,',',f5.2)
   90 format (g10.3,',')

c
      open(3,file='halfcells.out',status='unknown')
      write (3,*) ' Time (min)    V neg      V pos       Cur (A/m2)   '
     &,' Temp (C)   Heat Gen (W/m3)'
c     print *, 'Enter input file name, press return'
c     read *, filin
      open (1, FILE = 'dualfoil5.in', status = 'old')
c     open (1, FILE = filin, status = 'old')
c     print *, 'Enter output file name, press return'
c     read *, filout
      open (2, file = 'dualfoil5.out', status = 'unknown')
c     open (2, file = filout, status = 'unknown')
      open (4,file='profiles.out',status='unknown') ! for detailed profiles
      open (5,file='li-ion-ebar.in', status='old') ! activation energies
      open (6,file='resistances.out',status='unknown') ! for cell resistance
      open (7,file='solidprof.out',status='unknown') ! for profile inside solid
c
      data fc/96487.0d0/, r/8.314d0/, pi/3.141592653589d0/
      data ed/0.d0/, Vold/0.d0/

      npa=1 ! Number of particle sizes
      if (npa.ne.1.and.imp.eq.1) print *,'npa must be 1 for imp=1'

      ii2=3+npa ! equation number for current density i2
      ki2=3+npa ! variable number for current density i2
      i2div=4+npa ! equation number for current balance.
      kj=4+npa ! variable number for transfer current (flux density).
      ip1=5+npa ! equation number for Ohm's law in the matrix.
      kp1=5+npa ! variable number for PHI_1, matrix potential.
      ip2=2 ! equation number for Ohm's law in solution.
      kp2=2 ! variable number for PHI_2, solution potential.
      imb1=6+npa ! equation number (side reaction 1 matl balance)
      kS1=6+npa ! variable number (side reaction 1 matl balance)
      iSk1=7+npa ! equation number (side reaction 1 kinetics)
      kj1=7+npa ! variable number (side reaction 1 kinetics)
      iSk2=8+npa ! equation number (side reaction 2 kinetics)
      kj2=8+npa ! variable number (side reaction 2 kinetics)
      imb2=9+npa ! equation number (side reaction 2 matl balance)
      kS2=9+npa ! variable number (side reaction 2 matl balance)
      iSk3=10+npa ! equation number (side reaction 3 kinetics)
      kj3=10+npa ! variable number (side reaction 3 kinetics)
      lim2=100  ! limit for iterations to converge load/power/poten

c     Split of RG into its components. WHT
      RG=5.0d-3 !total resistance in foils, leads, and contacts, ohm-m2 WHT
      RG=0.0d0
      RGn=RG/3.0d0 ! resistance affecting negative half cell
      RGp=RG-RGn ! resistance affecting positive half cell
      RGext=0.0d0 ! resistance outside cell
c
c***********************************************************************
c     read in parameters and boundary conditions

c
      read (1,*) lim !limit on number of iterations
      read (1,*) h1  !thickness of negative electrode (m)
      read (1,*) h2  !thickness of separator (m)
      read (1,*) h3  !thickness of positive electrode (m)
      read (1,*) hcn !thickness of negative electrode current collector (m)
      read (1,*) hcp !thickness of positive electrode current collector (m)
      thk=h1+h2+h3
      read (1,*) n1  !number of nodes in negative electrode
c     If negative electrode is metal foil, let n1 = 0
      read (1,*) n2  !number of nodes in separator
      read (1,*) n3  !number of nodes in positive electrode
      read (1,*) n4  !number of nodes in solid particle
      read (1,*) mvdc1 !flag for variable solid diff coeff, anode
      read (1,*) mvdc3 !flag for variable solid diff coeff, cathode
      read (1,*) lims  !limit on number of iterations in solid phase
      read (1,*) t   !temperature (K)
      write (2,1101) lim,1.d6*h1,1.d6*h2,1.d6*h3,1.d6*hcn,1.d6*hcp
     &,n1,n2,n3,t
      n2=n2+1
      nj=n1+n2+n3
      nnj=n4
c
      read (1,*) xx(1,n1+2) ! initial concentration (mol/m3)
      read (1,*) csx  !initial stochiometric parameter for negative
      read (1,*) csy  !initial stochiometric parameter for positive
      read (1,*) tmmax !maximum time step size (s)
      read (1,*) dfs1 !diffusion coefficient in negative solid (m2/s)
      read (1,*) dfs3 !diffusion coefficient in positive solid (m2/s)
      dfs1save=dfs1 !save the original values because of temp dependence
      dfs3save=dfs3
      read (1,*) Rad1 !radius of negative particles (m)
      read (1,*) Rad3 !radius of positive particles (m)
      write (2,1102) xx(1,n1+2),csx,csy,tmmax,dfs1,dfs3,
     &1.d6*Rad1,1.d6*Rad3
c     If negative electrode is metal foil, let ep1=epp1=epf1=0.0
      read (1,*) ep1  !volume fraction of electrolyte in negative electrode
      read (1,*) epp1 !volume fraction of polymer phase in negative electrode
      read (1,*) epf1 !volume fraction of inert filler in negative electrode
      read (1,*) epg1 !volume fraction of gas in negative
      read (1,*) ep2  !volume fraction of electrolyte in separator
      read (1,*) epp2 !volume fraction of polymer phase in separator
      read (1,*) epg2 !volume fraction gas in separator
      read (1,*) ep3  !volume fraction of electrolyte in positive electrode
      read (1,*) epp3 !volume fraction of polymer phase in positive electrode
      read (1,*) epf3 !volume fraction of inert filler in positive electrode
      read (1,*) epg3 !volume fraction of gas in positive
      read (1,*) sig1 !conductivity of solid negative matrix (S/m)
      read (1,*) sig3 !conductivity of solid positive matrix (S/m)
      read (1,*) rka1 !reaction rate constant for negative reaction
      read (1,*) rka3 !reaction rate constant for positive reaction
      rka1save=rka1
      rka3save=rka3
      read (1,*) ranode  !anode film resistance (out of place)
      read (1,*) rcathde !cathode film resistance (out of place)
      ranodesave=ranode
      rcathdesave=rcathde
      read (1,*) cot1 !coulombic capacity of negative material (mAh/g)
      read (1,*) cot3 !coulombic capacity of positive material (mAh/g)
      write (2,1103) ep1,epp1,epf1,ep2,epp2,ep3,epp3,epf3,sig1,
     & sig3,cot1, cot3,rka1,rka3
      read (1,*) re  ! density of electrolyte (kg/m3)
      read (1,*) rs1 ! density of negative insertion material (kg/m3)
      read (1,*) rs3 ! density of positive insertion material (kg/m3)
      read (1,*) rf  ! density of inert filler (kg/m3)
      read (1,*) rpl ! density of polymer phase (kg/m3)
      read (1,*) rc  ! density of separator material (kg/m3)
      read (1,*) rcn ! density of negative current collector (kg/m3)
      read (1,*) rcp ! density of positive current collector (kg/m3)
      write (2,1104) re,rs1,rs3,rf,rpl,rc,rcn,rcp
      read (1,*) htc  !heat-transfer coefficient with external medium (W/m2K)
      read (1,*) Cp   !heat capacity of cell (J/kgK)
      read (1,*) Tam  !ambient temperature (K)
      read (1,*) ncell !number of cells in a cell stack
      read (1,*) lht  !0 uses htc, 1 calcs htc, 2 isothermal
      read (1,*) il1   !1 for long print-out 0 for short print-out
      read (1,*) il2   !1/il2 = fraction of nodes in long print-out
      read (1,*) il3   !1/il3 = fraction of time steps in long print-out
      read (1,*) lflag ! 0 for electrolyte in separator only, 1 for uniform
      read (1,*) imp  ! 0 for no impedance, 1 for impedance
      read (1,*) capp1 ! capacitance of negative material (F/m2)
      read (1,*) capp3 ! capacitance of positive material (F/m2)
      read (1,*) lpow  ! 0 for no power peaks, 1 for power peaks
      read (1,*) jsol  ! calculate solid profiles if 1<jsol<nj
      read (1,*) nside ! flag for side reaction
      read (1,*) rksa1 ! rate constant side reaction anode
      read (1,*) rksc1 ! rate constant side reaction cathode
      read (1,*) rksa2 ! rate constant side reaction anode
      read (1,*) rksc2 ! rate constant side reaction cathode
      read (1,*) rksa3 ! rate constant side reaction anode
      read (1,*) rksc3 ! rate constant side reaction cathode
      c1init=0.1d0   ! init conc of side reaction species 1 (mol/m^3)
      c2init=3.0d0     ! init conc of side reaction species 2 (mol/m^3)
      cn2=45.0d0       ! background N2 concentration for NiMH (mol/m^3)
      residm=0.0d0     ! Residual mass not in cell sandwich (kg/m^2)
      write (2,1105) ranode,rcathde,htc,dudt,Cp,residm,tam,ncell,lht
      write (2,1116) nside
      if (nside.ge.1) then
      write (2,1111) rksa1,rksc1,rksa2,rksc2,rksa3,rksc3
      endif
      read (1,*) nneg  ! designates negative electrode system
      read (1,*) nprop ! designates electrolyte system
      read (1,*) npos  ! designates positive electrode system
      read (1,*) lcurs ! number of current changes
      write (2,1106) il1,il2,il3,lflag,imp,capp1,capp3,lcurs
      read (5,*) EbarS1  !activation energy, anode solid diffusion
      read (5,*) EbarS3  !activiation energy, cathode solid diffusion
      read (5,*) Ebarkap !activation energy electrolyte conductivity
      read (5,*) EbarD   !activation energy electrolyte diffusion
      read (5,*) Ebarka  !activation energy negative kinetics
      read (5,*) Ebarkc  !activation energy positive kinetics
      read (5,*) Ebarks1a  !activation energy O2 side rxn
      read (5,*) Ebarks1c  !activation energy O2 side rxn
      read (5,*) Ebarks2a  !activation energy H2 side rxn
      read (5,*) Ebarks2c  !activation energy H2 side rxn
      read (5,*) Ebarks3a  !activation energy shuttle side rxn
      read (5,*) Ebarks3c  !activation energy shuttle side rxn
      read (5,*) Ebarr1  !activation energy, film resistance anode
      read (5,*) Ebarr3  !activation energy, film resistance cathode

      nn=1 ! used for variable solid-phase diffusion coefficient.
c     n is number of equations
      if (nside.eq.0) n=5+npa  !no side reaction, use six eqns
      if (nside.eq.1) n=7+npa  ! first side reaction, mass balance & kinetic
      if (nside.eq.2) n=9+npa  ! second side reaction, mass balance & kinetic
      if (nside.eq.3) n=10+npa  ! third side reaction, only kinetic

      if (nside.eq.1) kj2=kj2+1 !to prevent variable overlap
c
      do 555 i = 1, lcurs
      read (1,*) cu(i),tt(i),mc(i),vcutlo(i),vcuthi(i)
  555 continue
c     cu(i)      operating current density (A/m2) or cell potential (V)
c     tt(i)      time (min) or cutoff potential (V)
c     mc(i) The mode of discharge
c           0 for potentiostatic
c           1 for galvanostatic for a given time
c           2 for galvanostatic to a cutoff potential
c          -1 for galvanostatic for a given time with tapered
c             charge/discharge upon reaching cutoff potential
c          -2 for specified power (in W/m2)
c          -3 for specified load (in ohm-m2)
c     vcutlo(i)  lower cutoff potential
c     vcuthi(i)  upper cutoff potential

      print *, 'Now running DUAL...'
 1101 format (i7,'  lim, limit on number of iterations'
     &/1x,f6.2,'  h1,  thickness of negative electrode (microns)'
     &/1x,f6.2,'  h2,  thickness of separator (microns)'
     &/1x,f6.2,'  h3,  thickness of positive electrode (microns)'
     &/1x,f6.2,'  hcn, ',
     &'thickness of negative electrode current collector (microns)'
     &/1x,f6.2,'  hcp, thickness of positive electrode current'
     &,' collector (microns)'
     &/i7,'  n1,  number of nodes in negative electrode'
     &/i7,'  n2,  number of nodes in separator'
     &/i7,'  n3,  number of nodes in positive electrode'
     &/1x,f6.2, '  T,   temperature (K)')
 1102 format (/1x,f6.1,'  xx(1,n1+2), initial concentration (mol/m3)'
     &/1x,f6.4,'  csx,   initial stoichiometric parameter for negative'
     &/1x,f6.4,'  csy,   initial stoichiometric parameter for positive'
     &/1x,f6.1,'  tmmax, maximum time step size (s)'
     &/1x,e6.1,'  dfs1, diffusion coefficient in negative solid (m2/s)'
     &/1x,e6.1,'  dfs3, diffusion coefficient in positive solid (m2/s)'
     &/1x,f6.2,'  Rad1,  radius of negative particles (microns)'
     &/1x,f6.2,'  Rad3,  radius of positive particles (microns)')
 1103 format (/1x,f6.3,'  ep1,'
     &,'   volume fraction of electrolyte in negative electrode'
     &/1x,f6.3,'  epp1,'
     &,'  volume fraction of polymer phase in negative electrode'
     &/1x,f6.3,'  epf1,'
     &,'  volume fraction of inert filler in negative electrode'
     &/1x,f6.3,'  ep2,   volume fraction of electrolyte in separator'
     &/1x,f6.3,'  epp2,  volume fraction of polymer phase in separator'
     &/1x,f6.3,'  ep3,'
     &,'   volume fraction of electrolyte in positive electrode'
     &/1x,f6.3,'  epp3,'
     &,'  volume fraction of polymer phase in positive electrode'
     &/1x,f6.3,'  epf3,'
     &,'  volume fraction of inert filler in positive electrode'
     &/1x,g8.2,'  sig1,  conductivity of negative matrix (S/m)'
     &/1x,g8.2,'  sig3,  conductivity of positive matrix (S/m)'
     &/1x,f7.2,'  cot1,  coulombic capacity of negative material'
     &,' (mAh/g)'
     &/1x,f7.2,'  cot3,  coulombic capacity of positive material'
     &,' (mAh/g)'
     &/1x,e6.1,'  rka1,  reaction rate constant for negative reaction '
     &/1x,e6.1,'  rka3,  reaction rate constant for positive reaction')
 1104 format (/1x,f6.1,'  re,   density of electrolyte (kg/m3)'
     &/1x,f6.1,'  rs1,  density of negative insertion material (kg/m3)'
     &/1x,f6.1,'  rs3,  density of positive insertion material (kg/m3)'
     &/1x,f6.1,'  rf,   density of inert filler (kg/m3)'
     &/1x,f6.1,'  rpl,  density of polymer phase (kg/m3)'
     &/1x,f6.1,'  rc,   density of separator material (kg/m3)'
     &/1x,f6.1,'  rcn,  density of negative current collector (kg/m3)'
     &/1x,f6.1,'  rcp,  density of positive current collector (kg/m3)')
 1105 format (/1x,f10.6,'  ranode,   anode film resistance (ohm-m2)'
     &/1x,f10.6,'  rcathde,  cathode film resistance (ohm-m2)'
     &/1x,f6.2,'  htc,   heat-transfer coefficient with'
     &,' external medium (W/m2K)'
     &/1x,f10.6,'  dUdT,  temperature coefficient of EMF (V/K)'
     &/1x,f6.1,'  Cp,    heat capacity of cell (J/kg-K)'
     &/1x,f6.2,'  residm,    residual mass (kg/m2)'
     &/1x,f6.2,'  Tam,   ambient temperature (K)'
     &/i7,'  ncell, number of cells in a cell stack'
     &/i7,'  lht,   0 uses htc,  1 calcs htc,  2 isothermal')
 1106 format (/i7,'  il1,   1 for long print-out  0 for short print-out'
     &/i7,'  il2,   prints every il2 th node in long print-out'
     &/i7,'  il3,   prints every il3 th time step in long print-out'
     &/i7,'  lflag, 0 for electrolyte in separator only, 1 for uniform'
     &/i7,'  imp, 0 for no impedance, 1 for impedance'
     &/1x,f6.2,'  capp1,  capacitance of negative material'
     &,' (F/m2)'
     &/1x,f6.2,'  capp3,  capacitance of positive material'
     &,' (F/m2)'
     &/i7,'  lcurs, number of current changes')
 1107 format ('   Time     Util N  Util P  Cell Pot   Uocp      Curr',
     &'      pH2     pO2   Total P  Temp   heatgen')
 1108 format ('   (min)       x       y      (V)       (V)      (A/m2)',
     %'   (bar)   (bar)   (bar)   (C)     (W/m2)')
 1109  format ('   Time     Util N  Util P  Cell Pot   Uocp      Curr',
     &'      Temp   heatgen')
 1110 format ('   (min)       x       y      (V)       (V)      (A/m2)',
     %'    (C)    (W/m2)')
 1111 format (
     &g10.3,'  rksa1,   rate constant 1 for negative side reaction'
     &/1x,g10.3,'  rksc1,   rate constant 1 for positive side reaction'
     &/1x,g10.3,'  rksa2,   rate constant 2 for negative side reaction'
     &/1x,g10.3,'  rksc2,   rate constant 2 for positive side reaction'
     &/1x,g10.3,'  rksa3,  rate constant 3 for negative side reaction'
     &/1x,g10.3,'  rksc3r,  rate constant 3 for positive side reaction')
 1114 format ('   Time     Util N  Util P  Cell Pot   Uocp      Curr',
     &'    Temp    Heat Gen')
 1115  format ('   (min)       x       y      (V)       (V)      (A/m2)',
     %'  (C)    (W/m2)')
 1116  format(/i7, '    nside, side reaction flag ')
      write (2,*) ' '

      go to (131,132,133,134,135,136,137,138),nneg
  131 write (2,*) 'Li foil'
      go to 147
  132 write (2,*) 'Carbon (petroleum coke)'
      go to 147
  133 write (2,*) 'MCMB 2528 Graphite (Bellcore)'
      go to 147
  134 write (2,*) 'TiS2'
      go to 147
  135 write (2,*) 'Tungsten oxide (LixWO3 with 0<x<0.67)'
      go to 147
  136 write (2,*) 'Lonza KS6 graphite (Bellcore)'
      go to 147
  137 write (2,*) 'Metal Hydride'
      go to 147
  138 write (2,*) 'Add your own negative electrode'
  147 go to (101,102,103,104,105,106,107,108,109,110,111,112,
     & 113,114),nprop
  101 write (2,*) 'AsF6 in methyl acetate'
      go to 200
  102 write (2,*) 'Perchlorate in PEO'
      go to 200
  103 write (2,*) 'Sodium Triflate in PEO'
      go to 200
  104 write (2,*) 'LiPF6 in PC (Sony cell simulation)'
      go to 200
  105 write (2,*) 'Perchlorate in PC (West simulation)'
      go to 200
  106 write (2,*) 'Triflate in PEO'
      go to 200
  107 write (2,*) 'LiPF6 in EC/DMC and p(VdF-HFP) (Bellcore)'
      go to 200
  108 write (2,*) 'LiPF6 in EC/DMC and p(VdF-HFP) (Bellcore) cell #2'
      go to 200
  109 write (2,*) 'Ion exchange membrane, t+=1.0'
      go to 200
  110 write (2,*) 'LiTFSI in PEMO (data from Steve Sloop)'
      go to 200
  111 write (2,*) 'LiPF6 in EC:DMC'
      go to 200
  112 write (2,*) 'LiTFSI in PEO at 85 C (data from Ludvig Edman)'
      go to 200
  113 write (2,*) 'KOH in H2O'
      go to 200
  114 write (2,*) 'add your own electrolyte'
  200 go to (201,202,203,204,205,206,207,208,209,210,211,212,213),npos
  201 write (2,*) 'TiS2'
      go to 300
  202 write (2,*) 'Spinel Mn2O4 (lower plateau)'
      go to 300
  203 write (2,*) 'NaCoO2:  Sodium Cobalt Oxide'
      go to 300
  204 write (2,*) 'Spinel Mn2O4 (upper plateau)'
      go to 300
  205 write (2,*) 'Tungsten oxide (LixWO3 with 0<x<0.67)'
      go to 300
  206 write (2,*) 'LiCoO2 (Cobalt dioxide)'
      go to 300
  207 write (2,*) 'V2O5 (Vanadium oxide)'
      go to 300
  208 write (2,*) 'LiNi0.8Co0.2O2 (Nickel dioxide Gen 1)'
      go to 300
  209 write (2,*) 'Spinel Mn2O4 (Bellcore)'
      go to 300
  210 write (2,*) 'V6O13 (Vanadium oxide)'
      go to 300
  211 write (2,*) 'LiAl0.2Mn1.8O4F0.2 spinel from Bellcore'
      go to 300
  212 write (2,*) 'NiOOH'
      go to 300
  213 write (2,*) 'Add your own'
  300 continue

      if (mvdc1.eq.1) write (2,*)
     &'Variable solid-phase diff coeff, anode'
      if (mvdc3.eq.1) write (2,*)
     &'Variable solid-phase diff coeff, cathode'

c     Convert coulombic capacity to total concentrations in ekin
      ct1=3.6d03*cot1*rs1/fc
      ct3=3.6d03*cot3*rs3/fc
c      print *,'capacity ratio',(h1*ct1*(1-ep1-epp1-epf1-epg1))/
c     &(h3*ct3*(1-ep3-epp3-epf3-epg3))
c      print *,'Ah/m2, positive',(h3*cot3*(1-ep3-epp3-epf3-epg3)*rs3)
c Calculate the amount of gas in the headspace (1.4 factor for 2005 Prius battery)
c This is necessary only for a gas-phase material balance
      vol=(epg1*h1+epg2*h2+epg3*h3)/((epg1*h1+epg2*h2+epg3*h3)+
     &0.5d0*(h1+h2+h3))/1.4d0
      if (nside.ge.1.and.nprop.ne.13) vol=1.0d0 !vol=1 for liquid-phase side rxn
      lcount=1 !counter for max power

      shape1=3.0d0 ! 3. for spherical, 1. for planar, 2. for cylindrical
      shape3=3.0d0
      if (imp.eq.0) then 
      cap1=capp1 ! F/m2, capacitance for negative
      cap3=capp3 ! F/m2, capacitance for positive
      else ! For impedance set cap1 and cap3 to zero.
      cap1= 0.0d0 ! F/m2, capacitance for negative
      cap3= 0.0d0 ! F/m2, capacitance for positive
      endif

c     Keeps track of tapered discharge/charge mode
      do i=1,lcurs
         taper(i) = 0
         if(mc(i).eq.-1) then
            mc(i) = 1
            taper(i) = 1
         endif
      enddo
c
c     Convert times to seconds and sum up times of mode changes
      if (mc(1).lt.2) then
      tot(1)=6.0d01*tt(1)
      else
      tot(1)=0.0d0
      end if
      do 51 i=2,lcurs
      if (mc(i).lt.2) then
      tot(i)=tot(i-1)+6.0d01*tt(i)
      else
      tot(i)=tot(i-1)
      end if
   51 continue
c
c     ** This section defines the particle size and distribution **
c
c     Single particle in each electrode
c     Specific area calculated from geometry
      area3=shape3*(1.0d0-ep3-epf3-epp3-epg3)/Rad3
      if (n1 .gt. 0) then
      area1=shape1*(1.0d0-ep1-epf1-epp1-epg1)/Rad1
      else
      area1 = 1.0d0/h1
      endif

c Set up the radius and area of each particle.  The following six lines
c apply to electrodes with a single particle, or electrodes with
c multiple particles all of the same radius and area.
      if (npa.eq.1) then
      do mpa=1,npa
      Rad1pa(mpa)=Rad1 ! start with uniform particles
      area1pa(mpa)=area1/dble(npa)
      Rad3pa(mpa)=Rad3 ! start with uniform particles
      area3pa(mpa)=area3/dble(npa)
      enddo ! mpa
      vf1(1)=1.0d0
      vf3(1)=1.0d0
      else

c For nonuniform particles use this section.
c We use volume fractions to specify our distribution; the user should
c be able to calculate how to turn other distributions, such as number
c or area fraction, into volume fractions.

      vol1=(1.0d0-ep1-epf1-epp1-epg1) ! Solid vol frac in negative
      vol3=(1.0d0-ep3-epf3-epp3-epg3) ! Solid vol frac in positive

c Here input the radii of each particle.  This is set up for a two-particle
c distribution, but can be extended to as many as five particles by
c adding additional particle radii here.

      Rad1pa(1)=1.0d-6
      Rad1pa(2)=10.0d-6
      Rad1pa(3)=1.0d-6
      Rad1pa(4)=1.0d-6
      Rad1pa(5)=1.0d-6

      Rad3pa(1)=1.0d-6
      Rad3pa(2)=10.0d-6
      Rad3pa(3)=1.0d-6
      Rad3pa(4)=1.0d-6
      Rad3pa(5)=1.0d-6

c Here input the relative volume fraction of each particle.  Note that
c these need to add to 1.0, not to vol1 or vol3!  For example, a
c relative volume fraction of 0.2 means that 20% of the solid volume of an
c electrode is composed of particles of a given radius.

      vf1(1)=0.2d0
      vf1(2)=0.8d0
      vf1(3)=0.2d0
      vf1(4)=0.2d0
      vf1(5)=0.2d0

      vf3(1)=0.2d0
      vf3(2)=0.8d0
      vf3(3)=0.2d0
      vf3(4)=0.2d0
      vf3(5)=0.2d0

c Reset area1 and area3 for nonuniform particles.
      area1=0.0d0
      area3=0.0d0
c
      do mpa=1,npa
      area1pa(mpa)=vol1*vf1(mpa)/(4.d0/3.d0*pi*Rad1pa(mpa)**3.0d0)*
     &4.d0*pi*Rad1pa(mpa)**2.d0
      area3pa(mpa)=vol3*vf3(mpa)/(4.d0/3.d0*pi*Rad3pa(mpa)**3.0d0)*
     &4.d0*pi*Rad3pa(mpa)**2.d0
      area1=area1+area1pa(mpa)
      area3=area3+area3pa(mpa)
      enddo !mpa
      endif

      sig3=sig3*((1.0d0-ep3-epp3-epg3)**1.5d0)
      sig1=sig1*((1.0d0-ep1-epp1-epg3)**1.5d0)
c
      h2=h2/dble(n2-1)
      h3=h3/dble(n3)
      if (n1 .gt. 1) h1=h1/dble(n1)
      h=h2
      hha=rad1/dble(n4-1)
      hhc=rad3/dble(n4-1)
      frt=fc/(r*t)
c
      mpa=1 ! only one size for initialization
c     Find initial solid-phase potential guesses
c     from initial solid concentrations:
      call ekin(1,1,1,csx,mpa,mvdc1,mvdc3)
      write (2,*) 'open-circuit potential, negative ',g0
      write (2,*) 'delta S negative ', dudt
      Ua=g0
      call ekin(nj,1,1,csy,mpa,mvdc1,mvdc3)
      write (2,*) 'open-circuit potential, positive ',g0
      write (2,*) 'delta S positive ', dudt
      Uc=g0
      OCP=Uc-Ua ! open-circuit cell potential

      write (2,*) ' '
      write (2,*) '     DUAL INSERTION CELL VERSION 5.1'
      write (2,*) ' '
      if (imp.eq.1) then
      write (2,*) '         IMPEDANCE OPTION   '
      write (2,*) ' '
      write (2,*) '  Omega                Z real         Z imag  '
      write (2,*) ' (Rad/s)              (ohm cm2)      (ohm cm2)  '
      go to 54
      endif ! end of impedance section (output)

        if (nside.ne.0.and.nneg.eq.7) then
        write (2,1107) !This output includes pressures
        write (2,1108)
        else
        write (2,1109) !This output excludes pressures
        write (2,1110)
        endif
      write (2,*) ' '

      write (6,*) 'Resistances in each portion of the cell, from a',
     1' simple (delta V) / I'
      write (6,*) 'Time       Cur Den      R Anode       R Separator'
     &,'      R Cathode        R Total'
      write (6,*) '(min)      (A/m^2)     (ohm-m^2)       (ohm-m^2)'
     &,'       (ohm-m^2)      (ohm-m^2)'
      write (6,*) ' '
   54 continue
c
c***********************************************************************
c
c     initialize time counting variables
      k=1
      time=0.0d0
      time2=0.0d0
      rr=0.0d0
      coul=0.d0
      L=1
      ts(1)=0.0d0

      do j=1,nj
        xx(1,j)=xx(1,n1+2) ! initial concentration (mol/m3)
c     Uniform initial concentration if lflag=1.
c     Step-function initial concentration if lflag=0.
        if(lflag.eq.0 .and. (j.le.n1+1 .or. j.ge.n1+n2))
     &  xx(1,j)=1.0d-01
      enddo

c     assume current density linear in electrodes
      call prop(nj,n2,n1)
      curold=0.d0 ! open circuit
      vvold=OCP
      cur=1.0d0 ! guess current to get internal resistance
c     cur=1000.d0 ! If high currents are expected, use a high current
c                  here to get an appropriate resistance
      if(mc(1).ge.1.or.mc(1).eq.-1) cur=cu(1)
      if (cap1.ne.0.0d0.or.cap3.ne.0.0d0) cur=0.0d0

      do j=1,n1+1 ! negative electrode
        do mpa=1,npa
          xx(2+mpa,j)=ct1*csx
          utz(mpa,j)=csx
        enddo ! mpa
      enddo
      do j=n1+2,n1+n2-1 ! separator
        do mpa=1,npa
          xx(2+mpa,j)=0.d0
          utz(mpa,j)=0.d0
        enddo ! mpa
      enddo
      do j=n1+n2,nj ! positive electrode
        do mpa=1,npa
          xx(2+mpa,j)=ct3*csy
          utz(mpa,j)=csy
        enddo ! mpa
      enddo

c     n1+1 is the last node in the negative
c     n1+n2 is the first node in the positive
c     etaneg=xx(kj,j)*fc*ranode
c     etapos=xx(kj,j)*fc*rcathde

      xx(kp2,nj)=0.d0
      do 78 j=nj,n2+n1,-1 ! positive
        xx(ki2,j)=cur*dble(nj-j)/dble(n3)
        xx(kj,j)=-cur/fc/h3/area3/dble(n3)
        if(j.ne.nj)
     &  xx(kp2,j)=xx(kp2,j+1)+h3*xx(ki2,j)/cd(j) ! ohmic drop is not accurate
   78 xx(kp1,j)=Uc+xx(kp2,j)+xx(kj,j)*fc*rcathde

      do j=n1+n2-1,n1+2,-1 ! separator
        xx(ki2,j)=cur
        xx(kj,j)=0.d0
        xx(kp2,j)=xx(kp2,j+1)+h2*xx(ki2,j)/cd(j) ! ohmic drop is not accurate
        xx(kp1,j)=0.0d0
      enddo

      do 79 j=n1+1,1,-1 ! negative electrode
          if (n1.gt.0) then
          xx(ki2,j)=cur*dble(j-1)/dble(n1)
          xx(kj,j)=cur/fc/h1/area1/dble(n1)
          xx(kp2,j)=xx(kp2,j+1)+h1*xx(ki2,j)/cd(j) ! ohmic drop is not accurate
          xx(kp1,j)=xx(kp1,j)+xx(kp2,j)+xx(kj,j)*fc*ranode
          else
          xx(ki2,j)=cur*dble(j-1)/dble(n1+1)
          xx(kj,j)=cur/fc/h1/area1/dble(n1+1)
          xx(kp2,j)=xx(kp2,j+1)+h1*xx(ki2,j)/cd(j) ! ohmic drop is not accurate
          xx(kp1,j)=Ua+xx(kp2,j)+xx(kj,j)*fc*ranode
          endif
   79 continue
c
c initialization for variables kS1,kj1,kS2,kj3 ! side reactions, not impedance
      do j=1,nj
      if (nside.ge.1) then
      xx(kS1,j)=c1init
      xx(kS2,j)=c2init
      endif
      xx(kj1,j)=0.0d0
      xx(kj2,j)=0.0d0
      xx(kj3,j)=0.0d0
      enddo

      do j=1,nj
      do i=1,n
      xt(i,j,1)=xx(i,j)
      enddo
      enddo

c Assumed temperature dependence for solid-state diffusion
      dfs1=dfs1save*dexp((EbarS1)*(t-298.0d0)/(t*298.0d0))
      dfs3=dfs3save*dexp((EbarS3)*(t-298.0d0)/(t*298.0d0))

c Initialize solid-phase concentrations. PSA
      if (mvdc1.eq.1.or.mvdc3.eq.3) then
      mpa=1 ! only one size for variable solid-phase diffusion coefficient
        do jj=1,nj
          do jjj=1,nnj
            css(jj,jjj)=xx(2+mpa,jj)
            cssold(jj,jjj)=xx(2+mpa,jj)
c For a true variable solid-phase diffusion coefficient put in the
c functional dependence here and in subroutine vardc
            if (jj.le.n1+1) ds(jj,jjj)=dfs1
            if (jj.ge.n1+n2) ds(jj,jjj)=dfs3
          enddo
        enddo
      endif

      mcL=1
      if (imp.eq.0) call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 1 not converged' ! initial
         if(nconv.eq.1) go to 533
      vv=xx(kp1,nj)-xx(kp1,1)
      Rint=-(vv-vvold)/(cur-curold)+RG ! internal resistance
      OCP=vv+(Rint-RG)*cur
      pwrmax=OCP**2/Rint/4.d0
      if (pwrmax.le.cu(L) .and. L.eq.-2) then
      write (2,*) 'WARNING: Power is above estimated power maximum'
      write (2,*) 'Requested power= ',cu(L)
      write (2,*) 'Power maximum= ',pwrmax
      write (2,*) 'Lower the power or change the initial current guess.'
      endif

      cuL=cu(L)
      cur=cu(L) ! good for mc(L)=1, 2, and -1.
c     First guess for current for constant load.
      if(mc(L).eq.-3) cur=OCP/(cu(L)+Rint)
      if(mc(L).eq.0) cur=(OCP-cu(L))/Rint ! constant potential
      if(mc(L).eq.-2) then ! constant power
        cur=cu(L)/OCP
        do i=1,6 ! refine estimate of current for given power
          cur=cur-(cur*OCP-Rint*cur**2-cu(L))/(OCP-2.d0*Rint*cur)
        enddo
      endif
      
	avgeta=0.d0
      if(mc(L).eq.0) avgeta=(OCP-cuL)/2.d0
      if(mc(L).le.-2) avgeta=cur*(Rint-RG)/2.d0

c     Must activate lpow=1 in data file if you want peak powers:
      kkflag=0
      nflag=0

c     calculate mass (kg/m2) of the cell
      call mass(re,rs3,rs1,rf,rpl,rc,rcn,rcp)
      dens=tw/thk

c     Impedance section:
      if (imp.eq.1) then
c     initialize steady state to make impedance runs.
      imp=0
      call comp(n,lim,k,rr,0,nflag,0,jcount)

      imp=1
      ji=1
c     omega is the array of perturbation frequencies in Rad/s
      do 61 iii=1, 110

      omega(ji)=10.d6/(10.d0**0.1d0)**iii
c
      call comp(n,lim,k,rr,0,nflag,0,jcount)
c
   61 continue
      write (2,*)''
      write (2,*)'  Omega            Z real        Z imag  '
      write (2,*)' (Rad/s)          (ohmcm2)      (ohmcm2) '
      do 62 i=1,110
   62 write (2,57)omega(i),vreal(i)*10000,vimag(i)*10000
      write (2,*)''
      write (2,*)'  Omega           Z reala       Z imaga  '
      write (2,*)' (Rad/s)          (ohmcm2)      (ohmcm2) '
      do 63 i=1,110
   63 write (2,57)omega(i),vreala(i)*10000,vimaga(i)*10000
      write (2,*)''
      write (2,*)'  Omega           Z realc       Z imagc  '
      write (2,*)' (Rad/s)          (ohmcm2)      (ohmcm2) '
      do 64 i=1,110
   64 write (2,57)omega(i),vrealc(i)*10000,vimagc(i)*10000
      write (2,*)''
      write (2,*)'  Omega             Z mag         phi  '
      write (2,*)' (Rad/s)           (ohmcm2)      (deg) '
      do 65 i=1,110
   65 write (2,57)omega(i),zmag(i)*10000,phi(i)
      write (2,*)''
      write (2,*)'  Omega            Z maga        phia  '
      write (2,*)' (Rad/s)           (ohmcm2)      (deg) '
      do 66 i=1,110
   66 write (2,57)omega(i),zmaga(i)*10000,phia(i)
      write (2,*)''
      write (2,*)'  Omega             Z magc       phic  '
      write (2,*)' (Rad/s)           (ohmcm2)      (deg) '
      do 67 i=1,110
   67 write (2,57)omega(i),zmagc(i)*10000,phic(i)
      write (2,*)''
      go to 999
      endif ! end of impedance section (general)

c     Obtain better initial guess by running at constant current
c     which is close to the current for the chosen constraint.
      mcL=1
        call comp(n,lim,k,rr,0,nflag,0,jcount)
        if(nconv.eq.1) write (2,*) 'comp 2 not converged' ! refine
c     write (2,*) 'cellpot 1'
        call cellpot(k,vv,0,0,lflag)

      mcL=mc(L)
      cuL=cu(L)
      curold=0.d0
      vvold=OCP

      if(mc(L).eq.-2) then ! constant power
      LS=2 ! Bypass this for speed.
      if(LS.ne.0) go to 56
      mcL=mc(L)
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 4 not converged' ! power
      curold=xx(ki2,n1+1)
      vvold=xx(kp1,nj)-xx(kp1,1)
      cur=curold-1.d0
      mcL=1
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 5 not converged' ! get Rint
      vv=xx(kp1,nj)-xx(kp1,1)
      Rint=-(vv-vvold)/(cur-curold)+RG ! internal resistance
      OCP=vv+(Rint-RG)*cur
      pwrmax=OCP**2/Rint/4.d0
   56 continue

      mcL=mc(L)
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 6 not converged' ! repeat?
      call cellpot(k,vv,1,0,lflag)
c     write (2,*) 'cellpot 3'
      endif ! treatment of constant power
c
c     rr is the size of a time step.
      if(cap1.eq.0.d0 .and. cap3.eq.0.d0) then
      rr=0.02d0
c     tmmax !for constant time steps
      else
      rr=1.5d-13 !  initial time step is 1.5d-13 second
      endif

      ed=0.d0
      heat=0.d0
      qlosstot=0.d0
      iflag=0
      L=0
   53 L=L+1 ! new leg
      mcL=mc(L)
      cuL=cu(L)
      if(L.gt.1) rr=0.d0
  123 k=k+1 ! new time step
      nt=k-1

c     adjust time step to match time of change in current
      time=ts(k-1)+rr
      if(time .ge. tot(L) .and. mc(L).lt.2) then
      rr=tot(L)-ts(k-1)
      time=tot(L)
      iflag=1
      endif
c
  129 ts(k)=ts(k-1)+rr ! new time step
      call calca(k) ! should need only one call at each time step
      dtnow=rr
c
      if(mc(L).eq.-3 .or. mc(L).ge.0) then
c     mc is 1 or 2 so run galvanostatically
c     or mc is 0; run at constant potential
c     or mc is -3; run at constant load
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 7 not converged' ! cur,pot,load
        if (nconv.eq.1) go to 533
c
      else ! iterate on power, mc is -2
c     The strategy is to use the constant-power mode in comp if
c     the possible maximum power is greater than that called for.
c     Otherwise, one wants to report the maximum power.  (This
c     is a little like using a taper at constant potential when
c     a constant-current run reaches a cutoff potential.)
c     The maximum-power mode is abandoned for a while.  2-16-08
      do 610 j=1,nj
      term_s1s(j)=term_s1(j)
  610 terms(j)=term(j)

      LS=3 ! activate this section to get estimate of pwrmax
c     at each time step.
c     This does not get a true pwrmax, since it does not iterate
c     and take the current to the true curmax.  It gives an
c     approximation by calculating Rint from a small change in the
c     current.
      if(LS.ne.0) go to 58
      mcL=1 ! comp is run in constant-current mode
      call comp(n,lim,k,rr,0,nflag,0,jcount) ! at previous current
      if(nconv.eq.1) write (2,*) 'comp 8 not converged' ! power
         if (nconv.eq.1) go to 533
      curold=cur
      vvold=xx(kp1,nj)-xx(kp1,1)
      pwr=cur*vvold
      cur=cur-1.d0 ! decrement current to get internal resistance
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 9 not converged' ! to get Rint
         if (nconv.eq.1) go to 533
      vv=xx(kp1,nj)-xx(kp1,1)
      Rint=-(vv-vvold)/(cur-curold)+RG ! internal resistance
      OCP=vv+(Rint-RG)*cur
      pwrmax=OCP**2/Rint/4.d0
      curmax=OCP/Rint/2.d0 ! estimate of current at max power
   58 continue
      if(pwrmax.gt.cu(L)) then
      mcL=mc(L) ! specified power should be possible
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 10 not converged' ! power
      else

c     iterate to find maximum power ! really don't want to go here
      mcL=-4 ! try to let comp do it
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 11 not converged' ! max power
         if (nconv.eq.1) go to 533
      vv=xx(kp1,nj)-xx(kp1,1)
      endif ! end of treatment of constant power

      endif ! end time-step calculation for set mode

      if(rr.lt.0.9999*dtnow) iflag=0 ! time step decreased in comp
      mcL=mc(L)
      if(mc(L).eq.2) then
        call cellpot(k,vv,0,0,lflag)
      else
        call cellpot(k,vv,1,0,lflag)
      endif
      if(rr.ne.0.d0) coul=coul+rr*xx(ki2,n1+1)
      if (il1.eq.1.and.mod(k,il3).eq.0) call nucamb(k,il2)
      frt=fc/r/t

      sign=1.d0
      vcut = vcutlo(L)
      IF(VV.LT.VCUT.and.cur.eq.0.d0) then
        do i=L+1,Lcurs
          tot(i) = tot(i) + time - tot(L) !Correct time markers
        enddo
        go to 888
      endif
      if(cu(L).lt.0.d0) then
         sign=-1.d0
         vcut = vcuthi(L)
      endif
      IF(sign*VV.LT.sign*VCUT .and. cu(L).ne.0.d0) then
         if(taper(L).eq.1) then !Taper the current by holding VV at VCUT
      if(mc(L).eq.1) then
            mc(L) = 0
            mcL = 0
            cu(L) = vcut
            cuL=vcut
            k=k+1
            ts(k)=ts(k-1)
      call comp(n,lim,k,0.d0,0,nflag,0,jcount)
      call cellpot(k,vv,1,0,lflag)
      if(il1.eq.1) call nucamb(k,il2)
      endif
         else                   !End profile segment
            do i=L+1,Lcurs
               tot(i) = tot(i) + time - tot(L) !Correct time markers
            enddo
            go to 888
         endif
      endif
c     check to see if cutoff potential is exceeded if mc is 2
      if (mc(L).eq.2) then
        IF ((VV.LT.TT(L) .AND. CU(L).GT.0.0) .OR.
     &  (VV.GT.TT(L) .AND. CU(L).LT.0.0)) THEN
          if (dabs(vv-tt(l)) .gt. 1.0d-04) then
            rr=rr/2.0d0
            iflag=1
            go to 129
          else
            time2=time2+rr
c     write (2,*) 'cellpot 5'
            call cellpot(k,vv,1,0,lflag)
            frt=fc/r/t
            iflag=1
          endif
        else
          iflag=0
          time2=time2+rr
c     write (2,*) 'cellpot 6'
          call cellpot(k,vv,1,0,lflag)
          frt=fc/r/t
        end if
      end if
c
c     Increasing time steps: comment out for constant time steps
      rrmax=tmmax
      if(k.le.5) rrmax=0.5d0
      jlim=int(lim/5.0d0)
      if (jlim.lt.5) jlim=5
      if(jcount.lt.jlim.and. k.gt.2 .and. rr.lt.rrmax .and.
     1iflag.eq.0.and.lcount.le.1) then
      if(xx(ki2,n1+5).lt.1200.d0) rr=rr*1.1d0

c      print *,'next time step increased to ', rr,' (s)'
      end if

      if(k.GE.maxt-1) then
        write (2,*) 'kmax=',k,' a larger matrix needed for xt'
      endif
      if (k.GE.501) then ! trim stored solid concentrations
      do 92 kk=3,401,2
      kput=(KK+1)/2
      ts(kput)=ts(kk)
      do 92 j=1,nj
      do 92 i=1,n
   92 xt(i,j,kput)=xt(i,j,kk)
      do 93 kk=402,K
      ts(kk-200)=ts(kk)
      do 93 j=1,nj
      do 93 i=1,n
   93 xt(i,j,kk-200)=xt(i,j,kk)
      k=k-200
      endif

      if (iflag.eq.0 .and. rr.eq.0.d0) then
c     rr is the size of a time step.
      if(cap1.eq.0.d0 .and. cap3.eq.0.d0) then
      rr=0.02d0
c      tmmax !for constant time steps
      else
      rr=1.5d-13 !  initial time step is 1.5d-13 second
      endif
      endif
      if (iflag .eq. 0) go to 123 ! a new time step
  888 continue
      iflag=0
      if (mc(L).eq.2) then
      do 124 m=L,Lcurs
  124 tot(m)=tot(m)+time2
      time2=0.0d0
      end if
      IF(L.EQ.LCURS .AND. LCURS.GE.100) THEN 
c If Lcurs greater than 100 go into an endless loop
      L=0
      tot(1)=TOT(LCURS)
      if (mc(1).lt.2) TOT(1)=TOT(1)+60.0D0*TT(1)
      do 403 i=2,Lcurs
      if (mc(i).lt.2) then
      tot(i)=tot(i-1)+6.0d01*tt(i)
      else
      tot(i)=tot(i-1)
      end if
  403 continue
      ENDIF

  533 if(L.lt.Lcurs) then ! prepare for a new leg
      rr=0.d0
c       calculate zero-time solution for change in current
c     Need some initialization here to get started smoothly on
c     a new leg of the calculation.
      cuL=cu(L+1)
c     First guess for current
      cur=cu(L+1) ! constant current
      if(mc(L+1).eq.-3) cur=OCP/(cu(L+1)+Rint) ! constant load
      if(mc(L+1).eq.0) cur=(OCP-cu(L+1))/Rint ! constant potential

      if(mc(L+1).eq.-2) then ! constant power
      LS=1
      if(LS.eq.1) go to 55 ! bypass maximum-power calculations
c     first calculate power maximum
      curmax=OCP/Rint/2.d0
      pwrmax=OCP**2/Rint/4.d0

      do ii=1,4 ! could use fewer iterations here
      cur=curmax
      mcL=1
      call comp(n,lim,k,rr,0,nflag,0,jcount) 
      if(nconv.eq.1) write (2,*) 'comp 12 not converged'
      curold=cur ! save correctly calculated values at one cur
      vvold=xx(kp1,nj)-xx(kp1,1)
      cur=cur-1.d0 ! decrement current
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 13 not converged'
      vv=xx(kp1,nj)-xx(kp1,1)
      Rint=-(vv-vvold)/(cur-curold)+RG
      OCP=vv+(Rint-RG)*cur
      curmax=OCP/Rint/2.d0
      pwrmax=OCP**2/Rint/4.d0
      enddo ! ii
      if(pwrmax.lt.cuL) then ! can't deliver set power
      cur=curmax
      call comp(n,lim,k,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 14 not converged'
      mcL=mc(L)
      go to 53
      endif

   55 continue
      cur=cu(L+1)/OCP
      do ii=1,6 ! refine estimate of current for given power
c     could solve quadratic
      cur=cur-(cur*OCP-Rint*cur**2-cu(L+1))/(OCP-2.d0*Rint*cur)
      enddo
      endif ! end of power loop

      avgeta=cur*Rint/2.d0
      if(mc(L+1).eq.0) avgeta=(OCP-cuL)/2.d0

      Ua=xx(kp1,1)
      Uc=xx(kp1,nj)
      Ua=0.d0 ! need better open-circuit potentials
      Uc=4.d0 ! need better open-circuit potentials
      do 73 j=1,(n1+1) ! negative electrode
      do mpa=1,npa
      xt(2+mpa,j,k+0)=xx(2+mpa,j)
      enddo ! mpa
      xt(ki2,j,k+0)=cur*dble(j-1)/dble(n1)
      xt(kj,j,k+0)=cur/fc/h1/area1/dble(n1)
   73 xt(kp1,j,k+0)=xx(kp2,j)+Ua+0.60d0*avgeta

      do j=n1+1,n1+n2 ! separator
      xt(ki2,j,k+0)=cur
      enddo
c
      do 75 j=(n2+n1),nj ! positive electrode
      do mpa=1,npa
      xt(2+mpa,j,k+0)=xx(2+mpa,j)
      enddo ! mpa
      xt(ki2,j,k+0)=cur*dble(nj-j)/dble(n3)
      xt(kj,j,k+0)=-cur/fc/h3/area3/dble(n3)
   75 xt(kp1,j,k+0)=xx(kp2,j)+Uc-1.20d0*avgeta

c     straighten out the initialization
      mcL=1
      call comp(n,lim,k+1,rr,0,nflag,0,jcount)
      if(nconv.eq.1) write (2,*) 'comp 15 not converged'
      mcL=mc(L+1)

        if(mc(L+1) .gt. 0) then
      if (nconv.eq.0) k=k+1
          ts(k)=ts(k-1)
          rr=0.0d0 !run a zero time step at the start
c     mcL=1
          call comp(n,lim,k,rr,0,nflag,0,jcount) 
      if(nconv.eq.1) write (2,*) 'comp 16 not converged'
c     write (2,*) 'cellpot 7'
          call cellpot(k,vv,1,0,lflag)
      if (il1.eq.1.and.mod(k,il3).eq.0) call nucamb(k,il2)
        endif
      endif
c      rr = tmmax !uncomment and comment out next 5 lines for constant time steps
      if(cap1.eq.0.d0 .and. cap3.eq.0.d0) then
        rr=0.02d0   !  initial time step is 0.2 seconds
      else
        rr=1.5d-13 !  initial time step is 1.5d-13 second
      endif
      nconv=0
      if(L.le.Lcurs) then
      if (L.eq.1) pow=ed/tw/ts(k)
      if (L.ne.1) pow=ed/tw/(time-timesave)
      timesave=ts(k)
c     write (2,*) ' '
      if (L.eq.Lcurs) write (2,44) tw
      write (2,48) ed/tw/3.6d03
      write (2,49) pow
      write (2,47) heat
      write (2,*) ' '
      ed=0.0d0
      endif
      if(jsol.gt.0 .and. jsol.lt.nj) call sol(k,jsol)
      IF(L.LT.LCURS) GO TO 53 ! new leg

c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c     Additional features section:
c
c     peak-power subroutine:
      if(lpow.eq.1) then
      il1=0
      call peak(n,lim,cu(1), vcut)
      endif
c
c     Solid-phase concentration profiles at given time and position
      if(jsol.gt.0 .and. jsol.lt.nj) call sol(k,jsol)

      stop
  999 end
c
c***********************************************************************
      subroutine comp(n,lim,kk,tau,kkflag,nflag,lpow,jcount)
c     Main calculational subroutine, going through the
c     equations that are solved through the cell sandwich.
      implicit real*8(a-h,o-z)
      parameter(maxt=900)
      common /n/ tmmax,imp,ji,nx,nt,n1,n2,nj,n3,nconv,npa,iSk3,kj3
     &,ii2,ki2,kj,i2div,ip2,kp2,ip1,kp1,imb2,kS2,imb1,kS1,iSk1,kj1
     &,iSk2,kj2
      common/activ/ EbarD,Ebarkap,Ebarka,Ebarkc,Ebarr1,Ebarr3,
     &Ebars1,Ebars3,Ebarks1a,Ebarks1c,Ebarks2a,Ebarks2c,
     &Ebarks3a,Ebarks3c
      common /calc/ ai(maxt,5),ai2(maxt,5),ts(maxt),h,h1,h2,h3,hcn,
     1hcp,rr,rrmax,cuL,sumpa(5,221),Rad1pa(5),Rad3pa(5),area1pa(5)
     &,area3pa(5),mcL
      common/pindiv/pwfindiv(5,221),pwfindiv_d(5,221),vf1(5),vf3(5),izfl
      common/const/ fc,r,t,frt,cur,ep3,ep2,pi,ep1,epf3,epf1,
     &epp1,epp2,epp3,shape3,shape1,capp1,capp3,nneg,nprop,npos
      common/gas/ epg1,epg2,epg3
      common/power/ ed,Vold,ranodesave,rcathdesave,heat,qlosstot
      common/ssblock/ xp0(16),xx0(16,221),term(221)
      common/var/ xp(16),xx(17,221),xt(16,221,maxt)
      common/imp/ vreal(110),vimag(110),phi(110),zmag(110),omega(510),
     1vreala(110),vrealc(110),vimaga(110),vimagc(110),
     1phia(110),phic(110),zmaga(110),zmagc(110)
      common/cprop/ sig3,area3,rka3,rka3save,ct3,dfs3,Rad3,cap1,cap3,
     1sig1,area1,rka1,rka1save,ct1,dfs1,Rad1,tw,dfs1save,dfs3save
      common/tprop/df(221),cd(221),tm(221),ddo2(221),ddh2(221),
     1ddf(221),dcd(221),dtm(221),dfu(221),d2fu(221),do2(221),dh2(221)
      common/temp/ thk,htc,dudt,Cp,dens,tam,g0,qq,qloss,residm,ncell,lht
      common/side/rksc1,c1init,c2init,rksa1,term_s1(221),vol,
     &rksa2,rksa3,rksc3,rksc2,UsO2,UsH2,cn2,term_s2(221),nside
      common/mat/ b,d
      common/bnd/ a,c,g,x,y
      COMMON /vdc/ aa(1,1),bb(1,1),cc(1,150),dd(1,3),gg(1),xxx(1,1),
     1yy(1,1),hha,hhc,cssold(220,150),css(220,150),utz(5,221),
     &dsold(220,150),ds(220,150),time,nn,nnj,np,mvdc1,mvdc3,lims
      common /maxpow/ lcount
      common/RG/ RG,RGn,RGp,RGext
      save iijj
      dimension b(17,17),d(17,35),termn(221)
      dimension termn_s1(221), termn_s2(221)
      dimension a(17,17),c(17,221),g(17),x(17,17),y(17,17)
   99 format (1h ,//5x,'this run just did not converge'//)
  103 format(f13.5,'      ',f12.6,'     ',f12.6)
  104 format (i4,' ',g15.5,' ',g15.5,' ',
     &g15.5,' ',g15.5,' ',g15.5,' ',g15.5)
      ntot=n
      if(imp.eq.1) ntot=2*n
      nx=ntot
      mcLL=mcL
      if(mcLL.eq.-4) mcL=1 ! use const cur for maximum power
      Lcount=0 ! counter for maximum power
      curold=cur ! save correctly calculated values at one cur
      vvold=xx(kp1,nj)-xx(kp1,1)
      rrold=rr
      rr=tau
      if(mcLL.eq.-4) cur=cur-1.d0

      omi=dble(1-imp)
      exbrug=1.5d0
      if(nprop.eq.7 .or. nprop.eq.8) exbrug=3.3d0
c
  666 kadd=0
c     sets first guess to last time-step values
      if(rr.eq.0.d0 .and. lpow.eq.1) kadd=1
      if(kk .eq. 1) kadd = 1
      do 1 j=1,nj
      do 1 i=1,ntot
      if(rr.ne.0.d0) xx(i,j)=xt(i,j,kk-1+kadd)
    1 c(i,j)=xx(i,j)

      do 2 j=1,nj
      do 2 i=1,nnj
    2 css(j,i)=cssold(j,i)
c
c     initialize variables to begin each iteration (jcount is iteration #)
    3 Lcount=Lcount+1
      jcount=0
      do 4 i=1,n
    4 xp(i)=0.0d0
c
    8 j=0
      jcount=jcount+1
      jcountsum=jcountsum+1
      if(mcL.le.-2 .or. mcL.eq.0) cur=xx(ki2,n1+1)
c     calculate physical properties
      call prop(nj,n2,n1)
c
c     initialize x and y for iteration
      do 9 i=1,ntot
      do 9 k=1,ntot
      x(i,k)=0.0d0
    9 y(i,k)=0.0d0
c
c     store previous iteration of (xp in xp0)  &  (xx in xx0)
      do 6 i=1,n
      xp0(i)=xp(i)
      do 6 jj=1,nj
    6 xx0(i,jj) = xx(i,jj)  

c     for a given iteration, set up governing equations and bc's
c     start at the left interface and move across cell.
c     initialize a,b,d,g arrays for each node
c
   10 j=j+1
      hC=h2
      if(j.le.n1+1)  hC=h1
      if(j.le.n1+1)  area=area1
      if(j.le.n1+1)  capp=capp1
      if(j.le.n1+1)  cap=cap1
      if(j.ge.n1+n2) hC=h3
      if(j.ge.n1+n2) area=area3
      if(j.ge.n1+n2) capp=capp3
      if(j.ge.n1+n2) cap=cap3
      if(j.eq.1 .or. j.eq.n1+1 .or. j.eq.n1+n2 .or. j.eq.nj)
     &area=area/2.d0

c
      do 11 i=1,ntot+1
      g(i)=0.0d0
      xx(i,j)=c(i,j)
      do 11 k=1,ntot+1
      a(i,k)=0.0d0
      b(i,k)=0.0d0
   11 d(i,k)=0.0d0

c
c     Added if-then construct to account for impedance.  While
c     impedance will be taken at zero time, we must solve for
c     concentration profiles. (JPM - 29 Jun 95)
      if(imp.eq.1) then
      rr=1.0d0 ! set rr to 1.0 to avoid divide by 0
      else ! end of impedance section (preliminary)
      if(rr.le.0.0d0) then ! treat as a zero time step
      b(1,1)=1.0d0
      g(1)=xt(1,j,kk-1+kadd)-c(1,j) ! fix electrolyte concentration
      if (nside.ge.1) then
      b(imb1,kS1)=1.0d0
      g(imb1)=xt(kS1,j,kk-1+kadd)-c(kS1,j) !fix oxygen concentration
      if (nside.ge.2) then
      b(imb2,kS2)=1.0d0
      g(imb2)=xt(kS2,j,kk-1+kadd)-c(kS2,j) !fix hydrogen concentration
      endif
      endif
      go to 199
      endif
      endif
c     ________________________________________________
c     Equation 1, mass balance.
c
c     Adsorption behavior is fixed w/ dqn parameter
c     dqn is the dq-/dq derivative equal to either -1 (anion
c     adsorption) or 0 (cation adsorption)
c 111 dq1=0.0d0
c     dq3=0.0d0
  111 dq1=-1.0d0
      dq3=-1.0d0 ! Johnson and Newman, J. Electrochem. Soc., 118 (1971), 510-517.

      termn(j)=0.d0

c Set the stoichiometric coefficient according to the nature of the reacting ion.
c For a cation, set so=1.0d0, for an anion, set so=-1.0d0.

      so=1.0d0
      if (nprop.eq.13) so=-1.d0 ! for an anion in the intercalation reaction

      fac=1.d0
      if(j.eq.n1+2 .and. n1 .gt. 0)
     &fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2+1) fac=((ep3+epp3)/(ep2+epp2))**exbrug
      epn=ep2+epp2
      hn=h2
      acn=0.0d0
      dqn=0.0d0
      if(n1.gt.0 .and. j.le.n1+1) then
        epn=ep1+epp1
        hn=h1
        acn=area1*capp1
        dqn=dq1
      elseif(j.gt.n1+n2) then
        epn=ep3+epp3
        hn=h3
        acn=area3*capp3
        dqn=dq3
      endif

      if(j.gt.1) then ! deal with box to left of point.
        termn(j)=-(df(j)+fac*df(j-1))*(c(1,j)-c(1,j-1))/hn/2.d0
     &  -so*(1.d0-0.5d0*(tm(j)+tm(j-1)))*c(ki2,j-1)/fc
        a(1,1)=epn*hn*0.125d0/rr*omi
     &  -(df(j)+fac*df(j-1))/hn/4.d0+fac*ddf(j-1)*
     &  (c(1,j)-c(1,j-1))/hn/4.d0-so*dtm(j-1)*c(ki2,j-1)/4.d0/fc
        b(1,1)=epn*hn*0.375d0/rr*omi
     &  +(df(j)+fac*df(j-1))/hn/4.d0+ddf(j)*(c(1,j)-c(1,j-1))/hn/4.d0
     &  -so*dtm(j)*c(ki2,j-1)/4.d0/fc
        a(1,ki2)=so*(1.d0-0.5d0*(tm(j)+tm(j-1)))/2.d0/fc
        g(1)=-epn*hn*(0.375d0*(c(1,j)-xt(1,j,kk-1+kadd))
     &  +0.125d0*(c(1,j-1)-xt(1,j-1,kk-1+kadd)))/rr*omi
        if (imp.eq.1) then
          b(1,1+n)=-epn*hn*0.375d0*omega(ji)  /2.0d0
          a(1,1+n)=-epn*hn*0.125d0*omega(ji)  /2.0d0
          b(1+n,1)=epn*hn*0.375d0*omega(ji)   /2.0d0
          a(1+n,1)=epn*hn*0.125d0*omega(ji)   /2.0d0
        endif ! end of impedance section (Equation 1)
      endif

      fac=1.d0
      if(n1.gt.0 .and. j.eq.n1+1)
     &fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2) fac=((ep3+epp3)/(ep2+epp2))**exbrug
      epn=ep2+epp2
      hn=h2
      acn=0.0d0
      dqn=0.0d0
      if(n1.gt.0 .and. j.lt.n1+1) then
        epn=ep1+epp1
        hn=h1
        acn=area1*capp1
        dqn=dq1
      elseif(j.ge.n1+n2) then
        epn=ep3+epp3
        hn=h3
        acn=area3*capp3
        dqn=dq3
      endif
      if(j.ne.nj) then ! deal with box to right of point.
c     At the foil anode, only the box to the right should be used
        termn(j)=termn(j)-(fac*df(j)+df(j+1))*(c(1,j)-c(1,j+1))/hn/2.d0
     &  +so*(1.d0-0.5d0*(tm(j)+tm(j+1)))*c(ki2,j)/fc
        d(1,1)=epn*hn*0.125d0/rr*omi
     &  -(fac*df(j)+df(j+1))/hn/4.d0+ddf(j+1)*(c(1,j)-c(1,j+1))/hn/4.d0
     &  +so*dtm(j+1)*c(ki2,j)/4.d0/fc
        b(1,1)=b(1,1)+epn*hn*0.375d0/rr*omi
     &+(fac*df(j)+df(j+1))/hn/4.d0+fac*ddf(j)*(c(1,j)-c(1,j+1))/hn/4.d0
     &  +so*dtm(j)*c(ki2,j)/4.d0/fc
        b(1,ki2)=b(1,ki2)-so*(1.d0-0.5d0*(tm(j)+tm(j+1)))/2.d0/fc
        g(1)=g(1)-epn*hn*(0.375d0*(c(1,j)-xt(1,j,kk-1+kadd))
     &  +0.125d0*(c(1,j+1)-xt(1,j+1,kk-1+kadd)))/rr*omi
        if (imp.eq.1) then
          b(1,1+n)=b(1,1+n)-epn*hn*0.375d0*omega(ji) /2.0d0
          d(1,1+n)=-epn*hn*0.125d0*omega(ji)       /2.0d0
          b(1+n,1)=b(1+n,1)+epn*hn*0.375d0*omega(ji) /2.0d0
          d(1+n,1)=epn*hn*0.125d0*omega(ji)        /2.0d0
        endif ! end of impedance section (equation 1)
      endif
      g(1)=g(1)+(termn(j)+term(j))/2.d0

      if (imp.eq.1) g(1)=0.0d0

c     ________________________________________________
c Equation ii2, material balance in solid insertion material.
c This material balance includes the option of having side-reaction
c species 2 (hydrogen in the NiMH case) absorb or desorb from
c the solid phase in addition to the electrochemical intercalation
c or deintercalation of hydrogen. PSA 4/22/06

      if(imp.eq.0) then ! impedance is off

        if(j.le.n1+1) then ! in the anode
          if (n1 .eq. 0) then  ! foil anode
            mpa=1 ! only one size with foil anode
            g(ii2) = -c(2+mpa,j)
            b(ii2,2+mpa) = 1.0d0
          else ! nonfoil anode
            if (mvdc1.eq.0) then ! superposition
              if (nside.eq.3) then ! side reaction with hydrogen
                mpa=1 ! only one size with side reactions
c     If side reaction 3 does not involve absorption or desorption,
c     remove c(kj3,j) and modify the appropriate Jacobian.  This feature
c     is included to treat hydrogen absorption/desorption in the NiMH system.
                b(ii2,kj3)=1.0d0/Rad1pa(mpa)
                b(ii2,2+mpa)=ai2(1,mpa)/rr*omi
                b(ii2,kj)=1.0d0/Rad1pa(mpa)
                g(ii2)=ai2(1,mpa)*xt(2+mpa,j,kk-1)/rr*omi-sumpa(mpa,j)-
     &       ai2(1,mpa)/rr*omi*c(2+mpa,j)-(c(kj,j)+c(kj3,j))/Rad1pa(mpa)
              else !no side reaction 3
                b(ii2,kj)=1.d0
                g(ii2)=-c(kj,j)
                do mpa=1,npa ! Calculate g, no side reaction
                  b(ii2,2+mpa)=ai2(1,mpa)/rr*omi*Rad1pa(mpa)
     &            *area1pa(mpa)/area1
                  g(ii2)=g(ii2)+area1pa(mpa)/area1*(ai2(1,mpa)*(xt(2+mpa
     &            ,j,kk-1)-c(2+mpa,j))/rr*omi-sumpa(mpa,j))*Rad1pa(mpa)
                enddo ! mpa

              endif ! end side (and no side) reactions

c Special section for VDC -anode
            elseif (mvdc1.eq.1) then
              mpa=1 ! only one size for variable solid-phase diffusion coefficient
              call vardc(j,jcount,kk)
              b(ii2,2+mpa)=1.0D0
              g(ii2)=css(j,nnj)-xx(2+mpa,j)
            endif ! end calculate b's
          endif ! end nonfoil anode

        elseif (j.ge.n1+n2) then !cathode
          if (mvdc3.eq.0) then ! superposition
            if (nside.eq.3) then ! side reaction with hydrogen
              mpa=1 ! only one size with side reactions
c     If side reaction 3 does not involve absorption or desorption,
c     remove c(kj3,j) and modify the appropriate Jacobian.  This feature
c     is included to treat hydrogen absorption/desorption in the NiMH system.
              b(ii2,2+mpa)=ai(1,mpa)/rr*omi
              b(ii2,kj)=1.0d0/Rad3pa(mpa)
              b(ii2,kj3)=1.0d0/Rad3pa(mpa)
              g(ii2)=ai(1,mpa)*xt(2+mpa,j,kk-1)/rr*omi-sumpa(mpa,j)-
     &      ai(1,mpa)/rr*omi*c(2+mpa,j)-(c(kj,j)+c(kj3,j))/Rad3pa(mpa)
            else ! no side reaction 3
              b(ii2,kj)=1.d0
              g(ii2)=-c(kj,j)
              do mpa=1,npa
                b(ii2,2+mpa)=ai(1,mpa)/rr*omi*Rad3pa(mpa)
     &          *area3pa(mpa)/area3
                g(ii2)=g(ii2)+area3pa(mpa)/area3*(ai(1,mpa)*(xt(2+mpa,j
     &            ,kk-1)-c(2+mpa,j))/rr*omi-sumpa(mpa,j))*Rad3pa(mpa)
              enddo !mpa

            endif !end side reactions

c Special section of VDC - cathode
         elseif (mvdc3.eq.1) then !calculate b's
           mpa=1 ! only one size for variable solid-phase diffusion coefficient
           call vardc(j,jcount,kk)
           b(ii2,2+mpa)=1.0D0
           g(ii2)=css(j,nnj)-xx(2+mpa,j)
         endif ! end calculate b's

        else  !separator
          b(ii2,4+npa)=1.0d0
          g(ii2)=-c(4+npa,j) ! set solid flux = 0 in separator
        endif

      else ! Impedance section

      mpa=1 ! only one size for impedance
        if(j.eq.1 .and. n1.eq.0) then
c No solid concentration for foil electrode:
        b(ii2,2+mpa)=1.0d0
        g(ii2)=-c(2+mpa,j)
        b(ii2+n,2+mpa+n)=1.0d0
        g(ii2+n)=-c(2+mpa+n,j)
        g(ii2)=0.0d0
         elseif(j.gt.n1+1 .and. j.lt.n1+n2) then ! separator
        b(ii2,2+mpa)=1.0d0
        b(ii2+n,2+mpa+n)=1.0d0
        g(ii2)=0.d0
         else ! anode and cathode
        if(j.eq.1) then
        Rad=Rad1
        dfs=dfs1
        elseif(j.eq.n1+n2) then
        Rad=Rad3
        dfs=dfs3
        endif
c Solution of diffusion equation w/Laplace gives:
        b(ii2,kj)=-1.0d0
        b(ii2+n,kj3)=-1.0d0
c New section added April 4, 1998, JSN
        arg=Rad*(2.d0*omega(ji)/dfs)**0.5d0
        e1=expg(-arg)
        e2=e1**2
        if(arg.lt.0.22d0) then
        denom=2.d0*e1*arg**2*(1.d0+arg**4/3.6d2+arg**8/1.8144d6)
        diff2=2.d0*e1*(2.d0*arg+arg**5/6.d1+arg**9/1.99584d7*110.d0)
        diff3=2.d0*e1*(arg**3/3.d0+arg**7/2.52d3+arg**11/1.99584d7)
        Ysr=dfs/Rad*(arg**4/1.8d2+arg**8/4.536d5+arg**12/7.2648576d9)
     %  /1.d0/(1.d0+arg**4/3.6d2+arg**8/1.8144d6)
        else
        denom=1.d0+e2-2.d0*e1*dcos(arg)
        diff2=1.d0-e2+2.d0*e1*dsin(arg)
        diff3=1.d0-e2-2.d0*e1*dsin(arg)
        Ysr=-dfs/Rad*(1.d0-arg/2.d0*diff2/denom)
        endif
        Ysi=dfs/Rad*arg/2.d0*diff3/denom
        real1=-Ysr
        comp1=Ysi
        b(ii2,2+mpa)=real1
        b(ii2+n,2+mpa+n)=real1
        b(ii2,2+mpa+n)=comp1
        b(ii2+n,2+mpa)=-comp1
        g(ii2)=0.0d0
         endif
      g(ii2+n)=0.0d0
      endif ! impedance section (equation ii2)

c     ________________________________________________

c Equation imb2 for variable kS2
c Material balance on the second side-reaction species.
c This is set up for a balance on hydrogen for the NiMH system.

      if (nside.ge.2) then

      sf=1.0d0/2.0d0  !stoichiometric factor for H2 reaction

      termn_s2(j)=0.d0

c aream the area for use in mass balance

      if (j.le.n1+1) aream=area1
      if (j.ge.n1+n2) aream=area3

      fac=1.d0
      if(j.eq.n1+2 .and. n1 .gt. 0)
     &fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2+1) fac=((ep3+epp3)/(ep2+epp2))**exbrug
      epn=ep2+epp2

      hn=h2  !box to left of point.
      epn=epg2
      if (j.le.n1+1) then
      hn=h1
      epn=epg1
      elseif (j.gt.n1+n2) then
      hn=h3
      epn=epg3
      endif

      if (j.eq.n1+n2) sf=0.0d0 !half box w/o generation
      if (j.eq.n1+2) sf=0.0d0 !half box w/0 generation

      if (j.gt.1) then !box to left of point
      termn_s2(j)=-(dh2(j)+fac*dh2(j-1))*(c(kS2,j)-c(kS2,j-1))/hn/2.0d0
     &+sf*hn*aream*vol*(0.375d0*(c(kj3,j)+c(kj2,j))+0.125d0*
     &(c(kj3,j-1)+c(kj2,j-1)))
      g(imb2)=-epn*hn*(0.375d0*(c(kS2,j)-xt(kS2,j,kk-1+kadd))
     &+0.125d0*(c(kS2,j-1)-xt(kS2,j-1,kk-1+kadd)))/rr
      b(imb2,kS2)=epn*hn*0.375d0/rr+(dh2(j)+
     &fac*dh2(j-1))/hn/4.0d0+ddh2(j)*(c(kS2,j)-c(kS2,j-1))/hn/4.0d0
      a(imb2,kS2)=-(dh2(j)+fac*dh2(j-1))/hn/4.0d0+hn*epn*0.125d0/rr
     &+fac*ddh2(j-1)*(c(kS2,j)-c(kS2,j-1))/hn/4.0d0
      b(imb2,kj3)=-0.375d0*hn*sf*aream*vol/2.0d0
      b(imb2,kj2)=-0.375d0*hn*sf*aream*vol/2.0d0
      a(imb2,kj3)=-0.125d0*hn*sf*aream*vol/2.0d0
      a(imb2,kj2)=-0.125d0*hn*sf*aream*vol/2.0d0
      endif

      sf=1.0d0/2.0d0

      fac=1.d0
      if(j.eq.n1+1 .and. n1 .gt. 0)
     &fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2) fac=((ep3+epp3)/(ep2+epp2))**exbrug
      epn=ep2+epp2

      hn=h2 !deal with box to right of point.
      epn=epg2
      if (j.lt.n1+1) then
      hn=h1
      epn=epg1
      elseif (j.ge.n1+n2) then
      hn=h3
      epn=epg3
      endif

      if (j.eq.n1+1) sf=0.0d0 !half box w/o generation
      if (j.eq.n1+n2-1) sf=0.0d0 !half box w/o generation

      if (j.lt.nj) then !box to right of point
      termn_s2(j)=termn_s2(j)-(fac*dh2(j)+dh2(j+1))*
     &(c(kS2,j)-c(kS2,j+1))/hn/2.0d0+
     &sf*hn*aream*vol*(0.375d0*(c(kj3,j)+c(kj2,j))+0.125d0*
     &(c(kj3,j+1)+c(kj2,j+1)))
      g(imb2)=g(imb2)-epn*hn*(0.375d0*(c(kS2,j)-xt(kS2,j,kk-1+kadd))
     &+0.125d0*(c(kS2,j+1)-xt(kS2,j+1,kk-1+kadd)))/rr

      b(imb2,kS2)=b(imb2,kS2)+epn*hn*0.375d0/rr+(fac*dh2(j)+dh2(j+1))
     &/hn/4.0d0+fac*ddh2(j)*(c(kS2,j)-c(kS2,j+1))/hn/4.0d0
      b(imb2,kj3)=b(imb2,kj3)-sf*hn*aream*vol*0.375d0/2.0d0
      b(imb2,kj2)=b(imb2,kj2)-sf*hn*aream*vol*0.375d0/2.0d0
      d(imb2,kS2)=-(fac*dh2(j)+dh2(j+1))/hn/4.0d0+ddh2(j+1)*
     &(c(kS2,j)-c(kS2,j+1))/hn/4.0d0+epn*hn*0.125d0/rr
      d(imb2,kj3)=-sf*hn*aream*vol*0.125d0/2.0d0
      d(imb2,kj2)=-sf*hn*aream*vol*0.125d0/2.0d0
      endif

      g(imb2)=g(imb2)+(termn_s2(j)+term_s2(j))/2.0d0

      endif

c     ________________________________________________

c Equation imb1 for variable kS1
c Material balance on the first side-reaction species.
c This is set up for a balance on oxygen for the NiMH system.

      if (nside.ge.1) then
      sf=-1.0d0/4.0d0  !stoichiometric factor for O2 reaction

      termn_s1(j)=0.d0

      if (j.le.n1+1) aream=area1
      if (j.ge.n1+n2) aream=area3

      fac=1.d0
      if(j.eq.n1+2 .and. n1 .gt. 0)
     &fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2+1) fac=((ep3+epp3)/(ep2+epp2))**exbrug
      epn=ep2+epp2

      hn=h2  !box to left of point.
c For liquid-phase side reaction see User's Manual.
c     epn=ep2
      epn=epg2
      if (j.le.n1+1) then
      hn=h1
c     epn=ep1
      epn=epg1
      elseif (j.gt.n1+n2) then
      hn=h3
c     epn=ep3
      epn=epg3
      endif

      if (j.eq.n1+n2) sf=0.0d0 !half box w/o generation
      if (j.eq.n1+2) sf=0.0d0 !half box w/0 generation
      if (j.gt.1) then !box to left of point
      termn_s1(j)=-(do2(j)+fac*do2(j-1))*(c(kS1,j)-c(kS1,j-1))/hn/2.0d0
     &+sf*hn*aream*vol*(0.375d0*c(kj1,j)+0.125d0*c(kj1,j-1))
      g(imb1)=-epn*hn*(0.375d0*(c(kS1,j)-xt(kS1,j,kk-1+kadd))
     &+0.125d0*(c(kS1,j-1)-xt(kS1,j-1,kk-1+kadd)))/rr
      b(imb1,kS1)=epn*hn*0.375d0/rr+(do2(j)+fac*do2(j-1))/hn/4.0d0
     &+ddo2(j)*(c(kS1,j)-c(kS1,j-1))/hn/4.0d0
      a(imb1,kS1)=-(do2(j)+fac*do2(j-1))/hn/4.0d0+epn*hn*0.125d0/rr
     &+fac*ddo2(j-1)*(c(kS1,j)-c(kS1,j-1))/hn/4.0d0
      b(imb1,kj1)=-0.375d0*hn*sf*aream*vol/2.0d0
      a(imb1,kj1)=-0.125d0*hn*sf*aream*vol/2.0d0
      endif

      sf=-1.0d0/4.0d0

      fac=1.d0
      if(j.eq.n1+1 .and. n1 .gt. 0)
     &fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2) fac=((ep3+epp3)/(ep2+epp2))**exbrug
      epn=ep2+epp2

      hn=h2 !deal with box to right of point.
c For liquid-phase side reaction see User's Manual.
c     epn=ep2
      epn=epg2
      if (j.lt.n1+1) then
      hn=h1
c     epn=ep1
      epn=epg1
      elseif (j.ge.n1+n2) then
      hn=h3
c     epn=ep3
      epn=epg3
      endif

      if (j.eq.n1+1) sf=0.0d0 !half box w/o generation
      if (j.eq.n1+n2-1) sf=0.0d0 !half box w/o generation

      if (j.lt.nj) then !box to right of point
      termn_s1(j)=termn_s1(j)-(fac*do2(j)+do2(j+1))*(c(kS1,j)
     &-c(kS1,j+1))
     &/hn/2.0d0+sf*hn*aream*vol*(0.375d0*c(kj1,j)+0.125d0*c(kj1,j+1))
      g(imb1)=g(imb1)-epn*hn*(0.375d0*(c(kS1,j)-xt(kS1,j,kk-1+kadd))
     &+0.125d0*(c(kS1,j+1)-xt(kS1,j+1,kk-1+kadd)))/rr

      b(imb1,kS1)=b(imb1,kS1)+epn*hn*0.375d0/rr+(fac*do2(j)+do2(j+1))
     &/hn/4.0d0+fac*ddo2(j)*(c(kS1,j)-c(kS1,j+1))/hn/4.0d0
      b(imb1,kj1)=b(imb1,kj1)-sf*hn*aream*vol*0.375d0/2.0d0
      d(imb1,kS1)=-(fac*do2(j)+do2(j+1))/hn/4.0d0+epn*hn*0.125d0/rr-
     &ddo2(j+1)*(c(kS1,j)-c(kS1,j+1))/hn/4.0d0
      d(imb1,kj1)=-sf*hn*aream*vol*0.125d0/2.0d0
      endif

      g(imb1)=g(imb1)+(termn_s1(j)+term_s1(j))/2.0d0

      endif

  199 if(imp.eq.1) g(12)=0.0d0

c     ________________________________________________
c     Equation ip2, Ohm's law in the liquid phase
c     There are separate versions of this equation for a Li-ion
c     system and for the NiMH system.
      if(j.le.n1) then
      h=h1
      else if(j.lt.n1+n2) then
      h=h2
      else
      h=h3
      endif
      fac=1.0d0
      if(j.eq.n1+1 .and. n1 .gt. 0)
     1fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2)
     1fac=((ep3+epp3)/(ep2+epp2))**exbrug
c
         if(j.eq.nj) then
      b(ip2,kp2)=1.0d0
      g(ip2)=-c(kp2,j) !boundary condition on potential
      if (imp.eq.1) g(ip2)=0.0d0
      go to 12
         endif

      if (nprop.eq.13) then 
c Ohm's law in solution for NiMH different from Lithium
      dcf=(c(1,j+1)-c(1,j))/h 
      r1=(c(1,j+1)+c(1,j))/2.0d0
      r4=c(ki2,j)
      p1=1.0d0-(tm(j)+tm(j+1))/2.0d0 !different from used in mass balance
      p2=(fac*cd(j)+cd(j+1))/2.0d0
      p3=(fac*dcd(j)+dcd(j+1))/2.0d0
      H2OMW=18.016D0
      PKOHMW=56.11D0
      c1=r1/1.0d6
      fu= 1.004D0-36.23D0*c1**0.5+1374.3D0*c1-17850.7d0*c1**1.5d0+
     &55406.0D0*c1**2d0+7.16856D05*c1**2.5d0
      d_fu=-18.115D0*c1**-0.5+1374.3D0-2.6776D04*c1**0.5
     &+1.10812D05*c1 + 1.7921D06*c1**1.5
      d2_fu=9.0575D0/c1**1.5-13388.0D0/c1**0.5+1.10812D05+
     & 2.6882D06*c1**0.5
      d_fu=d_fu/1.0d6
      d2_fu=d2_fu/1.0d12
      densw=1.001D0+47.57D0*c1-776.22D0*c1**2
      dens_1d = 47.57D0 - 1552.44D0*c1
      dens_2d= -1552.44D0
      densw=densw*1.0d6
      dens_1d=dens_1d*1.0d0
      dens_2d=dens_2d/1.0d6

      c_solv=(densw-c(1,j)*PKOHMW)/H2OMW

      U1=densw-c(1,j)*PKOHMW
      U2=dens_1d-PKOHMW
      U3=-dcf*((p1*U1+H2OMW*c(1,j))*((fu*d2_fu-d_fu**2.0D0)/
     &(fu**2.0D0)-1.0D0/c(1,j)**2.0D0)+(d_fu/fu+1.0D0/xx(1,j))
     &*(p1*U2+H2OMW))
      U4=-p1*(d_fu*densw/fu +densw/xx(1,j)) -(H2OMW-p1*PKOHMW)*
     &(d_fu*c(1,j)/fu + 1.0D0)

      d(ip2,1)=-frt*c(ki2,j)/p2/p2*p3/2.0d0
      d(ip2,1)=d(ip2,1)+(p1+r1/c_solv)*(d_fu/fu+1/r1)*(1.0d0/h)
      d(ip2,1)=d(ip2,1)+(p1+r1/c_solv)*dcf*(-d_fu/fu/fu*d_fu*0.5d0
     &+0.5d0/fu*d2_fu-0.5d0/r1/r1)
      d(ip2,1)=d(ip2,1)+(d_fu/fu+1/r1)*dcf*(-0.5d0*r1*H2OMW/U1**2.d0*U2
     &+0.5*H2OMW/U1)
      d(ip2,1)=-d(ip2,1)

      b(ip2,1)=frt*c(ki2,j)/p2/p2*p3/2.0d0
      b(ip2,1)=b(ip2,1)+(p1+r1/c_solv)*(d_fu/fu+1/r1)*(1.0d0/h)
      b(ip2,1)=b(ip2,1)+(p1+r1/c_solv)*dcf*(-d_fu/fu/fu*d_fu*0.5d0
     &+0.5d0/fu*d2_fu-0.5d0/r1/r1)
      b(ip2,1)=b(ip2,1)+(d_fu/fu+1/r1)*dcf*(-0.5d0*r1*H2OMW/U1**2.d0*U2
     &+0.5*H2OMW/U1)
      b(ip2,1)=-b(ip2,1)
      d(ip2,kp2)=frt/h
      b(ip2,kp2)=-frt/h
      b(ip2,ki2)=frt/p2
      g(ip2)=-frt*(c(kp2,j+1)-c(kp2,j))/h-frt/p2*c(ki2,j)
     &-(p1+r1/c_solv)*(d_fu*1/fu+1/r1)*dcf

      else ! for Lithium chemistry

      dcf=(xx(1,j+1)-xx(1,j))/h *2.0d0 !factor of 2 added
      r1=(xx(1,j+1)+xx(1,j))/2.0d0
      r4=xx(ki2,j)
      p1=(tm(j)+tm(j+1))/2.0d0
      p2=(fac*cd(j)+cd(j+1))/2.0d0
      p4=(dfu(j)+dfu(j+1))/2.0d0
      d(ip2,1)=(1.0d0-p1)*(1.0d0/r1+p4)/h *2.0d0 !factor of 2 added
      b(ip2,1)=-d(ip2,1)+((1.0d0-p1)*(d2fu(j)-1.0d0/r1/r1)*dcf
     &-(1.0d0/r1+p4)*dcf*dtm(j)+frt*r4*fac*dcd(j)/p2/p2)/2.0d0
      d(ip2,1)=d(ip2,1)+((1.0d0-p1)*(d2fu(j+1)-1.0d0/r1/r1)*dcf
     &-(1.0d0/r1+p4)*dcf*dtm(j+1)+frt*r4*dcd(j+1)/p2/p2)/2.0d0
      d(ip2,kp2)=-frt/h
      b(ip2,kp2)=frt/h
      b(ip2,ki2)=-frt/p2
      g(ip2)=-(1.0d0-p1)*(1.0d0/r1+p4)*dcf+frt*(c(kp2,j+1)-c(kp2,j))/h
     &+frt/p2*c(ki2,j)
      endif

      if (imp.eq.1) g(ip2)=0.0d0
   12 if(imp.eq.1) g(ip2+n)=0.0d0

c     ________________________________________________
c     Equation 3, Butler-Volmer kinetics
c     Note that a zero and a nonzero time step are treated differently.
c     Also note that the equation is written differently depending on
c     whether superposition or a discretized version of the solid-phase
c     material balance is being used.

c Define the individual pore-wall fluxes.  These are necessary for
c multiple particle sizes in order to treat properly a film resistance,
c and for a rigorous energy balance.

      izfl=0
      if (rr.eq.0.0d0) izfl=1 ! need a special flag for zero timestep

      do ii=1,nj ! counter index loop
      do mpa=1,npa
      if (ii.le.n1+1) then
      pwfindiv(mpa,ii)=(ai2(1,mpa)*(xt(2+mpa,ii,kk-1+kadd)-c(2+mpa,ii))
     &/rr*omi-sumpa(mpa,ii))*Rad1pa(mpa)
      if (mvdc1.eq.1) pwfindiv(mpa,ii)=xx(kj,ii)
      pwfindiv_d(mpa,ii)=-Rad1pa(mpa)*ai2(1,mpa)/rr*omi
      elseif (ii.ge.n1+n2) then
      pwfindiv(mpa,ii)=(ai(1,mpa)*(xt(2+mpa,ii,kk-1+kadd)-c(2+mpa,ii))
     &/rr*omi-sumpa(mpa,ii))*Rad3pa(mpa)
      if (mvdc3.eq.1) pwfindiv(mpa,ii)=xx(kj,ii)
      pwfindiv_d(mpa,ii)=-Rad3pa(mpa)*ai(1,mpa)/rr*omi
      else
      pwfindiv(mpa,ii)=0.0d0
      pwfindiv_d(mpa,ii)=0.0d0
      endif
      enddo !mpa
      enddo

      if(rr.eq.0.d0) then ! do a sum of B-V equations to get j
        if(j.gt.n1+1 .and. j.lt.n1+n2) then
          b(ii2,kj)=1.0d0
          g(ii2)=-c(kj,j) ! zero transfer current in separator
        else
          do mpa=1,npa
            call ekin(j,kk,0,0.0d0,mpa,mvdc1,mvdc3)
            if(j.eq.n1+n2 .and. mpa.eq.1) Uc=g0
            if(j.eq.1 .and. mpa.eq.1)     Ua=g0
            if(j.le.n1+1) then
              g(ii2)=g(ii2)+g(2+mpa)*area1pa(mpa)/area1
                do ii=1,n
                  b(ii2,ii)=b(ii2,ii)+b(2+mpa,ii)*area1pa(mpa)/area1
                enddo ! ii
            elseif(j.ge.n1+n2) then
              g(ii2)=g(ii2)+g(2+mpa)*area3pa(mpa)/area3
              do ii=1,n
                b(ii2,ii)=b(ii2,ii)+b(2+mpa,ii)*area3pa(mpa)/area3
              enddo ! ii
            endif
          enddo ! mpa
        endif
        do mpa=1,npa
          do ii=1,n
            b(2+mpa,ii)=0.d0
          enddo
          b(2+mpa,2+mpa)=1.0d0
          g(2+mpa)=xt(2+mpa,j,kk-1+kadd)-c(2+mpa,j) ! fix solid concentration
        enddo ! mpa

      else !rr>0

      if((mvdc1.eq.1.and.j.le.n1+1).or.
     &(mvdc3.eq.1.and.j.ge.n1+n2)) then  ! Variable diffusion coefficient
      mpa=1 ! only one size for variable solid-phase diffusion coefficient
      if(j.gt.n1+1 .and. j.lt.n1+n2) then
      b(2+mpa,kj)=1.0d0
      g(2+mpa)=-c(kj,j) ! zero transfer current in separator
      else
      call ekin(j,kk,0,0.0d0,1,mvdc1,mvdc3)
       if(j.eq.n1+n2) Uc=g0
       if(j.eq.1)     Ua=g0
      endif

      else !using superposition

      do mpa=1,npa
        if(j.gt.n1+1 .and. j.lt.n1+n2) then !separator
          if (imp.eq.0) then
            b(2+mpa,2+mpa)=1.0d0
            g(2+mpa)=-c(2+mpa,j) ! zero solid concentration in separator
          else
            b(2+mpa,kj)=1.0d0
            g(2+mpa)=-c(kj,j)
          endif
        else
          call ekin(j,kk,0,0.0d0,mpa,mvdc1,mvdc3)
          if(j.eq.n1+n2 .and. mpa.eq.1) Uc=g0
          if(j.eq.1 .and. mpa.eq.1)     Ua=g0

          if(j.le.n1+1) then
            g(2+mpa)=g(2+mpa)+xx(kj,j)-Rad1pa(mpa)*(-sumpa(mpa,j)
     &      +ai2(1,mpa)/rr*omi*(xt(2+mpa,j,kk-1+kadd)-c(2+mpa,j)))
            b(2+mpa,2+mpa)=b(2+mpa,2+mpa)-Rad1pa(mpa)*ai2(1,mpa)/rr*omi
            if (imp.eq.0) b(2+mpa,kj)=b(2+mpa,kj)-1.d0
          elseif(j.ge.n1+n2) then
            g(2+mpa)=g(2+mpa)+xx(kj,j)-Rad3pa(mpa)*(-sumpa(mpa,j)
     &      +ai(1,mpa)/rr*omi*(xt(2+mpa,j,kk-1+kadd)-c(2+mpa,j)))
            b(2+mpa,2+mpa)=b(2+mpa,2+mpa)-Rad3pa(mpa)*ai(1,mpa)/rr*omi
            if (imp.eq.0) b(2+mpa,kj)=b(2+mpa,kj)-1.d0
          endif
        endif
      enddo ! mpa
      endif

      endif

      mpa=1 ! only one size for impedance
      if (imp.eq.1) g(2+mpa)=0.0d0
      if(imp.eq.1) g(2+mpa+n)=0.0d0

c ________________________________________________
c    Equation iSk1.  Side reaction 1 kinetics
c    This is set up for oxygen evolution (Tafel) and 
c    recombination (concentration driven) for the NiMH system.
c    We also have a sample side reaction for solvent reduction/oxidation
c    in a Li-ion system.  Tafel kinetics is used at each electrode.
c    For this reaction, we recommend setting rksa1=1.0d-9 and rksac1=1.0d-9
c    in the input file.

      if (nside.ge.1) then

      rksa1_save=rksa1
      rksc1_save=rksc1
      ti0n=dexp((Ebarks1a)*(t-298.0d0)/(t*298.0d0))
      ti0p=dexp((Ebarks1c)*(t-298.0d0)/(t*298.0d0))
      rksa1=rksa1*ti0n
      rksc1=rksc1*ti0p

      alphasO2=1.17d0
      UsO2=0.25d0

      if (j.le.n1+1) then  !recombination at the anode
      g(iSk1)=-c(kj1,j)+rksa1*c(kS1,j)
      b(iSk1,kj1)=1.0d0
      b(iSk1,kS1)=-rksa1

c Sample kinetics for side reaction in a Li-ion system: comment out
c the above and uncomment the below.  We use an alpha value of 0.5
c and use a value of 1.0 V for the "equilibrium" potential of
c solvent reduction

c      g(iSk1)=-c(kj1,j)-rksa1*expg(0.5d0*frt*(c(kp1,j)-c(kp2,j)-1.d0))
c      b(iSk1,kj1)=1.0d0
c      b(iSk1,kp1)=rksa1*0.5d0*frt*expg(0.5d0*frt*(c(kp1,j)-c(kp2,j)
c    &-1.d0))
c      b(iSk1,kp2)=-b(iSk1,kp1)

      elseif (j.ge.n1+n2) then !evolution at the cathode
      g(iSk1)=-c(kj1,j)-rksc1*expg(alphasO2*frt*(c(kp1,j)-c(kp2,j)
     &-UsO2))
      b(iSk1,kj1)=1.0d0
      b(iSk1,kp1)=rksc1*alphasO2*frt*expg(alphasO2*frt*
     &(c(kp1,j)-c(kp2,j)-UsO2))
      b(iSk1,kp2)=-b(iSk1,kp1)

c Sample kinetics for side reaction in a Li-ion system: comment out
c the above and uncomment the below.  We use an alpha value of 0.5
c and use a value of 4.4 V for the "equilibrium" potential of 
c solvent oxidation

c      g(iSk1)=-c(kj1,j)-rksc1*expg(0.5d0*frt*(c(kp1,j)-c(kp2,j)-4.d0))
c      b(iSk1,kj1)=1.0d0
c      b(iSk1,kp1)=rksc1*0.5d0*frt*expg(0.5d0*frt*(c(kp1,j)-
c     1c(kp2,j)-4.d0))
c      b(iSk1,kp2)=-b(iSk1,kp1)

      else !in the separator
      g(iSk1)=c(kj1,j)
      b(iSk1,kj1)=-1.0d0
      endif

      else
      g(iSk1)=0.0d0
      endif

      rksa1=rksa1_save
      rksc1=rksc1_save

c ________________________________________________
c    Equation iSk2.  Side reaction 2 kinetics 
c    This is set up for hydrogen evolution from the negative electrode
c    of the NiMH system.

      if (nside.ge.2) then
      alphasH2=1.0d0
      UsH2=-0.98d0

      rksa2_save=rksa2
      rksc2_save=rksc2
      ti0n=dexp((Ebarks3a)*(t-298.0d0)/(t*298.0d0))
      ti0p=dexp((Ebarks3c)*(t-298.0d0)/(t*298.0d0))
      rksa2=rksa2*ti0n
      rksc2=rksc2*ti0p

      if (j.le.n1+1) then  
c H2 evolution from the negative electrode only.
      g(iSk2)=c(kj2,j)-rksa2*expg(-alphasH2*frt*(c(kp1,j)
     &-c(kp2,j)-UsH2))
      b(iSk2,kj2)=-1.0d0
      b(iSk2,kp1)=-rksa2*alphasH2*frt*
     &expg(-alphasH2*frt*(c(kp1,j)-c(kp2,j)-UsH2))
      b(iSk2,kp2)=-b(iSk2,kp1)
      else !in the separator or the positive
      g(iSk2)=c(kj2,j) 
      b(iSk2,kj2)=-1.0d0
      endif

      else
      if(imp.eq.1) g(3+n)=0.0d0
      endif

      rksa2=rksa2_save
      rksc2=rksc2_save

c ________________________________________________
c    Equation iSk3.  Side reaction 3 kinetics
c    This is set up for hydrogen absorption/desorption in the
c    negative electrode of the NiMH system.
c Positive j is H leaving the solid phase

      if (nside.eq.3) then

      rksa3_save=rksa3
      rksc3_save=rksc3
      ti0e=dexp((Ebarks2a)*(t-298.0d0)/(t*298.0d0))
      ti0r=dexp((Ebarks2c)*(t-298.0d0)/(t*298.0d0))
      rksa3=rksa3*ti0e
      rksc3=rksc3*ti0r

      mpa=1 ! only one size with side reactions
      if (j.le.n1+1) then  
c H2 absorption, using an equilibrium isotherm from the literature

      g(iSk3)=c(kj3,j)+rksa3*(c(kS2,j)*t*8.206d-5-
     &(5.d0+4.9999d0*dtanh((c(2+mpa,j)/ct1-1.04d0)/0.11d0)))
      b(iSk3,kS2)=-rksa3*t*8.2d-5
      b(iSk3,2+mpa)=-rksa3*4.9999d0*
     &((dcosh((c(2+mpa,j)/ct1-1.04d0)/0.11d0))**(-2))*1.0d0/0.11d0/ct1
      b(iSk3,kj3)=-1.0d0

      elseif (j.ge.n1+n2) then !reaction at the positive electrode

      g(iSk3)=c(kj3,j)
      b(iSk3,kj3)=-1.0d0

      else !in the separator
      g(iSk3)=c(kj3,j)
      b(iSk3,kj3)=-1.0d0
      endif

      else
      mpa=1 ! only one size for impedance
      if(imp.eq.1) g(5+n)=0.0d0
      endif

      rksa3=rksa3_save
      rksc3=rksc3_save

c     ________________________________________________
c     Equation ip1, Ohm's law in solid
c
      if(j.eq.1) then
      if (n1 .eq. 0) then ! foil electrode
      g(ip1) = cur - fc*c(kj,j)
      b(ip1,kj) = fc
      if(imp.eq.1) then
      g(ip1) = 1.d0
c     area1 is mult by L here so aL is 1 in double-layer charging term:
      b(ip1,kp2+n)=omega(ji)*capp1
      b(ip1,kp1+n)=-omega(ji)*capp1
      b(ip1+n,kp2)=-omega(ji)*capp1
      b(ip1+n,kp1)=omega(ji)*capp1
      endif ! impedance section (j=1, a)
      else
      b(ip1,kp1)=-1.0d0/h1 ! not a foil electrode
      d(ip1,kp1)= 1.0d0/h1
      b(ip1,ki2)=-1.0d0/sig1
c     next statement requires that cur be the current in the separator
      g(ip1)=-cur/sig1+c(ki2,j)/sig1-(c(kp1,j+1)-c(kp1,j))/h1
      if(imp.eq.1) g(ip1)=-1.0d0/sig1
      endif
      elseif(j.lt.n1+1) then ! anode
      b(ip1,kp1)=-1.0d0/h1
      d(ip1,kp1)= 1.0d0/h1
      b(ip1,ki2)=-1.0d0/sig1
      g(ip1)=-cur/sig1+c(ki2,j)/sig1-(c(kp1,j+1)-c(kp1,j))/h1
      if(imp.eq.1) g(ip1)=-1.0d0/sig1
      elseif(j.eq.n1+1) then
c     This is the current boundary condition.
            if(mcL.eq.-3) then ! constant load
      g(ip1)=(xx(nx+1,j+1)-xx(nx+1,j))/xx(ki2,j)-cuL-RG !WHT 1-11-08
      d(ip1,nx+1)=-1.d0/xx(ki2,j)
      b(ip1,nx+1)=1.d0/xx(ki2,j)
      b(ip1,ki2)=(xx(nx+1,j+1)-xx(nx+1,j))/xx(ki2,j)**2
            elseif(mcL.eq.0) then ! constant potential
      g(ip1)=xx(nx+1,j+1)-xx(nx+1,j)-cuL-RG*xx(ki2,j) !WHT 1-11-08
      d(ip1,nx+1)=-1.d0
      b(ip1,nx+1)=1.d0
      b(ip1,ki2)=RG !WHT 1-11-08
            elseif(mcL.eq.-2) then ! constant power
      g(ip1)=(xx(nx+1,j+1)-xx(nx+1,j))*xx(ki2,j)-cuL-RG*xx(ki2,j)**2 !WHT 1-11-08
      d(ip1,nx+1)=-xx(ki2,j)
      b(ip1,nx+1)=xx(ki2,j)
      b(ip1,ki2)=-(xx(nx+1,j+1)-xx(nx+1,j))+2.d0*RG*xx(ki2,j) !WHT 1-11-08
            elseif(mcL.eq.-4) then ! maximum power
      g(ip1)=xx(nx+1,j+1)-xx(nx+1,j)-xx(ki2,j)*Rint !RG*something
      d(ip1,nx+1)=-1.d0
      b(ip1,nx+1)=1.d0
      b(ip1,ki2)=Rint
      cu0=(xx(nx+1,j+1)-xx(nx+1,j))*xx(ki2,j)
            else ! constant current
      b(ip1,ki2)=1.0d0
      g(ip1)=cur-xx(ki2,j)
            endif

      if(imp.eq.1) g(ip1)=1.0d0
      elseif(j.lt.n1+n2) then
      b(ip1,kp1)=1.0d0 ! separator, set phi1=0
      g(ip1)=-c(kp1,j)
      if(imp.eq.1) g(ip1)=0.0d0
      elseif(j.lt.nj) then ! cathode
      d(ip1,kp1)=1.0d0/h3
      b(ip1,kp1)=-1.0d0/h3
      b(ip1,ki2)=-1.0d0/sig3
      g(ip1)=-cur/sig3+c(ki2,j)/sig3-(c(kp1,j+1)-c(kp1,j))/h3
      if(imp.eq.1) g(ip1)=-1.0d0/sig3
      else ! j=nj
      b(ip1,ki2)=1.0d0
      g(ip1)=-c(ki2,j) !  i2 is no longer used at j=nj
      if(imp.eq.1) g(ip1)=0.d0
      endif
      if(imp.eq.1) g(12)=0.0d0

c     ________________________________________________

c     Equation n+1.  Carry end potentials.  This equation is added to
c     permit handling constant load, etc. without an extra loop.
      if (imp.eq.0) then
      if(j.eq.1) then
      g(nx+1)=xx(kp1,1)-xx(nx+1,1) ! pick up negative potential
      b(nx+1,kp1)=-1.d0
      b(nx+1,nx+1)=1.d0
      elseif(j.eq.nj) then
      g(nx+1)=xx(kp1,nj)-xx(nx+1,nj) ! pick up positive potential
      b(nx+1,kp1)=-1.d0
      b(nx+1,nx+1)=1.d0
      elseif(j.le.n1+1) then
      g(nx+1)=xx(nx+1,j)-xx(nx+1,j-1) ! carry negative potential
      b(nx+1,nx+1)=-1.d0
      a(nx+1,nx+1)=1.d0
      else ! if(j.ge.n1+n2) then
      g(nx+1)=xx(nx+1,j)-xx(nx+1,j+1) ! carry positive potential
      b(nx+1,nx+1)=-1.d0
      d(nx+1,nx+1)=1.d0
      endif
      endif

c     ________________________________________________
c
c     Equation i2div.  Current balance

      if(j.gt.n1+1 .and. j.lt.n1+n2) then ! separator
        g(i2div)=xx(ki2,j-1)-xx(ki2,j)
        if (imp.eq.0) a(i2div,ki2)=-1.d0
        b(i2div,ki2)=1.d0
        if (imp.eq.1) then ! impedance section
c     Set perturbation current = 1 purely real:
          b(i2div,ki2)=1.0d0
          g(i2div)=1.0d0
        endif ! end of impedance section (n1+1<j<n1+n2)
c
      elseif(j.eq.1 .and. n1.eq.0) then ! foil electrode
        g(i2div) = cur - c(ki2,j)
        b(i2div,ki2) = 1.0d0
        if(imp.eq.1) g(i2div)=1.d0 ! set perturbation current = 1 purely real

      else  !no impedance, no foil
        if(cap.eq.0.d0) then
          if (nside.ge.1) then ! Side reactions, both 1,2,3
            b(i2div,ki2)=-1.0d0/hC
            a(i2div,ki2)=1.0d0/hC
            b(i2div,kj)=area*fc
            b(i2div,kj1)=-area*fc
            b(i2div,kj2)=-area*fc
            g(i2div)=c(ki2,j)/hC-area*fc*(c(kj,j)-c(kj1,j)-c(kj2,j)) !not O(h2)
c written only for j=1
            if(j.gt.1) g(i2div)=(c(ki2,j)-c(ki2,j-1))/hC-
     &    area*fc*(c(kj,j)-c(kj1,j)-c(kj2,j))
            if (nside.eq.1) b(i2div,kj2)=0.0d0
          else ! No side reactions
            b(i2div,ki2)=-1.0d0/hC
            a(i2div,ki2)=1.0d0/hC
            b(i2div,kj)=area*fc
            g(i2div)=c(ki2,j)/hC-area*fc*c(kj,j) ! not order h2
c     need to include mpa ! problem
         if(j.gt.1) g(i2div)=(c(ki2,j)-c(ki2,j-1))/hC-area*fc*c(kj,j) !not O(h2)
          endif
        else !nonzero value of capacitances
        if (rr .eq. 0.d0) then ! no instantaneous change across double layer
          g(i2div)=c(kp1,j)-c(kp2,j)-xt(kp1,j,kk-1+kadd)+
     &    xt(kp2,j,kk-1+kadd)
          b(i2div,kp1)=-1.d0
          b(i2div,kp2)= 1.d0
        else !nonzero capacitance and nonzero time step
          b(i2div,ki2)=-1.0d0/hC
          a(i2div,ki2)=1.0d0/hC
          b(i2div,kj)=area*fc
          b(i2div,kj1)=area*fc
          b(i2div,kj2)=area*fc
          b(i2div,kp1)=area*cap/rr*omi*2.0d0
          b(i2div,kp2)=-area*cap/rr*omi*2.0d0
          g(i2div)=(c(ki2,j)+xt(ki2,j,kk-1+kadd))/hC
          if(j.gt.1) g(i2div)=g(i2div)+(-c(ki2,j-1)-xt(ki2,j-1,
     &  kk-1+kadd))/hC
          g(i2div)=g(i2div) -area*fc*(c(kj,j)+xt(kj,j,kk-1+kadd)
     &    -c(kj1,j)-xt(kj1,j,kk-1+kadd)-c(kj2,j)-xt(kj2,j,kk-1+kadd))
     &    -area*cap*(c(kp1,j)-c(kp2,j)-xt(kp1,j,kk-1+kadd)
     &    +xt(kp2,j,kk-1+kadd))/rr*omi*2.0d0
        endif
      endif
c
      if (imp.eq.1) then ! impedance section
c     Add double-layer charging:
      b(i2div,kp2+n)=omega(ji)*area*capp
      b(i2div,kp1+n)=-omega(ji)*area*capp
      b(kj+n,kp2)=-omega(ji)*area*capp
      b(kj+n,kp1)=omega(ji)*area*capp
      g(i2div)=0.0d0
      g(kj+n)=0.0d0
      endif ! end of impedance section (Equation i2div)
      endif

c     ________________________________________________
c
  121 if (imp.eq.1) then ! impedance section
c     Copy real matrices into imag matrices:
      do 122 nct1=1, n
      do 122 nct2=1, n
      if(j.eq.1) x(nct1+kp1,nct2+kp1)=x(nct1,nct2)
      if(j.gt.1) a(nct1+kp1,nct2+kp1)=a(nct1,nct2)
      b(nct1+kp1,nct2+kp1)=b(nct1,nct2)
      if(j.eq.nj) y(nct1+kp1,nct2+6)=y(nct1,nct2)
      if(j.lt.nj) d(nct1+6,nct2+6)=d(nct1,nct2)
  122 continue
      endif
c
      if(imp.eq.0 .and. mcL.ne.1) nx=nx+1 !carry end potentials only if needed

c     double the number of variables (real & imaginary parts)
      if(imp.eq.1) n=2*n
            call band(j)
      if(imp.eq.1) n=n/2

      if(imp.eq.0 .and. mcL.ne.1) nx=nx-1

      if(j.lt.nj) go to 10

      if (imp.eq.0) then
      do 607 jj=1,nj
      do 607 i=1,nx+1
  607 c(i,jj)=xx(i,jj)+c(i,jj)  
      endif ! make corrections for imp.ne.1

c     ________________________________________________
c
c     print impedance results
c
      if (imp.eq.1) then
      vreal(ji)=c(kp1,1)-c(kp1,nj)+RG !WHT 1-11-08
      vimag(ji)=c(kp1+n,1)-c(kp1+n,nj)
c     Half-cell impedance calculations:
      vreala(ji)=c(kp1,1)-c(kp2,n1+n2)+RGn !WHT 1-11-08
      vrealc(ji)=c(kp2,n1+n2)-c(kp1,nj)+RGp !WHT 1-11-08
      vimaga(ji)=c(kp1+n,1)-c(kp2+n,n1+n2)
      vimagc(ji)=c(kp2+n,n1+n2)-c(kp1+n,nj)
c     Z magnitude and phase angle calculation:
      zmag(ji)=(vreal(ji)*vreal(ji)+vimag(ji)*vimag(ji))**0.5d0
      zmaga(ji)=(vreala(ji)*vreala(ji)+vimaga(ji)*vimaga(ji))**0.5d0
      zmagc(ji)=(vrealc(ji)*vrealc(ji)+vimagc(ji)*vimagc(ji))**0.5d0
      phi(ji)=180.0d0/pi*datan(- vimag(ji)/vreal(ji))
      phia(ji)=180.0d0/pi*datan(- vimaga(ji)/vreala(ji))
      phic(ji)=180.0d0/pi*datan(- vimagc(ji)/vrealc(ji))

      write (2,103)omega(ji),(10**4)*vreal(ji),(10**4)*vimag(ji)  !Jun29

      rr=0.01 !to compensate for rr=1 set to avoid div by 0
      ji=ji+1
      return
      endif ! end of impedance section (printing)
c
      do 56 i=1,n
   56 xp(i)=(4.0d0*c(i,2)-3.0d0*c(i,1)-c(i,3))/2.0d0/h1
      nerr=0
      nconcflag = 0
      do 25 j=1,nj
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c     shoe horns:

      if(c(1,j).lt.xx(1,j)/100.d0) c(1,j)=xx(1,j)/100.d0
      if(c(1,j).gt.xx(1,j)*100.d0) c(1,j)=xx(1,j)*100.0d0
      if (c(kp2,j).lt.(xx(kp2,j)-0.02d0)) c(kp2,j)=xx(kp2,j)-0.02d0
      if (c(kp2,j).gt.(xx(kp2,j)+0.02d0)) c(kp2,j)=xx(kp2,j)+0.02d0
      if (c(kp1,j).lt.(xx(kp1,j)-0.02d0)) c(kp1,j)=xx(kp1,j)-0.02d0
      if (c(kp1,j).gt.(xx(kp1,j)+0.02d0)) c(kp1,j)=xx(kp1,j)+0.02d0
      if (c(kp2,j).gt. 9.9d0) c(kp2,j)= 9.9d0
      if (c(kp2,j).lt.-9.9d0) c(kp2,j)=-9.9d0
      if (c(kp1,j).gt. 9.9d0) c(kp1,j)= 9.9d0
      if (c(kp1,j).lt.-9.9d0) c(kp1,j)=-9.9d0
c
      do mpa=1,npa
      if (j .ge. n1+n2) then
      if(c(2+mpa,j).lt.xx(2+mpa,j)/100.d0) then
      nerr=nerr+1
      c(2+mpa,j)=xx(2+mpa,j)/100.d0 ! use cs min
      endif
      if(c(2+mpa,j).gt.xx(2+mpa,j)*100.d0) c(2+mpa,j)=xx(2+mpa,j)*100.d0 
      if(ct3-c(2+mpa,j).le.(ct3-xx(2+mpa,j))/100.d0) then
      nerr=nerr+1
      c(2+mpa,j)=ct3-(ct3-xx(2+mpa,j))/100.d0
      endif
      if(c(2+mpa,j).ge.ct3) c(2+mpa,j)=0.999999d0*ct3
c
      else if (j .le. n1+1 .and. n1 .gt. 0) then

      if(c(2+mpa,j).lt.xx(2+mpa,j)/100.d0) nerr=nerr+1
      if(c(2+mpa,j).lt.xx(2+mpa,j)/100.d0) c(2+mpa,j)=xx(2+mpa,j)/100.d0 ! use cs min
      if(c(2+mpa,j).gt.xx(2+mpa,j)*100.d0) c(2+mpa,j)=xx(2+mpa,j)*100.d0 
      if(ct1-c(2+mpa,j).le.(ct1-xx(2+mpa,j))/100.d0) nerr=nerr+1
      if(ct1-c(2+mpa,j).le.(ct1-xx(2+mpa,j))/100.d0)
     & c(2+mpa,j)=ct1-(ct1-xx(2+mpa,j))/100.d0
      if(c(2+mpa,j).ge.ct1) c(2+mpa,j)=0.999999d0*ct1
      endif
c
c     Not shoe-horns, but trips to force smaller time steps
      if (j. ge .n1+n2 .or. (j .le. n1+1 .and. n1 .gt. 0)) then
      if(c(2+mpa,j).lt.1.0d-10) c(2+mpa,j) = 1.d-10
      if(c(2+mpa,j).lt.1.0d-3) nconcflag = 1
      endif
      enddo ! mpa
c
      if(c(1,j).lt.1.0d-10) then
      c(1,j) = 1.d-10
c      c(kj,j) = 0.0d0
      endif
      if(c(1,j).lt.1.0d-03) nconcflag = 2
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
      c(nx+1,j)=c(ip1,nj)
      if(j.le.n1+1) c(nx+1,j)=c(ip1,1)
      do 25 i=1,n+1
   25 xx(i,j)=c(i,j) 

      if (jcount .gt. 3*lim .and. rr.eq.0.0d0) then
      write (2,99)

      stop
      endif
c     ________________________________________________
c
c     Decreasing time steps:
      if (jcount .gt. lim .and. rr.gt.1.0d-16) then
        tau=tau/2.0d0
        rr=tau
        Lcount=0
        if(mcLL.eq.-4) cur=curold
        ts(kk)=ts(kk-1)+tau
c     print *,'time step reduced to ', tau, ts(kk),jcount
        if (tau.lt.1.0d-4) then
      nconv=1  !trigger the flag for no convergence
          if (lpow.eq.1) then  !peak power activated
            nflag=1
            go to 66
          endif
          nt=nt-1
          pow=ed/tw/ts(nt+1)
      write (2,*) 'mass is ',tw
      write (2,*) 'energy is ',ed/tw/3.6d03
      write (2,*) 'power is ',pow
          write (2,*) kk-1,' this time step did not converge'
          if (nconcflag .eq. 1) then
            write (2,*), 'Solid concentration driven to zero.'
            write (2,*),'Consider changing vcut or the current density.'
            if (nconcflag.eq.3) then
              write (2,*) 'The solid utilization is above 1.'
            endif
          endif
          if (nconcflag .eq. 2) then
            write (2,*), 'Salt concentration driven to zero.'
            write (2,*),'Consider changing vcut or the current density.'
      endif
      write (2,*) ' '
      if (kk.gt.5) then
      write (2,*) 'dV/dt at point of failure is (V/s)', 
     &    ((xt(kp1,nj,kk-1)-xt(kp1,1,kk-1))-(xt(kp1,nj,kk-2)
     &    -xt(kp1,1,kk-2)))/(ts(kk-1)-ts(kk-2))
      write (2,*) 'Moving on to next leg if not the last'
      write (2,*) ' '
      endif
      return
          else
            iflag=0
            call calca(kk)
            go to 666
          end if
c
      else
c     ________________________________________________
c
c check for convergence

c
      if(nerr.ne.0) go to 8
      do jj=1,nj
      do ki=1,n
      errlim=1.d-10
      if(ki.eq.kj) errlim=1.d-16 !change to -14 if problems with convergence
      dxx=dabs( xx(ki,jj) - xx0(ki,jj) )
      enddo
      enddo
      do 55 ki=1,n
      errlim=1.d-10
      if(ki.eq.kj) errlim=1.d-16 !change to -14 if problems with convergence
      dxp=dabs( xp(ki)-xp0(ki) )
      n1hold = 1
      if (n1 .ge. 11) n1hold = n1-10
      dxx=dabs( xx(ki,n1hold) - xx0(ki,n1hold) )
      dxx2=dabs( xx(ki,n1+n2+10)-xx0(ki,n1+n2+10) )
      dxx3=dabs( xx(ki,2)-xx0(ki,2) )
      if(dxx .gt. 1.d-9*dabs(xx(ki,n1hold)).and.dxx.gt.errlim) then
      go to 8
      endif
      if(dxx2.gt.1.d-9*dabs(xx(ki,n1+n2+10)).and.dxx2.gt.errlim) then
      go to 8
      endif
      if(dxx3.gt.1.d-9*dabs(xx(ki,2)).and.dxx3.gt.errlim) then
      go to 8
      endif
   55 continue

      if(mcLL.eq.-4) then ! maximum power
      vvold=vv
      vv=xx(kp1,nj)-xx(kp1,1)
      pwr=cur*vv


      if(Lcount.eq.1) then
        curlow=0.d0
        curhigh=0.d0
        curmaxx=cur
        pwrmaxx=cur*vv
      else

      if(cur*vv.gt.pwrmaxx) then
        if(curmaxx.lt.cur) then
          curlow=curmaxx
        else
          if(Lcount.gt.2) curhigh=curmaxx
        endif
        curmaxx=cur
        pwrmaxx=cur*vv

      else ! no new power maximum
        if(cur.lt.curmaxx) then
          curlow=cur
        elseif(cur.lt.curhigh) then
          curhigh=cur
        endif
      endif

c     May have a problem with Rint if the difference between
c     cur and curold becomes too small.  However, the bounds
c     from curlow and curhigh should still give valid answers.
      Rint=-(vv-vvold)/(cur-curold)+RG ! internal resistance
      OCP=vv+(Rint-RG)*cur
      curmax=OCP/Rint/2.d0 ! estimate of current at max power
      pwrmax=OCP**2/Rint/4.d0 ! estimate of max power
      if(Lcount.eq.2) curhigh=2.d0*curmax
      endif
      curold=cur

      cur=cur-1.d0
      if(Lcount.gt.1) then
      if(curmax.gt.1.2d0*cur) curmax=1.2d0*cur
      cur=curmax
      if(cur.lt.curlow .or. cur.gt.curhigh) cur=0.5d0*(curlow+curhigh)
      endif
      if(Lcount.lt.18 .and. dabs(cur-curold).gt.1.d-5*dabs(cur)) go to 3
      if(Lcount.lt.5) go to 3
      cur=curold
      endif

      if ((xx(1,nj)/xt(1,nj,kk-1)).lt.0.75d0.and.xx(1,nj).lt.100.0d0)
     1rr=rr/2.0d0

c     ________________________________________________
c
c     if(lpow.ne.1) write (2,*) jcount,' iterations required'
c
      do 60 ll=1, nj  ! save present time results in xt()
      do 60 lk=1,n
   60 xt(lk,ll,kk)=xx(lk,ll)

      do 61 ll=1,nj
      do 61 lk=1,nnj
   61 cssold(ll,lk)=css(ll,lk)


c
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if(rr.ne.0.0d0) then
      do 58 j=1,nj  ! nonzero time step
      term_s2(j)=termn_s2(j)
      term_s1(j)=termn_s1(j)
   58 term(j)=termn(j)
      else ! zero time step
      do 65 j=1,nj
      term(j)=0.d0
      term_s1(j)=0.0d0
      term_s2(j)=0.0d0

      fac=1.d0
      if(j.eq.n1+2 .and. n1 .gt. 0)
     &fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2+1) fac=((ep3+epp3)/(ep2+epp2))**exbrug

      hn=h2 
      if(n1.gt.0 .and. j.le.n1+1) hn=h1
      if(n1.gt.0 .and. j.le.n1+1) area=area1
      if(j.gt.n1+n2) hn=h3
      if(j.gt.n1+n2) area=area3

      if(j.gt.1) term(j)=
     &-(df(j)+fac*df(j-1))*(c(1,j)-c(1,j-1))/hn/2.d0
     &-(1.d0-0.5d0*(tm(j)+tm(j-1)))*c(ki2,j-1)/fc
      if (j.gt.1) term_s1(j)=
     &-(do2(j)+fac*do2(j-1))*(c(kS1,j)-c(kS1,j-1))/hn/2.0d0
     &+hn*area*vol*(0.375d0*c(kj1,j)+0.125d0*c(kj1,j-1))
      if (j.gt.1) term_s2(j)=
     &-(dh2(j)+fac*dh2(j-1))*(c(kS2,j)-c(kS2,j-1))/hn/2.0d0
     &+hn*area*vol*(0.375d0*(c(kj3,j)+c(kj2,j))+0.125d0*
     &(c(kj3,j-1)+c(kj2,j-1)))

      fac=1.d0
      if(j.eq.n1+1 .and. n1 .gt. 0)
     &fac=((ep2+epp2)/(ep1+epp1))**exbrug
      if(j.eq.n1+n2) fac=((ep3+epp3)/(ep2+epp2))**exbrug

      hn=h2
      if(n1.gt.0 .and. j.lt.n1+1) hn=h1
      if(n1.gt.0 .and. j.lt.n1+1) area=area1
      if(j.ge.n1+n2) hn=h3
      if(j.ge.n1+n2) area=area3

      if(j.lt.nj) term(j)=term(j)
     &-(fac*df(j)+df(j+1))*(c(1,j)-c(1,j+1))/hn/2.d0
     &+(1.d0-0.5d0*(tm(j)+tm(j+1)))*c(ki2,j)/fc
      if (j.lt.nj) term_s2(j)=term_s2(j)
     &-(fac*dh2(j)+dh2(j+1))*(c(kS2,j)-c(kS2,j+1))/hn/2.0d0
     &+hn*area*vol*(0.375d0*(c(kj3,j)+c(kj2,j))+0.125d0*
     &(c(kj3,j+1)+c(kj2,j+1)))
   65 if (j.lt.nj) term_s1(j)=term_s1(j)
     &-(fac*do2(j)+do2(j+1))*(c(kS1,j)-c(kS1,j+1))/hn/2.0d0
     &+hn*area*vol*(0.375d0*c(kj1,j)+0.125d0*c(kj1,j+1))
      endif
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
      end if
c
   66 continue

      return
      end
c
c***********************************************************************
      subroutine calca(kk)
c     subroutine to calculate diffusion in solid particles
c     of active material, based on a superposition integral.
      implicit real*8(a-h,o-z)
      parameter(maxt=900)
      common /n/ tmmax,imp,ji,nx,nt,n1,n2,nj,n3,nconv,npa,iSk3,kj3
     &,ii2,ki2,kj,i2div,ip2,kp2,ip1,kp1,imb2,kS2,imb1,kS1,iSk1,kj1
     &,iSk2,kj2
      common/activ/ EbarD,Ebarkap,Ebarka,Ebarkc,Ebarr1,Ebarr3,
     &Ebars1,Ebars3,Ebarks1a,Ebarks1c,Ebarks2a,Ebarks2c,
     &Ebarks3a,Ebarks3c
      common /calc/ ai(maxt,5),ai2(maxt,5),ts(maxt),h,h1,h2,h3,hcn,
     1hcp,rr,rrmax,cuL,sumpa(5,221),Rad1pa(5),Rad3pa(5),area1pa(5)
     &,area3pa(5),mcL
      common/pindiv/pwfindiv(5,221),pwfindiv_d(5,221),vf1(5),vf3(5),izfl
      common/const/ fc,r,t,frt,cur,ep3,ep2,pi,ep1,epf3,epf1,
     &epp1,epp2,epp3,shape3,shape1,capp1,capp3,nneg,nprop,npos
      common/var/ xp(16),xx(17,221),xt(16,221,maxt)
      common/gas/ epg1,epg2,epg3
      common/cprop/ sig3,area3,rka3,rka3save,ct3,dfs3,Rad3,cap1,cap3,
     1sig1,area1,rka1,rka1save,ct1,dfs1,Rad1,tw,dfs1save,dfs3save
      dimension ar(4,maxt),bz(6)
c
      do mpa=1,npa ! different particle sizes
      do 319 l=1,nt
      ai2(l,mpa)=0.0d0
  319 ai(l,mpa)=0.0d0
c
      do 70 i=1,kk-1
      ar(1,i)=dfs3*(ts(kk)-ts(i))/Rad3pa(mpa)/Rad3pa(mpa)
      ar(2,i)=dfs3*(ts(kk)-ts(i+1))/Rad3pa(mpa)/Rad3pa(mpa)
      ar(3,i)=dfs1*(ts(kk)-ts(i))/Rad1pa(mpa)/Rad1pa(mpa)
      ar(4,i)=dfs1*(ts(kk)-ts(i+1))/Rad1pa(mpa)/Rad1pa(mpa)
c
      do 69 m=1,2
      t1=ar(m,i)
      t2=ar((m+2),i)
c
      a1=0.0d0
      a12=0.0d0
c
      s=1.644934066848d0
c
c     Bessel's function zeros:
      bz(1)=2.4048255577d0
      bz(2)=5.5200781103d0
      bz(3)=8.6537281103d0
      bz(4)=11.7915344391d0
      bz(5)=14.9309177086d0
c
      if (shape3.gt.2.0d0) then
c     spherical particles:
      if (t1 .gt. 0.06d0) then
c
      do 59 j=1,5
      y1=dble(j*j)*pi*pi*t1
   59 if (y1 .le. 1.5d02) a1=a1+(expg(-y1))/dble(j*j)
      a1=2.0d0*(s-a1)/pi/pi
c
      else
c
      if (t1.LE.0.0d0) then
      a1=0.0d0
      else
      do 60 j=1,3
      z=dble(j)/dsqrt(t1)
      call erfcg(z,e)
      y2=dble(j*j)/t1
      if(y2 .ge. 1.5d02) then
      da=-dble(j)*dsqrt(pi/t1)*e
      else
      da=expg(-y2)-dble(j)*dsqrt(pi/t1)*e
      end if
   60 a1=a1+da
      a1=-t1 + 2.0d0*dsqrt(t1/pi)*(1.0d0+2.0d0*a1)
      end if
c
      end if
      else
c
      if (shape3.lt.2.0d0) then
c     planar particles:
      if(t1 .gt. 0.06d0) then
c
      do 61 j=1,5
      da=((-1.0d0)**j)*(1.0d0 - expg(-dble(2*j+1)**2
     1*pi*pi*t1))/dble(2*j+1)**2
   61 a1=a1+da
      a1=4.0d0*a1/pi/pi
c
      else
c
      do 62 j=1,3
      z=dble(j)/2.0d0/dsqrt(t1)
      call erfcg(z,e)
      da=((-1.0d0)**j)*(expg(-dble(j*j)/4.0d0/t1)
     &-dble(j)/2.0d0*dsqrt(pi/t1)*e)
   62 a1=a1+da
      a1=2.0d0*dsqrt(t1/pi)*(1.0d0+2.0d0*a1)
c
      end if
      else ! shape = 2
c     cylindrical particles:
      if (t1.gt.0.06d0) then
c
      do 63 j=1,5
      da=(1.0d0-expg(-bz(j)*bz(j)*t1))/bz(j)/bz(j)
   63 a1=a1+da
      a1=2.0d0*a1
c
      else
c
      a1=2.0d0*dsqrt(t1/pi)-t1/4.0d0-5.0d0*(t1**1.5d0)/96.0d0/dsqrt(pi)
     1-31.0d0*t1*t1/2048.0d0
c
      end if
      end if
      end if
c
      if (n1 .eq. 0) go to 36
c     (skip calculations of Li diffusion in the solid
c     if the negative electrode is metal foil)
c
      if (shape1.gt.2.0d0) then
c     spherical particles:
      if(t2 .gt. 0.06d0) then
c
      do 64 j=1,5
      y2=dble(j*j)*pi*pi*t2
   64 if(y2 .le. 1.5d02) a12=a12+(expg(-y2))/dble(j*j)
      a12=2.0d0*(s-a12)/pi/pi
c
      else
c
      if (t2.eq.0.0d0) then
      a12=0.0d0
      else
      do 65 j=1,3
      z=dble(j)/dsqrt(t2)
      call erfcg(z,e)
      y2=dble(j*j)/t2
      if(y2 .gt. 1.5d02) then
      da=-dble(j)*dsqrt(pi/t2)*e
      else
      da=expg(-y2)-dble(j)*dsqrt(pi/t2)*e
      end if
   65 a12=a12+da
      a12=-t2 + 2.0d0*dsqrt(t2/pi)*(1.0d0+2.0d0*a12)
      end if
      end if
c
      else
      if (shape1.lt.2.0d0) then
c     planar particles:
      if(t2 .gt. 0.06d0) then
c
      do 66 j=1,5
      da=((-1.0d0)**j)*(1.0d0 - expg(-dble(2*j+1)**2
     1*pi*pi*t2))/dble(2*j+1)**2
   66 a12=a12+da
      a12=4.0d0*a12/pi/pi
c
      else
c
      do 67 j=1,3
      z=dble(j)/2.0d0/dsqrt(t2)
      call erfcg(z,e)
      da=((-1.0d0)**j)*(expg(-dble(j*j)/4.0d0/t2)
     &-dble(j)/2.0d0*dsqrt(pi/t2)*e)
   67 a12=a12+da
      a12=2.0d0*dsqrt(t2/pi)*(1.0d0+2.0d0*a12)
c
      end if
      else
c     cylindrical particles:
      if (t2.gt.0.06d0) then
c
      do 68 j=1,5
      da=(1.0d0-expg(-bz(j)*bz(j)*t2))/bz(j)/bz(j)
   68 a12=a12+da
      a12=2.0d0*a12
c
      else
c
      a12=2.0d0*dsqrt(t2/pi)-t2/4.0d0-5.0d0*(t2**1.5d0)/96.0d0
     &/dsqrt(pi)-31.0d0*t2*t2/2048.0d0
c
      end if
      end if
      end if
c
   36 continue
c
      ar(m,i)=a1
   69 ar((m+2),i)=a12
c
      ai(kk-i,mpa)=ar(1,i)-ar(2,i)
   70 ai2(kk-i,mpa)=ar(3,i)-ar(4,i)
c
      do j=1,nj
      sumpa(mpa,j)=0.d0
      enddo
      do j=1,n1+1
      if (kk.gt.2) then  
      do  54 i=1, kk-2
   54 if(ts(i+1)-ts(i).ne.0.0d0) sumpa(mpa,j)=sumpa(mpa,j)
     &+(xt(2+mpa,j,i+1)- xt(2+mpa,j,i))*ai2(kk-i,mpa)/(ts(i+1)-ts(i))
      endif
      enddo
      do j=n1+n2,nj
      if (kk.gt.2) then  
      do 95 i=1, kk-2
   95 if(ts(i+1)-ts(i).ne.0.0d0) sumpa(mpa,j)=sumpa(mpa,j) ! problem
     &+(xt(2+mpa,j,i+1)- xt(2+mpa,j,i))*ai(kk-i,mpa)/(ts(i+1)-ts(i))
      endif
      enddo
      enddo ! mpa
c
      return
      end
c***********************************************************************
      subroutine erfcg(z,e)
c   an error function complement subroutine.
      implicit real*8(a-h,o-z)
      common/const/ fc,r,t,frt,cur,ep3,ep2,pi,ep1,epf3,epf1,
     &epp1,epp2,epp3,shape3,shape1,capp1,capp3,nneg,nprop,npos
      common/gas/ epg1,epg2,epg3
c
      a1=0.254829592d0
      a2=-0.284496736d0
      a3=1.421413741d0
      a4=-1.453152027d0
      a5=1.061405429d0
      if(z .lt. 2.747192d0) then
      t3=1.0d0/(1.0d0+0.3275911d0*z)
      e=(a1*t3+a2*t3*t3+a3*(t3**3)+a4*(t3**4)
     1+a5*(t3**5))*expg(-z*z)
      else
c
      if(z .gt. 25.0d0) then
      e=0.0d0
      else
c
      sum=0.0d0
      max=z*z + 0.5d0
      fac=-0.5d0/z/z
      sum=fac
      tl=fac
      n=1
   10 n=n+1
      if(n .gt. max) go to 15
      tn=tl*dble(2*n-1)*fac
      sum=sum + tn
      if(tn .lt. 1.0d-06) go to 15
      tl=tn
      go to 10
   15 e=(expg(-z*z))*(1.0d0+sum)/dsqrt(pi)/z
      end if
      end if
c
      return
      end
c***********************************************************************
      subroutine band(j)
c     subroutine to solve coupled, linear differnece equations
      implicit real*8(a-h,o-z)
      common /n/ tmmax,imp,ji,nx,nt,n1,n2,nj,n3,nconv,npa,iSk3,kj3
     &,ii2,ki2,kj,i2div,ip2,kp2,ip1,kp1,imb2,kS2,imb1,kS1,iSk1,kj1
     &,iSk2,kj2
      common/mat/ b,d
      common/bnd/ a,c,g,x,y
      dimension b(17,17),d(17,35)
      dimension a(17,17),c(17,221),g(17),x(17,17),y(17,17)
      dimension e(17,18,221)
      save e, np1
  101 format (15h determ=0 at j=,i4)
      n=nx
      if (j-2)  1,6,8
    1 np1 = n + 1
      do 2 i=1,n
      d(i,2*n+1)= g(i)
      do 2 l=1,n
      lpn= l + n
    2 d(i,lpn)= x(i,l)
      call matinv(n,2*n+1,determ)
      if (determ)  4,3,4
    3 print 101, j
    4 do 5 k=1,n
      e(k,np1,1)= d(k,2*n+1)
      do 5 l=1,n
      e(k,l,1)= - d(k,l)
      lpn= l + n
    5 x(k,l)= - d(k,lpn)
      return
    6 do 7 i=1,n
      do 7 k=1,n
      do 7 l=1,n
    7 d(i,k)= d(i,k) + a(i,l)*x(l,k)
    8 if (j-nj)  11,9,9
    9 do 10 i=1,n
      do 10 l=1,n
      g(i)= g(i) - y(i,l)*e(l,np1,j-2)
      do 10 m=1,n
   10 a(i,l)= a(i,l) + y(i,m)*e(m,l,j-2)
   11 do 12 i=1,n
      d(i,np1)= - g(i)
      do 12 l=1,n
      d(i,np1)= d(i,np1) + a(i,l)*e(l,np1,j-1)
      do 12 k=1,n
   12 b(i,k)= b(i,k) + a(i,l)*e(l,k,j-1)
      call matinv(n,np1,determ)
      if (determ)  14,13,14
   13 print 101,j
   14 do 15 k=1,n
      do 15 m=1,np1
   15 e(k,m,j)= - d(k,m)
      if (j-nj)  20,16,16
   16 do 17 k=1,n
   17 c(k,j)= e(k,np1,j)
      do 18 jj=2,nj
      m= nj - jj + 1
      do 18 k=1,n
      c(k,m)= e(k,np1,m)
      do 18 l=1,n
   18 c(k,m)= c(k,m) + e(k,l,m)*c(l,m+1)
      do 19 l=1,n
      do 19 k=1,n
   19 c(k,1)= c(k,1) + x(k,l)*c(l,3)
   20 return
      end
c***********************************************************************
      subroutine matinv(n,m,determ)
c     matrix inversion program used by BAND(j).
      implicit real*8(a-h,o-z)
      common/mat/ b,d
      dimension b(17,17),d(17,35)
      dimension id(17)
      determ=1.0d0
      do 1 i=1,n
    1 id(i)=0
      do 18 nn=1,n
      bmax=1.1d0
      do 6 i=1,n
      if(id(i).ne.0) go to 6
      bnext=0.0d0
      btry=0.0d0
      do 5 j=1,n
      if(id(j).ne.0) go to 5
      if(dabs(b(i,j)).le.bnext) go to 5
      bnext=dabs(b(i,j))
      if(bnext.le.btry) go to 5
      bnext=btry
      btry=dabs(b(i,j))
      jc=j
    5 continue
      if(bnext.ge.bmax*btry) go to 6
      bmax=bnext/btry
      irow=i
      jcol=jc
    6 continue
      if(id(jc).eq.0) go to 8
      determ=0.0d0
      return
    8 id(jcol)=1
      if(jcol.eq.irow) go to 12
      do 10 j=1,n
      save=b(irow,j)
      b(irow,j)=b(jcol,j)
   10 b(jcol,j)=save
      do 11 k=1,m
      save=d(irow,k)
      d(irow,k)=d(jcol,k)
   11 d(jcol,k)=save
   12 f=1.0d0/b(jcol,jcol)
      do 13 j=1,n
   13 b(jcol,j)=b(jcol,j)*f
      do 14 k=1,m
   14 d(jcol,k)=d(jcol,k)*f
      do 18 i=1,n
      if(i.eq.jcol) go to 18
      f=b(i,jcol)
      do 16 j=1,n
   16 b(i,j)=b(i,j)-f*b(jcol,j)
      do 17 k=1,m
   17 d(i,k)=d(i,k)-f*d(jcol,k)
   18 continue
      return
      end
c***********************************************************************
      subroutine nucamb(kk,il2)
c     subroutine to print detailed profiles as go along
      implicit real*8(a-h,o-z)
      parameter(maxt=900)
      common /n/ tmmax,imp,ji,nx,nt,n1,n2,nj,n3,nconv,npa,iSk3,kj3
     &,ii2,ki2,kj,i2div,ip2,kp2,ip1,kp1,imb2,kS2,imb1,kS1,iSk1,kj1
     &,iSk2,kj2
      common/activ/ EbarD,Ebarkap,Ebarka,Ebarkc,Ebarr1,Ebarr3,
     &Ebars1,Ebars3,Ebarks1a,Ebarks1c,Ebarks2a,Ebarks2c,
     &Ebarks3a,Ebarks3c
      common /calc/ ai(maxt,5),ai2(maxt,5),ts(maxt),h,h1,h2,h3,hcn,
     1hcp,rr,rrmax,cuL,sumpa(5,221),Rad1pa(5),Rad3pa(5),area1pa(5)
     &,area3pa(5),mcL
      common/pindiv/pwfindiv(5,221),pwfindiv_d(5,221),vf1(5),vf3(5),izfl
      common/const/ fc,r,t,frt,cur,ep3,ep2,pi,ep1,epf3,epf1,
     &epp1,epp2,epp3,shape3,shape1,capp1,capp3,nneg,nprop,npos
      common/gas/ epg1,epg2,epg3
      common/var/ xp(16),xx(17,221),xt(16,221,maxt)
      common/cprop/ sig3,area3,rka3,rka3save,ct3,dfs3,Rad3,cap1,cap3,
     1sig1,area1,rka1,rka1save,ct1,dfs1,Rad1,tw,dfs1save,dfs3save
      common/tprop/df(221),cd(221),tm(221),ddo2(221),ddh2(221),
     1ddf(221),dcd(221),dtm(221),dfu(221),d2fu(221),do2(221),dh2(221)
      common /vdc/ aa(1,1),bb(1,1),cc(1,150),dd(1,3),gg(1),xxx(1,1),
     1yy(1,1),hha,hhc,cssold(220,150),css(220,150),utz(5,221),
     &dsold(220,150),ds(220,150),time,nn,nnj,np,mvdc1,mvdc3,lims
      dimension zz(221)
  109 format(f6.1,', ',f11.1,', ',f8.4,', ',g10.3,', ',f8.3,', ',g10.4
     1,', ',g10.4,', ',g10.4,', ',g10.4,', ',g10.4)
  310 format('Distance   C Elec    C Sol Surf  Liq Pot    Solid Pot  ',
     &'Liq Cur   j main     j side 1    j side 2    j side 3')
  311 format('(microns)    (mol/m3)   x or y     (V)         (V) ',
     1'    (A/m2)     (A/m2)       (A/m2)      (A/m2)      (A/m2)')
   44 format(' t = ',1pe18.6,' min')
c
      do 5 i=1,n1+1
      w=dble(i-1)-1.d-8 ! helps rounding when diff'ing files from diff comps
    5 zz(i) = w*h1*1.0d06
      do 71 i=n1+2,n2+n1
      w=dble(i-(n1+1))-1.d-8
   71 zz(i)=zz(n1+1)+w*h2*1.0d06
      do 72 i=n1+n2+1,nj
      w=dble(i-(n1+n2))-1.d-8
   72 zz(i)=zz(n2+n1)+w*h3*1.0d06
c
      write (4,*) ' '
      write (4,310) 
      write (4,311)

      write (4,44) ts(kk)/60.0d0
c
      do 10 j=1,nj,il2
      if (j .le. n1+1) then
      csol=ct1
      else
      csol=ct3
      end if
      if(j.le.n1+1) then
      curden_main=area1*fc*xx(kj,j)
      curden_side=area1*fc*xx(kj3,j)
      else if(j.ge.n1+n2) then
      curden_main=area3*fc*xx(kj,j) !corrected mixed-up area variables 1/15/99
      curden_side=area3*fc*xx(kj3,j)
      else
      curden_main=0.0d0
      curden_side=0.0d0
      endif
      mpa=1 ! only one size for detailed profiles

   10 write (4,109) zz(j),xx(1,j),xx(2+mpa,j)/csol,xx(kp2,j),xx(kp1,j),
     &xx(ki2,j),xx(kj,j)*fc,xx(kj1,j)*fc,xx(kj2,j)*fc,xx(kj3,j)*fc
      if (kk.eq.1) write (4,*) 'that was the zero-time solution!'
      write (4,*) ' '
      write (4,*) ' '

      return
      end
c***********************************************************************
      subroutine peak(n,lim,curr, vcut)
c     subroutine to calculate the peak-power capability at
c     a particular time in the discharge.
      implicit real*8(a-h,o-z)
      parameter(maxt=900)
      common /n/ tmmax,imp,ji,nx,nt,n1,n2,nj,n3,nconv,npa,iSk3,kj3
     &,ii2,ki2,kj,i2div,ip2,kp2,ip1,kp1,imb2,kS2,imb1,kS1,iSk1,kj1
     &,iSk2,kj2
      common/activ/ EbarD,Ebarkap,Ebarka,Ebarkc,Ebarr1,Ebarr3,
     &Ebars1,Ebars3,Ebarks1a,Ebarks1c,Ebarks2a,Ebarks2c,
     &Ebarks3a,Ebarks3c
      common /calc/ ai(maxt,5),ai2(maxt,5),ts(maxt),h,h1,h2,h3,hcn,
     1hcp,rr,rrmax,cuL,sumpa(5,221),Rad1pa(5),Rad3pa(5),area1pa(5)
     &,area3pa(5),mcL
      common/pindiv/pwfindiv(5,221),pwfindiv_d(5,221),vf1(5),vf3(5),izfl
      common/const/ fc,r,t,frt,cur,ep3,ep2,pi,ep1,epf3,epf1,
     &epp1,epp2,epp3,shape3,shape1,capp1,capp3,nneg,nprop,npos
      common/gas/ epg1,epg2,epg3
      common/var/ xp(16),xx(17,221),xt(16,221,maxt)
      common/cprop/ sig3,area3,rka3,rka3save,ct3,dfs3,Rad3,cap1,cap3,
     1sig1,area1,rka1,rka1save,ct1,dfs1,Rad1,tw,dfs1save,dfs3save
  311 format(f8.5,', ',f7.3,',  ',f8.3,', ',f8.5)
c
c     Peak power current ramp section:
      write (2,*) ' '
      write (2,*) '    PEAK POWER '
      write (2,*) ' '
      write (2,*) 'avg. cell pot, current, power, min pot,'
      write (2,*) '  (V)  ,  (A/m2),  (W/m2),     (V)'
c
c     Duration of current pulse is 30 seconds.
c
      nend = 1
      curpmax = 0.0d0
      curtry = 0.0d0
      pwrmax = 0.0d0
      vvpmax = 0.0d0
      curmin=0.d0
      pwrpmax=0.d0
      rrmax=30.0d0
      cur=curr
      curmax=0.d0
  127 kcount=0
      fact=10.0d0
      vfmax=0.d0
      k=nt+2
      do 126 j=1,nj
      do 126 i=1,n
  126 xt(i,j,k)=xt(i,j,k-1)
      ppow=0.0d0
      ii=0
c     Ramp current:
c
  128 continue
      if (ii.gt.100) return
c
      energ=0.0d0
  130 ii=ii+1
      cur=cur+fact  ! start a new current density
         if(curmax.ne.0.d0) then
         cur=0.5d0*(curpmax+curmax)
         if(mod(ii,2).eq.1) cur=0.5d0*(curpmax+curmin)
         if(cur.eq.0.d0) cur=0.5d0*(curmax+curmin)
         if(vfmax.gt.0.d0 .and. vfmax.lt.vcut) then
         curtry=curpmax+(vcut-vvpmax)/(vfmax-vvpmax)*(curmax-curpmax)
     &+0.01d0*(curmax-curpmax)*dble(mod(ii,3)-1)
         if(curtry.lt.curmax .and. curtry.gt.curpmax) cur=curtry
         endif
      if(vfmax.gt.vcut .and. curmin.gt.0.0d0) then
      v2=pwrmax/curmax
      cur2=curmax
      v1=pwrmin/curmin
      cur1=curmin
      vm=pwrpmax/curpmax
      curm=curpmax
      resis=-((pwrpmax-pwrmin)/(curpmax-curmin)
     &-(pwrpmax-pwrmax)/(curpmax-curmax))/(curmin-curmax)
      Uop=(pwrpmax-pwrmin)/(curpmax-curmin)+resis*(curpmax+curmin)
      curtry=Uop/2.d0/resis+0.1d0*(curmax-curmin)*dble(mod(ii,3)-1)
      if(curtry.lt.curmax .and. curtry.gt.curmin) cur=curtry
      endif
         endif
      kkflag=0
      iflag=0
      nflag=0
      k=nt+2
      timpk=0.0d0
      rr=0.0d0
      ts(k)=ts(k-1)
      call comp(n,lim,k,rr,kkflag,nflag,1,jcount)
      call cellpot(k,vv,0,1,lflag)
      vlast=vv
      rr=0.02d0
  129 kkflag=kkflag+1
      k=k+1
      ts(k)=ts(k-1)+rr
      call calca(k)
c
      call comp(n,lim,k,rr,kkflag,nflag,1,jcount)
c
      if (nflag.eq.1.and.kcount.lt.20) then
      write (2,*) 'Peak current decreased',kcount,fact
      if(cur.lt.curpmax) then
      write (2,*) 'Convergence on power failed; already converged at a
     & higher current'
      write (2,*) 'Best results obtained are:'
      write (2,311) pwrpmax/curpmax,curpmax,pwrpmax,vvpmax
      return
      endif
      curmax=cur
      vfmax=0.d0
      go to 128
      endif
      if (kcount.ge.10) return
c
      call cellpot(k,vv,0,1,lflag)
      energ=energ+(vlast+vv)*(ts(k)-ts(k-1))*cur/2.0d0
c
      timpk=timpk+rr
      if (dabs(timpk-30.0d0).gt.0.1d0) then
c
      if (timpk.lt.30.0d0) then
      vlast=vv
c     Increasing time steps:
      if(jcount.lt.12 .and. kkflag.gt.5 .and. (2.0d0*rr
     1+timpk).lt.30.0d0 .and. iflag.eq.0) then
         rr=rr*2.0d0
      print *,'should not be here'
        print *, 'next time step increased to ', rr,'(s)'
      end if
      if(timpk+rr.gt.30.0d0) iflag=1
      if(timpk+rr.gt.30.0d0) rr=30.0d0-timpk
      go to 129
      end if
c
      end if
      ppow=energ/30.0d0
      write (2,311) ppow/cur,cur,ppow,vv
         if(ppow.gt.pwrpmax .and. vv.gt.vcut) then
         if(curpmax.lt.cur) then
            curmin=curpmax
            pwrmin=pwrpmax
            vfmin=vvpmax
         else
            curmax=curpmax
            pwrmax=pwrpmax
            vfmax=vvpmax
         endif
         curpmax=cur
         pwrpmax=ppow
         vvpmax=vv
         endif
         if(vv.lt.vcut) then
         if(curmax.eq.0.d0 .or. cur.lt.curmax) then
            curmax=cur
            pwrmax=ppow
            vfmax=vv
         endif
         else
         if(cur.gt.curmin .and. cur.lt.curpmax) then
            curmin=cur
            pwrmin=ppow
            vfmin=vv
         endif
         if(cur.gt.curpmax) then
            curmax=cur
            pwrmax=ppow
            vfmax=vv
         endif
         endif
      if(curmax.eq.0.0d0) go to 128
      if(curmin. lt. 0.999d0*curmax) go to 128
      write (2,311) pwrpmax/curpmax,curpmax,pwrpmax,vvpmax
      if(nend .eq. 0) return
      nend = 0
      go to 127
c
      end
c***********************************************************************
      subroutine cellpot(kk,v,li,lpow,lflag)
c     subroutine to calculate the cell potential at the end
c     of a time step and to print results.
      implicit real*8(a-h,o-z)
      parameter(maxt=900)
      COMMON /vdc/ aa(1,1),bb(1,1),cc(1,150),dd(1,3),gg(1),xxx(1,1),
     1yy(1,1),hha,hhc,cssold(220,150),css(220,150),utz(5,221),
     &dsold(220,150),ds(220,150),time,nn,nnj,np,mvdc1,mvdc3,lims
      common /n/ tmmax,imp,ji,nx,nt,n1,n2,nj,n3,nconv,npa,iSk3,kj3
     &,ii2,ki2,kj,i2div,ip2,kp2,ip1,kp1,imb2,kS2,imb1,kS1,iSk1,kj1
     &,iSk2,kj2
      common/activ/ EbarD,Ebarkap,Ebarka,Ebarkc,Ebarr1,Ebarr3,
     &Ebars1,Ebars3,Ebarks1a,Ebarks1c,Ebarks2a,Ebarks2c,
     &Ebarks3a,Ebarks3c
      common /calc/ ai(maxt,5),ai2(maxt,5),ts(maxt),h,h1,h2,h3,hcn,
     1hcp,rr,rrmax,cuL,sumpa(5,221),Rad1pa(5),Rad3pa(5),area1pa(5)
     &,area3pa(5),mcL
      common/pindiv/pwfindiv(5,221),pwfindiv_d(5,221),vf1(5),vf3(5),izfl
      common/const/ fc,r,t,frt,cur,ep3,ep2,pi,ep1,epf3,epf1,
     &epp1,epp2,epp3,shape3,shape1,capp1,capp3,nneg,nprop,npos
      common/gas/ epg1,epg2,epg3
      common/power/ ed,Vold,ranodesave,rcathdesave,heat,qlosstot
      common/var/ xp(16),xx(17,221),xt(16,221,maxt)
      common/cprop/ sig3,area3,rka3,rka3save,ct3,dfs3,Rad3,cap1,cap3,
     1sig1,area1,rka1,rka1save,ct1,dfs1,Rad1,tw,dfs1save,dfs3save
      common/temp/ thk,htc,dudt,Cp,dens,tam,g0,qq,qloss,residm,ncell,lht
      common/side/rksc1,c1init,c2init,rksa1,term_s1(221),vol,
     &rksa2,rksa3,rksc3,rksc2,UsO2,UsH2,cn2,term_s2(221),nside
      common/RG/ RG,RGn,RGp,RGext
      save ut1, ut3
  306 format (f8.5,',   ',f6.2,',   ',g12.5,',   ',g12.5,',   '
     &,g12.5,',   ',g12.5)
  307 format(f9.3,', ',f6.4,', ',f6.4,', ',f8.5,', ',f8.5,', ',f8.2,
     &', ', f6.2, ', ', f6.2, ', ', f6.2, ', ', f6.2, ', ', f6.2)
  309 format(f6.3,', ',f8.5,', ',f5.2,', ',f12.3,',',f9.5,', ',f6.2,
     & ',',g10.3, ',', g10.3)

      if (kk.le.1) then
      ut3=utz(1,nj)
      if (n1.ne.0) ut1=utz(1,1)
      endif

c Define the individual pore-wall fluxes.  These are necessary for
c multiple particle sizes in order to treat properly a film resistance,
c and for a rigorous energy balance.

      do j=1,nj
c     This makes the fluxes the same for different particle sizes
c     when a new leg of the driving profile is encountered.
      if (rr.eq.0.0d0) then ! time step check
      do mpa=1,npa
      if (j.le.n1+1) then
      pwfindiv(mpa,j)=xx(kj,j)*area1pa(mpa)/area1
      elseif (j.ge.n1+n2) then
      pwfindiv(mpa,j)=xx(kj,j)*area3pa(mpa)/area3
      else
      pwfindiv(mpa,j)=0.0d0
      pwfindiv_d(mpa,j)=0.0d0
      endif
      enddo !mpa

      else !rr>0

        if (li.eq.1) then
      kadd=0
      do mpa=1,npa
      if (j.le.n1+1) then
      pwfindiv(mpa,j)=(ai2(1,mpa)*(xt(2+mpa,j,kk-1+kadd)-xx(2+mpa,j))
     &/rr-sumpa(mpa,j))*Rad1pa(mpa)
      if (mvdc1.eq.1) pwfindiv(mpa,j)=xx(kj,j)
      pwfindiv_d(mpa,j)=-Rad1pa(mpa)*ai2(1,mpa)/rr
      utz(mpa,j)=utz(mpa,j)-pwfindiv(mpa,j)*rr*area1pa(mpa)
     &/(1.0d0-ep1-epp1-epf1-epg1)/ct1/vf1(mpa)
      elseif (j.ge.n1+n2) then
      pwfindiv(mpa,j)=(ai(1,mpa)*(xt(2+mpa,j,kk-1+kadd)-xx(2+mpa,j))
     &/rr-sumpa(mpa,j))*Rad3pa(mpa)
      if (mvdc3.eq.1) pwfindiv(mpa,j)=xx(kj,j)
      pwfindiv_d(mpa,j)=-Rad3pa(mpa)*ai(1,mpa)/rr
      utz(mpa,j)=utz(mpa,j)-pwfindiv(mpa,j)*rr*area3pa(mpa)/
     &(1.0d0-ep3-epp3-epf3-epg3)/ct3/vf3(mpa)
      else
      pwfindiv(mpa,j)=0.0d0
      pwfindiv_d(mpa,j)=0.0d0
      endif
      enddo ! mpa
        endif ! li
      endif ! time step check
      enddo ! j

c     Material and current balance tests:
c Current balance - integrate the current balance equation
      sum=0.0d0 !mat balance on electrolyte
      sum1=0.0d0 !current neg elec, main
      sum1_s1=0.0d0 !side current neg elec, side reaction 1
      sum3=0.0d0 !main current pos elec
      sum3_s1=0.0d0 !side current pos elec, side reaction 1
      sum1_s2=0.0d0 !side current neg elec, side reaction 2
      sum3_s2=0.0d0 !side current pos elec, side reaction 2
      sum1_s3=0.0d0 !side current neg elec, side reaction 3
      sum3_s3=0.0d0 !side current pos elec, side reaction 3
      sum4=0.0d0 !mat balance oxygen 
      sum5=0.0d0 !mat balance hydrogen

      if (n1 .gt. 2) then
      do 85 j=2,n1
      sum1_s1=sum1_s1+area1*fc*h1*xx(kj1,j)
      sum1_s2=sum1_s2+area1*fc*h1*xx(kj2,j)
      sum1_s3=sum1_s3+area1*h1*xx(kj3,j)
      sum1=sum1+area1*fc*h1*xx(kj,j) !j main
      sum4=sum4+xx(kS1,j)*h1*epg1 !conc O2
      sum5=sum5+xx(kS2,j)*h1*epg1 !conc H2
   85 sum=sum+xx(1,j)*(ep1+epp1)*h1 !conc elec
      endif
      sum=sum+(xx(1,1)+xx(1,n1+1))*(ep1+epp1)*h1/2.0d0
      sum1=sum1+area1*fc*h1*(xx(kj,1)+xx(kj,n1+1))/2.0d0
      sum1_s1=sum1_s1+area1*fc*h1*(xx(kj1,1)+xx(kj1,n1+1))
     & /2.0d0
      sum1_s2=sum1_s2+area1*fc*h1*(xx(kj2,1)+xx(kj2,n1+1))
     & /2.0d0
      sum1_s3=sum1_s3+area1*h1*(xx(kj3,1)+xx(kj3,n1+1))/2.0d0
      sum4=sum4+h1*(xx(kS1,1)+xx(kS1,n1+1))/2.0d0*epg1
      sum5=sum5+h1*(xx(kS2,1)+xx(kS2,n1+1))/2.0d0*epg1

      do 86 j=n1+2,n2+n1-1
      sum4=sum4+h2*xx(kS1,j)*epg2
      sum5=sum5+h2*xx(kS2,j)*epg2
   86 sum=sum+xx(1,j)*(ep2+epp2)*h2
      sum=sum+(xx(1,n1+1)+xx(1,n2+n1))*(ep2+epp2)*h2/2.0d0
      sum4=sum4+h2*(xx(kS1,n1+1)+xx(kS1,n1+n2))/2.0d0*epg2
      sum5=sum5+h2*(xx(kS2,n1+1)+xx(kS2,n1+n2))/2.0d0*epg2

      do 87 j=n2+n1+1,nj-1
      sum3=sum3+area3*fc*xx(kj,j)*h3
      sum3_s1=sum3_s1+area3*fc*xx(kj1,j)*h3
      sum3_s2=sum3_s2+area3*fc*xx(kj2,j)*h3
      sum3_s3=sum3_s3+area3*fc*xx(kj3,j)*h3
      sum4=sum4+h3*xx(kS1,j)*epg3
      sum5=sum5+h3*xx(kS2,j)*epg3
   87 sum=sum+xx(1,j)*(ep3+epp3)*h3
      sum3=sum3+area3*fc*(xx(kj,n1+n2)+xx(kj,nj))
     & *h3/2.0d0
      sum3_s1=sum3_s1+area3*fc*(xx(kj1,n1+n2)+xx(kj1,nj))
     & *h3/2.0d0
      sum3_s2=sum3_s2+area3*fc*(xx(kj2,n1+n2)+xx(kj2,nj))
     & *h3/2.0d0
      sum3_s3=sum3_s3+area3*fc*(xx(kj3,n1+n2)+xx(kj3,nj))
     & *h3/2.0d0
      sum=sum+(xx(1,n1+n2)+xx(1,nj))*h3*(ep3+epp3)/2.0d0
      sum4=sum4+h3*(xx(kS1,n1+n2)+xx(kS1,nj))/2.0d0*epg3
      sum5=sum5+h3*(xx(kS2,n1+n2)+xx(kS2,nj))/2.0d0*epg3
      sum4=sum4/thk !adjusts to average concentration
      sum5=sum5/thk
c     calculate total salt in cell from initial profile:
      w=xt(1,n1+2,1)*(dble(n2-1)*(ep2+epp2)*h2+dble(n1)*(ep1+epp1)*h1
     1+dble(n3)*(ep3+epp3)*h3)
      if(lflag.eq.0) w=w-(xt(1,n1+2,1)-xt(1,1,1))*(dble(n1)*
     1(ep1+epp1)*h1+dble(n3)*(ep3+epp3)*h3)
      if(lflag.eq.0) w=w-(xt(1,n1+2,1)-xt(1,1,1))*(ep2+epp2)*h2
c     material-balance parameter should be ca=1.00
      ca=sum/w

c     Calculate cell potential from dif of solid-phase potentials:
      v=xt(kp1,nj,kk)-xt(kp1,1,kk)-RG*cur !WHT 1-11-08

c Calculate the resistances in the cell, R=V/I
      mpa=1 ! only one size for resistance calculation
      xxsave=xx(2+mpa,n1+1)
      call ekin(n1+1,kk,0,0.0d0,mpa,mvdc1,mvdc3)
      xx(2+mpa,n1+1)=xxsave
      Ua_end=g0
      xxsave=xx(2+mpa,n1+n2)
      call ekin(n1+n2,kk,0,0.0d0,mpa,mvdc1,mvdc3