c234567
c --------------- example: adsorp_wave.dat -----
c 1.26d-5 7.67d-6 ! A 4mm-phi,Gamma_tot !electrode area m^-2, max amount of ad. mol m-2 0.74micCm-2=7.67d-6 molm-2 0.25Cm-2=2.6d-6 uchida exp
c 1.0d0 1.0d0 ! z:valnce of ox, n:electron transfer number
c 1100.0d0 ! ekt0:rate constant k_t^0 / s-1 1100 reversible
c 0.5d0 ! alpha:trans. coeffi.
c 0.00d0 ! e0: formal potential E^0' / V
c 0.360d0 0.05d0 12.0d0 ! ei,v,th ! CV parameter E_i / V, v / Vs-1, th: return time
c -0.0d0 ! EDL phi_ad / V
c 6.0d0 -771.88d0 0.0d0 0.0d0    ! znn, eoo,
c---------------------------------------------------------

c234567 reversible surface wave program with interaction
      implicit real*8 (a-h,o-z)i
      dimension ttt(20000),eee(20000),thth(20000),cur(20000)
      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,*) 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
      read (1,*) znn, eoo, err, eor !# of nearest, interaction energy J mol-1
      write (6,*) v
      fa=96485.3399d0   ! faraday constant
      r=8.314472d0   ! gas constant
      temp=300.0d0   ! absolute temperature
      theta0=1.0d0   ! initial value of theta
      t=0.0d0   ! stat time
      tf=2.0d0*th  ! total scan time
      mesh=10000 ! for huge number is required.
      delt=tf/real(mesh)   ! time for one mesh

      do it=1, mesh
       t=t+delt
       if (t <= th) then
	ee=ei-v*t
       else
	ee=ei-2.0d0*v*th+v*t
       endif
c      write (6,*) t,ee,cur,theta
       pjud0=ee-e0-phiad-2.0d0*znn/en/fa*(eor-err)
     &  -r*temp/en/fa*log(theta0/(1.0d0-theta0))
     & -2.0d0*znn/en/fa*(eoo+err-2.0d0*eor)*theta0

c ----- theta search routine ----
       do i=1, 100001i
        theta=theta0-real(i-1)/100001.0d0
        pjud=ee-e0-phiad-2.0d0*znn/en/fa*(eor-err)
     &  -r*temp/en/fa*log(theta/(1.0d0-theta))
     &  -2.0d0*znn/en/fa*(eoo+err-2.0d0*eor)*theta
        pppp=pjud0*pjud
        if (pppp <= 0.0d0) goto 100
        pjud0=pjud
       enddo
c100   write (6,*) t,ee,theta
100    ttt(it)=t
       eee(it)=ee
       thth(it)=theta
      enddo

      do it=1, mesh-1
       cur(it)=en*fa*a*gamma_tot*(thth(it+1)-thth(it))/
     & (ttt(it+1)-ttt(it))
       write (6,*) ttt(it),eee(it),thth(it),cur(it)
      enddo
      stop
      end
c234567