c *******btwaterlk3da.for  August 3, 2012
c include the basic parameters for the forward scattering phase noise
c modified for double precision on the variables that are summed and differenced frequently
c modified to set dydx =0 when the pressure goes negative
c modified July 27 to keep adsorption time from becoming negative and also to keep
c occupation probability less or equal to 1
c also fix an error in writing the surf vector
c Version uses simple Runge Kutta solution technique without adaptive stepsize
c
c This version writes and reads a data file with the pressure, outgassing rate and occupation probability
c to allow continuation of a previous run or to enable a subsequent run.
c
c Program calculates the long term effects of a water leak in the beamtube
c after it has been baked out. It operates from a time vs temperature schedule which
c includes the activation of different pumps, the opening and closing of  a leak, the effect of
c a secondary bake.
c
c The program uses the surface physics method of the waterbake programs but also allows the
c the tube to have gradients and different pumping arrangements at different times. The leak can
c be located anywhere in the tube and is open and shut at designated times. The pumping arrangement can
c be changed at designated times. The pressure distribution is determined from the outgassing rate in each section
c determined by the surface physics routines and the pressure in each section is determined by a finite difference 
c Runge Kutta solutions of the free molecule diffusion.
c  The surface properties need to be calculated along with the pressure distribution
c The results of the surface calculation are the outgassing rate in each section aj(n)
c
c The calculation is made in mixed cgs units with the pressure in torr, surface area
c in cm^2, volumes in cc, velocities in cm/sec, surface loading particles/cm^2, 
c outgassing in torr cc/cm^2
c
c conversion factor are used from input engineering values to the cgs values
c
c      particle/cm^2 = 1.0e15 * monolayers
c      particles/cm^3 = (296.0/T)*3.0e16*torr
c      v = velocity of gas cm/sec = 5.26e4*sqrt(18.0/amu)*sqrt(T/296)
c
c  symbols
c     a = tube radius = 62.0 cm
c     al0 = tube length = 4.0e5 cm
c     an = number of sections 
c     al = length/section = al0/an
c     np= number pump ports = 17 evenly spaced at 250meters (16gaps)
c     ftrap= pumping speed of traps at ends = 1.0e5 liters/sec
c     fpump= pumping speed of bakeout cryo pumps along tube = 6.25e3 liters/sec
c     fionp = pumping speed of ion pumps along tube = 100 liters/sec
c     amu = atomic mass of the gas units of atomic hydrogen
c
c
c     tdr0 = DR temperature peak = 1.0e4 K
c     sigma0= initial surface water loading = 150 monolayers
c     rp = fractional repulsive potential = 0.7
c     alpha = water accomodation coefficient = 0.5
c     tau0 = molecule oscillation period on surface = 1.0e-12 sec
c     tresmin = minimum residence time on surface=1.0e-6 sec
c     tresmax = maximum residence time on surface=1.0e10 sec
c     apressrec = maximum pressure to begin calculation = 20.0 torr
c     at(k) = the activation temperature of the surface site k
c     ap(n,k) = the probability that site k is occupied in tube segment n
c     temp = the surface temperature
c
c      aj(n) = outgassing rate torr liters/sec in tube section n
c      f(n)  = pumping speed at section n  liters/sec
c     qlk(n) = leak at section n torr liters/sec
c     surf(n) = surface loading in section n monolayers
c     ql = leak size in torr liters/sec
c     nsecwlk = the section number of the leak
c     tlkon = time leak turns on days
c     tlkoff = time leak turns off days
c
c
c
c
c Difference equation to solve
c
c Section 1:  dp/dt = ((2*v*a)/(3*al**2))*(p(2)-p(1)) + (2/a)*aj(1)
c end pt            - (f*(1)*p(1))/(pi*al*a**2) + qlk(1)/(pi*al*a**2)
c                  
c
c
c Section k : dp/dt = ((2*v*a)/(3*al**2))*(p(k-1)-2*p(k) +p(k+1)) 
c                     + (2/a)*aj(k)  
c                    -(f(k)*p(k)/(pi*al*a**2)+ qlk(k)/(pi*al*a**2)
c                     
c
c last section dp/dt = ((2*v*a)/(3*al**2))*(p(nsec-1)-*p(nsec)) 
c                     + (2/a)*aj(nsec)            
c                     -(f(nsec)*p(nsec))/(pi*al*a**2)+ qlk(nsec)/(pi*al*a**2)
c
c
c    Dimensions
c    number of  output points for file = 9000
c    number of sections = 500
c    number of activation site energies = 1024
c
c
c The program uses the Runge-Kutta algorithms given in Press et al 
c Numerical Recipes 2 for the pumping dynamics and a simple integration
c for the probability of occupancy evolution
c 
        use winteracter
        character fileout*150, fileout1*150
        character filein*150
        real*8 ystart,ap,ajj
        dimension p(500,9000),ystart(500),time(9000)
        dimension qlk(500),aj(500,9000),at(1024),ap(500,1024),w(1024)
        dimension f(500),surf(500,9000),isecpump(20),ajj(500)
        dimension pumpvon(20),pumpvec(20,20),wg(500)
        dimension timec(100),tempc(100),hstr(500,9000)
        common /deriv/ aa,bb,cc,dd,alen0,nsec,al
        common /deriv1/ tdr0,sigma0,rp,alpha,tau0,temp,v
        common/deriv2/nsmin,nsmax
c index to indicate if input data is used
          iread = 0
1        write(unit=*,fmt=2)
2        format( ' enter operation: ' /
     & ' 1 = choose program parameters and defaults ' /
     & ' 2 = choose leak and  pumping parameters '/
     & ' 3 = choose time increments for computation and recording ' /
     & ' 4 = choose schedule of temperature and times ' /
     & ' 5 = read/write data file '/
     & ' 6 = review parameters '/
     & ' 7 = run computation ' /
     &'  8 = write output files ' /
     & ' 9 = estimate phase noise' /
     & ' 10 = end program '/
     & ' mode = : '$)
           read(unit=*,fmt=3)mode
3         format(i3)
           goto(100,200,300,400,500,600,700,800,900,1000)mode
c
c choose program parameters
100     write(unit=*,fmt=150)
150     format(' enter start time in days :'$)
        read(unit=*,fmt=104)timein
        write(unit=*,fmt=101)
101     format(' enter # of sections :' $)
            read(unit=*,fmt=102)nsec
102     format(i6)
            an = real(nsec)
            write(unit=*,fmt=103)
103     format(' enter tube length cm : '$)
            read(unit=*,fmt=104)al0
104      format(e15.6)
            if(al0.eq.0.0)al0 = 4.0e5
            al = al0/an	
            write(unit=*,fmt=105)
105      format(' enter tube diameter cm : '$)
            read(unit=*,fmt=104)a
            if(a.eq.0.0) a= 62.0
            write(unit=*,fmt=106)
106      format(' enter amu of gas : '$)
            read(unit=*,fmt=104)amu
            if(amu.eq.0.0)amu = 18.0
            write(unit=*,fmt=107)
107     format(' enter init H20 load (monolayers), R-D pk Tact (K)'/
     &  'fractional well depth : '$)
            read(unit=*,fmt=108)sigma0,tdr0,rp
108     format(3e15.6)
            if(sigma0.eq.0.0) then
            sigma0= 150.0
            tdr0 = 1.0e4
            rp = 0.7
            end if
            write(unit=*,fmt=109)
109     format(' enter accom. coeff.,turn on press.(torr): '$)
            read(unit=*,fmt=110)alpha,apressrec
110     format(2e15.6)
            if(alpha.eq.0.0.and.apressrec.eq.0.0)then
            alpha = 0.5
            apressrec=20.0
            end if
            write(unit=*,fmt=111)
111      format(' enter min and max res.times (sec) : '$)
            read(unit=*,fmt=112)tresmin,tresmax
112      format(2e15.6)
             if(tresmin.eq.0.0.and.tresmax.eq.0.0)then
             tresmin=1.0e-6
             tresmax=1.0e20
             end if
             write(unit=*,fmt=113)
113       format(' enter molecule oscillation time sec : '$)
             read(unit=*,fmt=104)tau0
             if(tau0.eq.0.0)tau0=1.0e-12
             emtmin=alog(tresmin/tau0)
             emtmax=alog(tresmax/tau0)	
c thermal velocity need to modify for temperature
             vt = (5.263e4)*sqrt(18.0/amu)
               bb = 2.0/a
               dd = 1.0/(3.14159*al*a**2)
c optical parameters of the arm cavity
             write(unit=*,fmt=119)
119       format(' enter w0 for cavity cm : '$)
             read(unit=*,fmt=104)wgg0
             if(wgg0.eq.0.0)wgg0=1.22
             write(unit=*,fmt=120)
120       format(' enter Rayleigh length cm : '$)
             read(unit=*,fmt=104)rayl
             if(rayl.eq.0.0)rayl=4.27806e4
             write(unit=*,fmt=121)
121       format(' pos of waist from lvea cm: '$)
             read(unit=*,fmt=104)zw
             if(zw.eq.0.0)zw=1.83422e5
             write(unit=*,fmt=122)
122       format(' molecular polarizability at 1 micron cm**3/mol : '$)
             read(unit=*,fmt=104)alpol
             if(alpol.eq.0.0)alpol=1.32e-24
c value for water, x(f) and h(f) for one arm (assume other is perfect)
c  x**2(f)= (8*pi**2*alpol**2*(3e+16)/v)*sum(p(torr)*seclen(cm)/w(cm))
c           phcoeff=sqrt((8*pi**2*alpol**2*3e+16)/v)
c            <h(f)> = sqrt(x**2(f))/al0(cm)
              do 125 n=1,nsec
              za = (real(n-1)*al)-zw
              wg(n) = wgg0*sqrt(1.0+(za/rayl)**2)
125       continue            
              go to 1


c leak and pumping parameters

200     write(unit=*,fmt=201)
201     format('  enter leak size t l/s and location m from lvea : '$)
            read(unit=*,fmt=202)ql,xlk
202     format(2e15.6)
c convert torr liters/sec to torr cc/sec
            ql=1000.0*ql
c establish the segment with the leak
            nsecwlk =int(an*100.0*xlk/al0)
            write(unit=*,fmt=203)
203     format(' enter time start and end lk days : '$)
            read(unit=*,fmt=202)tlkon,tlkoff
c       establish the segments with pump ports
            isecpump(1) = 1
            do 205 k=2,17
            ak= real(k)
            isecpump(k) = int(250.0*(ak-1.0)*an/4000.0)
205     continue
            ip=0
210     ip=ip+1
            write(unit=*,fmt=211)ip
211      format(' enter start (days) pump vect '  i3,' : ' $)
            read(unit=*,fmt=212)pumpvon(ip)
212     format(e15.6)
            do 206 k=1,17
            write(unit=*,fmt=213)k
213     format(' enter pumping speed l/s, pump at loc ' i4, ' : ' $)
            read(unit=*,fmt=214)speed
c  convert to cc/sec
            pumpvec(ip,k)=1000.0*speed
206     continue
214     format(e15.6)
            write(unit=*,fmt=215)
215     format(' enter 1 for another pumpvec, cr> to go on : ' $)
            read(unit=*,fmt=216)ig
216     format(i3)
            if(ig.eq.1)go to 210
            ipmax = ip
c set all leaks and pumping speeds to zero
            do 220 k=1,nsec
            qlk(k) = 0.0
            f(k) = 0.0
220     continue
c setup surface physics
            nsite = 1024
c initial occupation probability for each site = 1.0
            do 225 n=1,nsec
            ajj(n)=0.0 ! first time only NOTE Inconsistent with later program
            do 226 k=1,1024
            ap(n,k)=1.0
226     continue
225     continue
c assign activation temperatures 0 to 3*tdr0
            deltat=3.0*tdr0/1024.0
            do 227 k=1,1024
            at(k)=real(k)*deltat
227     continue
c establish weighting function w(k)
           sum =0.0
           do 228 k=1,1024
           w(k)=(2.0*at(k)*deltat/tdr0**2)*exp(-(at(k)/tdr0)**2)
           sum = sum +w(k)
228    continue
           do 229 k =1,1024
c the normalized weighting function
           w(k)=w(k)/sum
229    continue 
            go to 1

c set program times and scaling
300     write(unit=*,fmt=301)
301     format(' enter record step time days, #comp steps/rec : '$)
            read(unit=*,fmt=302)tstep1,intstep
302     format(e15.6,i8)
c convert computation step size to seconds
            tstep = tstep1*3600.0*24.0/real(intstep)
            write(unit=*,fmt=310)
310     format(' enter RK step size sec commensurate wcomp step : '$)
            read(unit=*,fmt=311)rkst
311     format(e15.6)
             go to 1

c make schedule times and temperatures
c
400      do 401 k=1,100
             tempc(k)=0.0
             timec(k)=0.0
401      continue
             write(unit=*,fmt=402)
402      format(' enter filename for read/write of time vs temp : ' $)
             read(unit=*,fmt=403)fileout1
403      format(a150)
             write(unit=*,fmt=404)
404      format(' enter 1 if time and temp to be read from file : '$)
             read(unit=*,fmt=405)ittread
405      format(i3)
             if(ittread.eq.0) go to 411
             open(unit=3,file=fileout1)
             read(unit=3,fmt=406)maxstep
406      format(i5)
             write(unit=*,fmt=407)maxstep
407       format(' #steps = 'i4)
              do 408 k=1,maxstep
              read(unit=3,fmt=409)tempc(k),timec(k)
409       format(2e15.6)
              write(unit=*,fmt=410)k,tempc(k),timec(k)
410        format(' step # = 'i4,' temp = 'e15.6,' time = 'e15.6)
408        continue
              close(3)
              go to 1
411       write(unit=*,fmt=412)
412     format( 'enter the schedule of temperature and times,'/
     &  ' first point must be initial temperature and zero time')
            jc=0
            do 420 k=1,100
            jc = jc + 1
            write(unit=*,fmt=413)
413     format(' temp (K) = , time (days) = : '$)
            read(unit=*,fmt=414)tempc(jc),timec(jc)
414     format(e15.6,e15.6)
            if(tempc(jc).eq.0.0)then
            maxstep=jc-1
            open(unit=3,file=fileout1)
            write(unit=3,fmt=415)maxstep
415     format(i5)
           do 416 kl=1,maxstep
           write(unit=3,fmt=414)tempc(kl),timec(kl)
           write(unit=*,fmt=421)kl, tempc(kl),timec(kl)
416     continue
            close(3)
            go to 1
            end if
420      continue
421      format(i5,' temp= '1pe15.6,' time=' 1pe15.6)
            go to 1

c read and write data files
500     write(unit=*,fmt=501)
501     format(' read data file = 1, write data file = 2 :'$)
            read(unit=*,fmt=502)iread
502     format(i5)
            if(iread.eq.2)then
            write(unit=*,fmt=503)
503     format(' enter name of data file : '$)
            read(unit=*,fmt=504)fileout
504     format(a150)
            open(unit=3,file=fileout)
            write(unit=3,fmt=505)nsec
505     format(i5)
           do 506 n=1,nsec
           do 507 k=1,1024
           write(unit=3,fmt=508)ap(n,k)
507    continue
506    continue
508    format(1pe15.6)
          do 509 n=1,nsec
          write(unit=3,fmt=510)p(n,its),aj(n,its)
509    continue
510    format(1pe15.6,1pe15.6)
          close(3)
          end if
          if(iread.eq.1)then
          write(unit=*,fmt=511)
511     format('enter name of data file : '$)
           read(unit=*,fmt=512)filein
512    format(a150)
          open(unit=3,file=filein)
          read(unit=3,fmt=505)nsec
          do 513 n=1,nsec
          do 514 k=1,1024
          read(unit=3,fmt=508)ap(n,k)
514     continue
513     continue
          do 515 n=1,nsec
          read(unit=3,fmt=510)ystart(n),ajj(n)
          ajj(n)=ajj(n)*1000.0
515    continue
         close (3)
         end if
         go to 1




c review parameters
600    write(unit=*,fmt=650)timein
650    format('  timein(days) = '1pe15.6)
       write(unit=*,fmt=601)nsec
601    format(' nsec =' i5)
           write(unit=*,fmt=602)al0
602    format(' al0 =' 1pe15.6)
           write(unit=*,fmt=603)a
603    format(' a =' 1pe15.6)
           write(unit=*,fmt=604)amu
604    format(' amu =' 1pe15.6)
           write(unit=*,fmt=605)sigma0,tdr0,rp
605    format(' sigma0 =' 1pe15.6,' tdr0= '1pe15.6,' rp = '1pe15.6)
           write(unit=*,fmt=606)alpha,apressrec
606    format(' alpha =' 1pe15.6,' apressrec = '1pe15.6)
           write(unit=*,fmt=607)tresmin,tresmax,tau0
607    format(' tresmin = ' 1pe15.6,' tresmax = '1pe15.6, ' tau = '
     &  1pe15.6)
           write(unit=*,fmt=608)emtmin,emtmax
608    format(' emtmin = ' 1pe15.6, ' emtmax = ' 1pe15.6)
           write(unit=*,fmt=609)ql,xlk,nsecwlk
609    format(' ql= ' 1pe15.6, ' xlk = '1pe15.6, ' nsecwlk =' i6)
           write(unit=*,fmt=610)tlkon,tlkoff
610    format(' tlkon = ' 1pe15.6, ' tlkoff = '1pe15.6)
           do 612 k=1,17
           write(unit=*,fmt=611)k,isecpump(k)
611    format(' k = ' i6,' isecpump = 'i6)
612    continue
          do 613 n= 1,ipmax
           write(unit=*,fmt=614)n,pumpvon(n)
614    format('n = ' i5, ' pumpvon =' 1pe15.6)
          do 615 k = 1,17
           write(unit=*,fmt=616)n,k,pumpvec(n,k)
616    format(' n = 'i6,' k = 'i6, ' pumpvec = '1pe15.6)
615    continue
613    continue
          do 617 n=1,nsec
          write(unit=*,fmt=618)n,wg(n)
618    format(' n = 'i5, 'spot radius cm = '1pe15.6)
617    continue          
           write(unit=*,fmt=621)rkst
621    format(' Runge_Kutta stepsize (sec) = '1pe15.6)
            go to 1

c run computation time in hours sets all events
700       itrec=1
              time(1)=timein
              its=1
              kk=1  
             x2 = timein*3600.0*24.0
c begin loop set by schedule of events
701     if(itrec.eq.maxstep)then
            write(unit=*,fmt=730)
730     format(' ended by last entry of times and temperatures' )
            go to 1
            end if
       if(time(its).ge.timec(itrec).and.time(its).lt.
     & timec(itrec+1))then
            temp=tempc(itrec)
            end if
        if(time(its).ge.timec(itrec+1))then
           itrec=itrec+1
           temp=tempc(itrec)
           end if
             v = vt*sqrt(temp/296.0)
           phcoeff=(3.14159*alpol*sqrt(8.0*3.0e16/v))/al0
            aa = (2.0*v*a)/(3.0*al**2)
c set the range of activation sites used in the calculation
            nsmin=int((1024.0*temp*emtmin)/(3.0*tdr0))
            if(nsmin.lt.1) nsmin=1
            nsmax=int((1024.0*temp*emtmax)/(3.0*tdr0))
            if(nsmax.gt.1024) nsmax=1024
            do 7025 n=1,nsec
            do 7024 k=1,nsmin
            ap(n,k)=0.0
7024   continue
            do 7026 k=nsmax,1024
            ap(n,k)=1.0
7026   continue
7025   continue
c begin loop set by sampling time
c check for pumping arrangement
7045            if(kk.eq.ipmax)go to 6010
                     if(kk.lt.ipmax)then
       if(time(its).ge.pumpvon(kk).and.time(its).lt.pumpvon(kk+1))then
              do 7022 k=1,17
              n = isecpump(k)
              f(n) = pumpvec(kk,k)
7022     continue
             go to 6020
              end if
              if(time(its).ge.pumpvon(kk+1))then
              kk=kk+1
              go to 7045
              end if
             end if
6010     do 7029 k=1,17
              n = isecpump(k)
              f(n) = pumpvec(kk,k)
7029     continue
c establish if there is a leak 
6020         qlk(nsecwlk)=0.0
        if(time(its).ge.tlkon.and.time(its).lt.tlkoff)qlk(nsecwlk)=ql
c        write(unit=*,fmt=4002)time(its),tlkon,tlkoff,nsecwlk,
c     &   qlk(nsecwlk)
c4002   format('time='1pe12.3,'on='1pe12.3,'off='1pe12.3,'nlk='i5,
c     &   'qlk ='1pe12.3)
c begin loop of the internal calculation
703            do 799 j=1,intstep
                   x1 = x2
                   x2 = x1 + tstep
                   if(j.eq.1.and.its.eq.1.and.iread.eq.0)then
                   do 704 n=1,nsec
                   ystart(n)= apressrec
                   p(n,1)=apressrec
704             continue
                   p(nsec+1,1)=apressrec
                   end if
                  if(its.eq.1.and.j.eq.1.and.iread.eq.1)go to 9999
c establish the outgassing rate for each section
        call probev(ystart,tstep,ap,at,ajj,w)
c produces ajj(n) for this tstep, apply to pressure calc
9999    call rkdumb(ystart,x1,x2,dydx,rkst,qlk,f,ajj)
 799           continue
                its=its +1
                time(its) = x2/(3600.0*24.0)
                sump =0.0
                sumj=0.0
                sums=0.0
                sumh=0.0
                do 705 n=1,nsec
                p(n,its)=ystart(n)
                sump = sump+ystart(n)/an
c outgassing rate written to vector as torr liter/sec
                aj(n,its)=ajj(n)/1000.0
                sumj=sumj + aj(n,its)/an
               sumtr=0.0
                do 740 k=1,1024
                 sumtr = sumtr + sigma0*w(k)*ap(n,k)
740        continue
                surf(n,its)=sumtr
                sums=sums + sumtr/an
c need to avoid negative pressures
                if(p(n,its).lt.0.0)p(n,its)=0.0
                hstr(n,its) =  p(n,its)*al/wg(n)
c the running hstr(n,its) values are the pressure*section l/gaussian spot
                sumh = sumh +hstr(n,its)
705         continue
                p(nsec+1,its)=sump
               aj(nsec+1,its) =sumj
                surf(nsec+1,its) = sums
                hstr(nsec+1,its)=phcoeff*sqrt(sumh)
c       the hstr(nsec+1,its) value is the <h(f)> along the tube
         write(unit=*,fmt=706)time(its), p(nsec+1,its),
     & aj(nsec+1,its),surf(nsec+1,its),hstr(nsec+1,its)
706    format(' time= '1pe12.4, ' <p>= '1pe12.4,' <j>= '1pe12.4,
     &  ' <surf>= '1pe12.4, ' h(f) = ' 1pe12.4)
                go to 701

c write output files
800    write(unit=*,fmt=801)
801    format(' enter output file : ' $)
           read(unit=*,fmt=802)fileout
802    format(a150)
           open(unit=3,file=fileout)
           write(unit=*,fmt=803)
803    format(' enter 1= p,2=j ,3=sigma, 4 = h(f) : '$)
           read(unit=*,fmt=804)igo
           write(unit=*,fmt=812)
812     format(' enter 1 for param vs t, 2 param vs n : '$)
           read(unit=*,fmt=804)igg
           if(igg.eq.1)then
           write(unit=*,fmt=806)
806    format(' enter section # nsec+1 = <> ' :$)
           read(unit=*,fmt=804)n
           write(unit=3,fmt=804)its-1
804     format(i5)
          do 805 k=2,its
          if(igo.eq.1)write(unit=3,fmt=810)time(k),p(n,k)
          if(igo.eq.2)write(unit=3,fmt=810)time(k),aj(n,k)
          if(igo.eq.3)write(unit=3,fmt=810)time(k),surf(n,k)
          if(igo.eq.4)write(unit=3,fmt=810)time(k),hstr(n,k)
805   continue
          close (3)
          end if
          if(igg.eq.2)then
          write(unit=3,fmt=804)nsec
          do 820 nn=1,nsec
          bn = al*real(nn)/100.0
          if(igo.eq.1)write(unit=3,fmt=810)bn,p(nn,its)
          if(igo.eq.2)write(unit=3,fmt=810)bn,aj(nn,its)
          if(igo.eq.3)write(unit=3,fmt=810)bn,surf(nn,its)
          if(igo.eq.4)write(unit=3,fmt=810)bn,hstr(nn,its)
820   continue
          close (3)
          end if
          write(unit=*,fmt=807)
807    format(' enter 1 for another file, cr> to go on :' $)
          read(unit=*,fmt=804)ig
          if(ig.eq.1) go to 800
810   format(1pe15.6,1pe15.6)
         go to 1


c estimate phase noise
900         continue
                go to 1 

1000    continue
	end

        subroutine derivs(x,y,dydx,f,ajj,qlk)
        real*8 y,dydx,ajj
        dimension y(500),dydx(500)
        dimension f(500),ajj(500),qlk(500)
        common /deriv/ aa,bb,cc,dd,alen0,nsec,al
        do 100 n=2,nsec-1
        dydx(n) = aa*(y(n-1)+y(n+1)-2.0*y(n))
     & +bb*ajj(n)+dd*(qlk(n)-f(n)*y(n))
        if(y(n-1).lt.0.0.or.y(n).lt.0.0.or.y(n+1).lt.0.0)dydx(n)=0.0
100     continue
        dydx(1) = aa*(y(2)-y(1)) + bb*ajj(1)
     &  +dd*(qlk(1)-f(1)*y(1))
        if(y(1).lt.0.0.or.y(2).lt.0.0)dydx(1)=0.0 
        dydx(nsec)=aa*(y(nsec-1)-y(nsec))+bb*ajj(nsec)
     &  +dd*(qlk(nsec)-f(nsec)*y(nsec))
        if(y(nsec-1).lt.0.0.or.y(nsec).lt.0.0)dydx(nsec)=0.0
          return
          end

       subroutine probev(ystart,tstep,ap,at,ajj,w)
       real*8 ap,ystart,ajj,suma,sumb,rt,tg,tt,temit,tads
       real*8 tau,pequil,ax
       dimension at(1024),ap(500,1024),w(1024),ajj(500)
       dimension ystart(500)
       common /deriv/ aa,bb,cc,dd,alen0,nsec,al
       common /deriv1/ tdr0,sigma0,rp,alpha,tau0,temp,v
       common/deriv2/nsmin,nsmax
c loop over the sections
       do 100 n=1,nsec
c loop over the useful binding sites
        sumb=0.0
        suma=0.0                 
        do 200 k=nsmin,nsmax
        sumb = sumb + ap(n,k)*w(k)
c estimate the the changes in ap
        rt=at(k)/temp
        tg=(1.0-rp)*rt
        tt=(1.0+tg)*dexp(-tg)
        temit=tau0*dexp(rt)
c the 30 is the #/cm**2 per monolayer/#cc/torr
        tads=4.0*sigma0/(30.0*alpha*ystart(n)*v*tt)
        tau=(temit*tads)/(temit+tads)
        pequil = temit/(temit+tads)
        if(tads.gt.1.0e12.or.tads.le.0.0)then
        tau= temit
        pequil = 0.0
        end if
        ax = dexp(-tstep/tau)
        ap(n,k)=ap(n,k)*ax + pequil*(1.0-ax)
        if(ap(n,k).gt.1.0)ap(n,k)=1.0
        if(ap(n,k).lt.0.0)ap(n,k)=0.0
        suma=suma+ap(n,k)*w(k)
200    continue
c 30 converts monolayers/sec to torr cc/sec/cm^2
        ajj(n)=(sumb-suma)*sigma0/(30.0*tstep) 
100    continue
       return
       end

      SUBROUTINE rk4(v,dydx,x,h,ystart,f,ajj,qlk)
      INTEGER nsec,NMAX
      REAL*8 dydx,yout,ajj,ystart,y,v,dv
      real*4 x,xx,h6,hh,xh
        dimension f(500),ajj(500),qlk(500),ystart(500),dydx(500)
        dimension y(500,500),v(500),xx(500),dv(500)
      PARAMETER (NMAX=500)
      INTEGER i
      REAL*8 dym(NMAX),dyt(NMAX),yt(NMAX)
        common /deriv/ aa,bb,cc,dd,alen0,nsec,al
      hh=h*0.5
      h6=h/6.
      xh=x+hh
      do 11 i=1,nsec
        yt(i)=v(i)+hh*dydx(i)
11    continue
      call derivs(xh,yt,dyt,f,ajj,qlk)
      do 12 i=1,nsec
        yt(i)=v(i)+hh*dyt(i)
12    continue
      call derivs(xh,yt,dym,f,ajj,qlk)
      do 13 i=1,nsec
        yt(i)=v(i)+h*dym(i)
        dym(i)=dyt(i)+dym(i)
13    continue
      call derivs(x+h,yt,dyt,f,ajj,qlk)
      do 14 i=1,nsec
        ystart(i)=v(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
14    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 7%W3.


       SUBROUTINE rkdumb(ystart,x1,x2,dydx,rkst,qlk,f,ajj)
       INTEGER nstep,nsec,NMAX,NSTPMX
       PARAMETER (NMAX=500,NSTPMX=500)
       REAL*4 x1,x2,xx,x,h
        real*8 ystart,y,v,ajj,dydx,dv
        dimension f(500),ajj(500),qlk(500),ystart(500),dydx(500)
        dimension y(500,500),v(500),xx(500),dv(500)
CU    USES rk4
       INTEGER i,k
        common /deriv/ aa,bb,cc,dd,alen0,nsec,al
        nstep = int((x2-x1)/rkst)
        do 11 i=1,nsec
        v(i)=ystart(i)
        y(i,1)=v(i)
11    continue
      xx(1)=x1
      x=x1
      h=(x2-x1)/real(nstep)
      do 13 k=1,nstep
        call derivs(x,v,dv,f,ajj,qlk)
        call rk4(v,dv,x,h,ystart,f,ajj,qlk)
        if(x+h.eq.x)pause 'stepsize not significant in rkdumb'
        x=x+h
        xx(k+1)=x
        do 12 i=1,nsec
          y(i,k+1)=v(i)
12      continue
13    continue
         return
        END
C  (C) Copr. 1986-92 Numerical Recipes Software 7%W3.