c234567 Runge-Kutta-Gill program for quasi-reversible CV of adsorbed species
      implicit real*8 (a-h,o-z)
      real(8), external :: func
      open (1,file='adsorp_wave.dat')
      read (1,*) a,gamma_tot !electrode area m^-2, max amount of adsorp.
mol m-2
      read (1,*) z,en ! z:valnce of ox, n:electron transfer number
      read (1,*) ekt0 ! rate constant k_t^0
      read (1,*) alpha ! trans. coeffi.
      read (1,*) e0 ! formal potential E^0' / V
      read (1,*) ei,v,th ! CV parameter E_i / V, v / Vs-1, th: return time
      read (1,*) phiad ! EDL phi at ad. /V

      write (6,*) ekt0,v
      fa=96485.3399d0
      r=8.314472d0
      temp=300.0d0

      theta=0.9999908750d0 ! start theta < 1
      t=0.0d0
      tf=2.0d0*th
      mesh=10000000 ! for convergence RKG huge number is required.
      meshw=1000
      delt=tf/real(mesh)
      dtheta_dt=func(t,theta,ekt0,alpha,en,phiad,ei,v,th,e0)

      dtheta_dt=func(t,theta,ekt0,alpha,en,phiad,ei,v,th,e0)
      cur=en*fa*a*gamma_tot*dtheta_dt

      if (t <= th) then
       ee=ei-v*t
      else
       ee=ei-2.0d0*v*th+v*t
      endif
      write (6,*) t,ee,cur,theta

c --RKG routine
      do i=1, mesh
c      write (6,*) t,theta,ekt0,alpha,en,phiad,ei,v,th,e0
       ek1=delt*func(t,theta,ekt0,alpha,en,phiad,ei,v,th,e0)
       t2=t+delt/2.0d0
       theta2=theta+ek1/2.0d0

c      write (6,*) t2,theta2,ekt0,alpha,en,phiad,ei,v,th,e0
       ek2=delt*func(t2,theta2,ekt0,alpha,en,phiad,ei,v,th,e0)
       t3=t2
       theta3=theta+ek1*(-0.5d0+1.0d0/sqrt(2.0d0))
     & +(1.0d0-1.0d0/sqrt(2.0d0))*ek2

       ek3=delt*func(t3,theta3,ekt0,alpha,en,phiad,ei,v,th,e0)
       t4=t+delt
       theta4=theta-ek2/sqrt(2.0d0)+ek3*(1.0d0+1.0d0/sqrt(2.0d0))

       ek4=delt*func(t4,theta4,ekt0,alpha,en,phiad,ei,v,th,e0)

c      write (6,*) ek1,ek2,ek3,ek4
       theta_up=theta+(ek1+2.0d0*(1.0d0-1.0d0/sqrt(2.0d0))*ek2
     & +2.0d0*(1.0d0+1.0d0/sqrt(2.0d0))*ek3+ek4)/6.0d0

       t=t4
       theta=theta_up
       dtheta_dt=func(t,theta,ekt0,alpha,en,phiad,ei,v,th,e0)
       cur=en*fa*a*gamma_tot*dtheta_dt

       if (t <= th) then
        ee=ei-v*t
       else
        ee=ei-2.0d0*v*th+v*t
       endif

       if (mod(i,meshw) == 0 ) then
        write (6,*) t,ee,cur,theta
       endif
      enddo
      stop
      end

      function func(t,theta,ekt0,alpha,en,phiad,ei,v,th,e0) result(ff)

      real(8) :: ff,fa,r,temp,ee
      real(8) :: t,theta,ekt0,alpha,en,phiad,ei,v,th,e0

      fa=96485.3399d0
      r=8.314472d0
      temp=300.0d0

      if (t <= th) then
       ee=ei-v*t
      else
       ee=ei-2.0d0*v*th+v*t
      endif

      ff=-ekt0*exp(alpha*en*fa*phiad/r/temp)
     & *exp(-alpha*en*fa*(ee-e0)/r/temp)
     & *theta + ekt0*exp(-(1.0d0-alpha)*en*fa*phiad/r/temp)
     & *exp((1.0d0-alpha)*en*fa*(ee-e0)/r/temp)*(1.0d0-theta)

      end
c234567