c ---------------------------------------------------------------------------
c
c     program scimms       08/29/1993   Simulation of size-exclusion
c                          chromatography properties of molecules capable of 
c                             homologous (self) association, or
c                             heterologous association       or
c                             non-covalent polymerization.
c 
c                          Heterologous association may be 
c                             mono- or bi-valent.
c          
c                          Heterologous association may be competitive.
c                         
c                          The column may be saturated with one component.
c
c                          Sample size may be set at values that correspond
c                             to small zone or large zone experimental
c                             protocols.
c
c                          Simulated elution profiles may be output
c                             as a function of time, volume, or 
c                             relative volume.
c
c
c                 SCIMMS :  Simulated
c                           Chromatography of
c                           Interactive
c                           MacroMolecular
c                           Systems
c
c
c-----------------------------------------------------------------


        character*4 fil1
        character*5 position_label
        character*12 recfl
        character*20 infile_pro,infile_col
        character*25 infiles(100), infile

        integer ss,cl,mqt,ndif,prfr,prfrs,
     x   prn, rl,prl, sprd, saturate, valence,cell2

        real gs(705),mq(0:200,0:15),eluted(0:15),label(10),
     x   prload,kf(0:3),kr(0:3),mw(0:15),dif(0:15),
     x   pc(0:15),vl(0:15),dis(0:15),ec(0:15)

        double precision sum,p1(-10:6000,0:15),p2(0:6000,0:15),
     x    prel(0:15), t_inc

        data infile_pro /'proteins_inpt.lst'/
        data infile_col /'column_inpt.dat'/
        data fil1 /'ksm_'/

c---------------------------------------------------------------
c       Set Up:
c         load Gaussian distribution table   (GAUSS)
c         input column parameters            (COLUMN)
c         input protein input file list      (PRO_FILE)
c         input protein parameters           (PROTEINS)       
c         calculate weighted extinctions     (EXTINCTION)
c         create unique output file index    (FILE_NUM)
c         document run                       (LOG_DAT)
c         saturate column if chosen          (SATURAT)
c         apply sample to column             (LOADER)
c         pre-run equilibration              (ITER_KIN)
c         setup full output                  (TERM_WRITE)
c         setup plot output                  (PLOT_WRITE)
c----------------------------------------------------------------


         call GAUSS(gs)
         call COLUMN(cl,cv,mqt,ss,rl,prl,t_inc,t_inc_set,delta_max,
     x               t_tot,rtime,ndif,sprd,flo_rt,saturate,
     x               position_label, infile_col)

         call PRO_FILE(infile_pro, infiles,max_files)

         
         do n_run = 1, max_files

              infile = infiles(n_run)
              t_inc = t_inc_set

         call PROTEINS(pc,mw,vl,ec,dif,dis,kf,kr,prn,infile,
     x                 valence,mess_t_inc)
         call EXTINCTION(mw,ec,valence)
         call FILE_NUM(n_rec,fil1,label)
         call SATURAT(pc,mw,vl,p1,p2,saturate,prn,sprd,cl,ss,
     x                dis,gs,mqt,mq,ndif,dif)
         call LOADER(ss,cl,mq,mqt,pc,saturate,mw,p1,valence)
         call TERM_WRITE(cl,rl,mqt,ss,sprd,ndif,vl,dis,t_inc,
     x        delta_max,t_tot,pc,mw,ec,dif,kf,kr,
     x        n_rec,recfl,fil1,06,valence)
         call TERM_WRITE(cl,rl,mqt,ss,sprd,ndif,vl,dis,t_inc,
     x        delta_max,t_tot,pc,mw,ec,dif,kf,kr,
     x        n_rec,recfl,fil1,07,valence)
         call PLOT_WRITE(n_rec,fil1,08)
         call ITER_KIN(prl,p1,p2,pc,mw,kf,kr,delta_max,t_inc,
     x        t_tot,1,ss,valence,mess_t_inc)
         call LOG_DAT(infile,sprd,prn,mqt,cl,cv,rl,ss,ndif,vl,dis,pc,
     x                mw,ec,dif,fil1,n_rec,kf,kr,t_tot,t_inc,valence,
     x                prl,delta_max,flo_rt,saturate,recfl)


c--------------------------------------------------------
c      Column Run:
c        calculate compositions           (ITER_KIN)
c        advance proteins                 (MIGRATE)
c        passage monitor queue            (MONITOR)
c        calculate terminus absorbance    (ABSORBANCE)
c        output elution status            (REPORT)
c        diffusion cycles                 (DIFFUSION_1)
c        output sums; close files         (SUMS)   
c----------------------------------------------------------       
 


      do 100  ncycl = 1,rl
        
         call ITER_KIN(1,p1,p2,pc,mw,kf,kr,delta_max,t_inc,
     x                 t_tot,1,cl,valence,mess_t_inc)
         call MIGRATE(prn,p1,p2,sprd,cl,ss,vl,dis,gs,
     x                saturate,pc,mw)
         call MONITOR(prn,p1,p2,mqt,cl,mq)
         call ABSORBANCE(prn,mqt,prel,prel_t,eluted,eluted_t,
     x                  mw,ec,mq,kf)
         call REPORT(ncycl,prel,prel_t,mqt,cl,ss,p1,pc,mw,
     x               count,rl,t_tot,prn,eluted,eluted_t,
     x               cv, flo_rt,
     x               position_label,saturate,valence)
         call DIFFUSION_1(ndif,p1,p2,ncycl,prn,cl+ss,dif)


 100  continue                    

         call SUMS (eluted,eluted_t,mess_t_inc,t_inc)


         end do         ! n_run
c-----------------------------------------------------------------
                    STOP
                    END


c------------------------------------------------------------------
c------------------------------------------------------------------
c------------------------------------------------------------------


       subroutine ABSORBANCE(prn,mqt,prel,prel_t,eluted,eluted_t,
     x                       mw,ec,mq,kf)

            !
            !  Subroutine to convert molar concentration to absorbance
            !  units at elution
            !
            !  k:      component i.d.
            !  mqt:    length of monitor queue
            !  prel:   protein eluted         
            !  mw:     molecular weights
            !  ec:     extinction coefficients
            !  mq:     sample concentrations in monitor queue
            !
            !     01/24/93
            !
            !

       integer mqt,prn
       real mq(0:200,0:15),mw(0:15),ec(0:15),kf(0:3),eluted(0:15)
       double precision prel(0:15)


       do 100 k = 0, prn
        prel(k)=mq(mqt,k)*mw(k)*ec(k)
        prel(k)=int(1000000*prel(k)+.5)/1000000.
 
 100   continue


     
         prel_t = 0.
         do 155 k = 0, prn 
155         prel_t = prel_t + prel(k)

           do 160 j=0,prn
            eluted(j) = eluted(j)+prel(j)
160        continue

           eluted_t = eluted_t + prel_t  

                end
   


c-----------------------------------------------------------------------------


         subroutine COLUMN(column_len,column_vol,mon_que,samp_size,
     x                  run_cycles,pre_cycles,t_inc,t_inc_set,
     x                  delta_max,t_tot,run_time,num_dif_cyc,sprd,
     x                  flo_rate,saturate,position_label,infile_col)
        
               !
               ! Routine opens file (column_in.dat) and obtains values
               ! which are fixed as experimental parameters  throughout 
               ! a series of iterative simulations
               ! 
               !         07/12/93
               !
               !

          character*5 position_label
          character*20 infile_col
          integer column_len,samp_size,run_cycles,pre_cycles,sprd
          double precision  t_inc
                    

          open(unit=03,file=infile_col,status='old')

             read(03,100) column_len  !number of cells in column
             read(03,110) column_vol  !column volume (ml)
             read(03,120) mon_que     !number of cells in monitor queu
             read(03,130) samp_size   !sample size in cells          
             read(03,140) run_cycles  !number of cycles in run
             read(03,140) pre_cycles  !number of cycles pre-run
             read(03,180) t_inc       !micro time step
             read(03,180) delta_max   !max % change during micro time step
             read(03,190) t_tot       !seconds per trans cycle
             read(03,160) num_dif_cyc !number of diffusion cycles/iter
             read(03,200) sprd        !gaussian spread (cells) translation
             read(03,170) flo_rate    !ml/min
             read(03,120) saturate   !if 1, then saturate column with B2

            column_len = column_len + samp_size

          close(unit=03)


          position_label = 'Mins.'
          if (flo_rate .ne. 0)   position_label = 'Vol. '
          if (column_vol .ne. 0) position_label = 'Ve/Vt'
          if ((column_vol .ne. 0) .and. 
     x         (flo_rate .eq. 0)) position_label = 'Mins.'
          t_inc_set = t_inc

          
100        format(i4)            !column_len
110        format(f7.4)          !column_vol
120        format(i3)            !mon_que
130        format(i3)            !samp_size
140        format(i5)            !run_cycles
150        format(f7.3)          !run_time
160        format(i2)            !num_dif_cyc
170        format(f8.4)          !flow_rate
180        format(f7.3)          !time increment
190        format(f7.3)          !time/cycle
200        format(i2)            !gaussian spread


                      end





c-----------------------------------------------------------------------------


          subroutine DIFFUSION_1(ndif,p1,p2,ncycl,prn,prfr,dif)

       !    routine to simulate diffusion of protein between
       !    adjacent cells in column

       !     ncycle:     number of cycle
       !     ndif:       number of diffusion cycles
       !     prn:        number of species
       !     prfr:       position of front edge
       !     dif:        diffusion factors
       !
       !         09/11/91
       !
       !
      
        double precision p1(-10:6000,0:15),p2(0:6000,0:15)
        integer ncycl,prn,ndif,prfr,itop
        real dif(0:15)

        do 100 kd = 1,ndif
       ndif1=ndif+1
       itop=prfr+ndif

       if (ncycl .ge. ndif) then

        do  k=0, prn

          do  j=1,itop
            p2(j,1)=p1(j+1,k)-p1(j,k)
            p1(j,k)=p1(j,k)+dif(k)*p2(j,1)-dif(k)*p2(j-1,1)
            p2(j-1,1)=0.
          end do 
 
         p2(itop,1)=0.
        end do
 
       end if

 100   continue

              end



c------------------------------------------------------------------------------


        subroutine EXTINCTION(mw,ec,valence)
  
c       
c        Calculates extinction coefficients of complexes
c        Calculates molecular weights of complexes
c       
c                Molecule 0:   A:A   (if valence of A for B1 = 1)
c                Molecule 0:   B1:A:B1  (if valence = 2)
c                Molecule 1:   A
c                Molecule 2:   B1
c                Molecule 3:   B2
c                Molecule 4:   A:B1
c                Molecule 5:   A:B2
c                Molecule 6:   (A:A)n
c
c       
c              06/14/93

       integer valence
       real mw(0:15),ec(0:15)
       

       if (valence .eq. 1) then  
           mw(0) = 2.*mw(1)
          else
           mw(0) = mw(1) + 2.*mw(2)
       end if
         mw(4) = mw(1) + mw(2)
         mw(5) = mw(1) + mw(3)
         mw(6) = 2.*mw(0)

       if (valence .eq. 2) then
         ec(0) = (mw(1)*ec(1) + 2*mw(2)*ec(2))
     x             / (mw(1) + 2.*mw(2))
         else
            ec(0) = ec(1)
       end if
         ec(4) = (mw(1)*ec(1)+mw(2)*ec(2))/mw(4)
         ec(5) = (mw(1)*ec(1)+mw(3)*ec(3))/mw(5)
         ec(6) = ec(0)
         
           do j=7,15
                mw(j) = float(j-4)*mw(0)
                ec(j) = ec(0)
           end do

                 end




c-----------------------------------------------------------------------------


        subroutine FILE_NUM(n_rec,fil1,label)

               !
               !
               !  Generates a unique output file name by
               !  reading the file 'ksym_rec.dat' which contains
               !  a single integer.  The value of this integer
               !  is increased by 1, and the new value is written
               !  to the file.  This number is used as part of the
               !  name of files containing output data from this
               !  run.
               !
               !      09/09/91
               !
 
       character*16 rec_file      !record file name
       character*4 fil1           !plot file prefix
       character*13 pltfl         !plotting file
       character*12 recfl         !documentation file   
       real label(10)
      
        data rec_file/'ksym_rec.dat'/
       open (unit=04, file = rec_file, status = 'old')
              
          read(04,*) n_rec 
       close (unit=04)
          n_rec = n_rec + 1
       open (unit=04, file = rec_file, status = 'old')
          write(04,*) n_rec
       close (unit=04)
          label(n_run) = n_rec

       

            return
             end
 


c----------------------------------------------------------------------------


        subroutine GAUSS(gs)
  
          !
          !  Routine loads gaussian distribution values to be
          !  used by translation routine
          !
          !     08/26/91
          !


        real gs(705),valu1(180),valu2(171),value(351)

       data valu1/.5000,.5040,.5080,.5120,.5160,
     x .5199,.5239,.5279,.5319,.5359,
     x .5398,.5438,.5478,.5517,.5557,.5596,.5636,.5675,.5714,.5753,
     x .5793,.5832,.5871,.5910,.5948,.5987,.6026,.6064,.6103,.6141,
     x .6179,.6217,.6255,.6293,.6331,.6368,.6406,.6443,.6480,.6517,
     x .6554,.6591,.6628,.6664,.6700,.6736,.6772,.6808,.6844,.6879,
     x .6915,.6950,.6985,.7019,.7054,.7088,.7123,.7157,.7190,.7224,
     x .7257,.7291,.7324,.7357,.7389,.7422,.7454,.7486,.7517,.7549,
     x .7580,.7611,.7642,.7673,.7704,.7734,.7764,.7794,.7823,.7852,
     x .7881,.7910,.7939,.7967,.7995,.8023,.8051,.8078,.8106,.8133,
     x .8159,.8186,.8212,.8238,.8264,.8289,.8315,.8340,.8365,.8389,
     x .8413,.8438,.8461,.8485,.8508,.8531,.8554,.8577,.8599,.8621,
     x .8643,.8665,.8686,.8708,.8729,.8749,.8770,.8790,.8810,.8830,
     x .8849,.8869,.8888,.8907,.8925,.8944,.8962,.8980,.8997,.9015,
     x .9032,.9049,.9066,.9082,.9099,.9115,.9131,.9147,.9162,.9177,
     x .9192,.9207,.9222,.9236,.9251,.9265,.9279,.9292,.9306,.9319,
     x .9332,.9345,.9357,.9370,.9382,.9394,.9406,.9418,.9429,.9441,
     x .9452,.9463,.9474,.9484,.9495,.9505,.9515,.9525,.9535,.9545,
     x .9554,.9564,.9573,.9582,.9591,.9599,.9608,.9616,.9625,.9633/
c
         data valu2
     x /.9641,.9649,.9656,.9664,.9671,.9678,.9686,.9693,.9699,.9706,
     x .9713,.9719,.9726,.9732,.9738,.9744,.9750,.9756,.9761,.9767,
     x .9772,.9778,.9783,.9788,.9793,.9798,.9803,.9808,.9812,.9817,
     x .9821,.9826,.9830,.9834,.9838,.9842,.9846,.9850,.9854,.9857,
     x .9861,.9864,.9868,.9871,.9875,.9878,.9881,.9884,.9887,.9890,
     x .9893,.9896,.9898,.9901,.9904,.9906,.9909,.9911,.9913,.9916,
     x .9918,.9920,.9922,.9925,.9927,.9929,.9931,.9932,.9934,.9936,
     x .9938,.9940,.9941,.9943,.9945,.9946,.9948,.9949,.9951,.9952,
     x .9953,.9955,.9956,.9957,.9959,.9960,.9961,.9962,.9963,.9964,
     x .9965,.9966,.9967,.9968,.9969,.9970,.9971,.9972,.9973,.9974,
     x .9974,.9975,.9976,.9977,.9977,.9978,.9979,.9979,.9980,.9981,
     x .9981,.9982,.9982,.9983,.9984,.9984,.9985,.9985,.9986,.9986,
     x .9987,.9987,.9987,.9988,.9988,.9989,.9989,.9989,.9990,.9990,
     x .9990,.9991,.9991,.9991,.9992,.9992,.9992,.9992,.9993,.9993,
     x .9993,.9993,.9994,.9994,.9994,.9994,.9994,.9995,.9995,.9995,
     x .9995,.9995,.9995,.9996,.9996,.9996,.9996,.9996,.9996,.9997,
     x .9997,.9997,.9997,.9997,.9997,.9997,.9997,.9997,.9997,.9998,
     x 1.000/



       do 300 i=1,180
300        value(i)=valu1(i)
 
        j=180
       do 301 i=1,171
        j=j+1
        value(j)=valu2(i)
301     continue
701    format(10f7.4)


       j=350
       do 302 i=1,351
        j=j+1
       gs(j)=value(i)
302    continue
 
        k=702
        do 3 i=1 , 350
            m=k-i
            gs(i)=1-gs(m)
3       continue
 
           end


c----------------------------------------------------------------------------

       subroutine ITER_KIN(n_cycles,p1,p2,pc,mw,kf,kr,delta_max,
     x         t_inc,t_tot,cell1,cell2,valence,mess_t_inc)


           !
           ! p1(j,0) == A:A         cell j
           ! p1(j,1) == A           cell j
           ! p1(j,2) == B1          cell j
           ! p1(j,3) == B2          cell j
           ! p1(j,4) == A:B1        cell j
           ! p1(j,5) == A:B2        cell j
           ! p1(j,6) == (A:A)(A:A)  cell j 
           !

           integer cell1,cell2,valence,mess_t_inc
           real kf(0:3),kr(0:3),mw(0:15),pc(0:15),p_old(0:15),
     x          c_ref,dt_max,delta_max

           double precision p1(-10:6000,0:15),del_p1(15),fr_free,
     x        del_complx,complx,bound_frac,p1_tot,total_a,t_inc



c
c        Check for appropriate t_inc:  change in concentration to be less
c        than   delta_max    of lowest concentration  component A, B1, or B2.
c

            if (n_cycles .gt. 1) then

               dt_max = t_inc

              if( (kf(0) .ne. 0) .and. (p1(1,1) .ne. 0) ) then
                 dt_max = delta_max*( 1./(kf(0)*p1(1,1)) )
              end if

              if ( (kf(1) .ne. 0) .and. (p1(1,1)*p1(1,2) .ne. 0) )
     x             then
                 c_ref = p1(1,1)
                   if (p1(1,2) .lt. c_ref) c_ref = p1(1,2)
                 dt_max = delta_max*( c_ref/(kf(0)*p1(1,1)*p1(1,2)) ) 
              end if

              if ( (kf(2) .ne. 0) .and. (p1(1,1)*p1(1,3) .ne. 0) )
     x             then
                 c_ref = p1(1,1)
                   if (p1(1,3) .lt. c_ref) c_ref = p1(1,3)
                 dt_max = delta_max*( c_ref/(kf(0)*p1(1,1)*p1(1,3)) ) 
              end if

              if (dt_max .lt. t_inc) then 
                 mess_t_inc = 1
                 t_inc = dt_max
                 write(6,850) t_inc
              end if

            end if

           n_reps = t_tot/t_inc




          do 200 n = 1,n_cycles
            do 150 j = cell1,cell2


             do 100 k = 1, n_reps


              p_tot = mw(0)*p1(j,0) + mw(1)*p1(j,1)
     x              +  mw(2)*p1(j,2) + mw(3)*p1(j,3)
     x              + (mw(4))*p1(j,4)
     x              + (mw(5))*p1(j,5)


c
              if (p_tot .lt. .0001*(pc(1)+pc(2))) then 
                  goto 105
                  end if


c         Calculate interaction between  molecule A [p1(j,1)] with
c         itself to form dimer A:A [p1(j,0)]
         
c

           if (kf(0) .gt. 0)  then
            del_complx = t_inc*( kf(0)*p1(j,1)*p1(j,1) 
     x                     - kr(0)*p1(j,0) )

               p1(j,0) = p1(j,0) + del_complx
               p1(j,1) = p1(j,1) - 2.*del_complx

           end if







c         Calculate interaction between  molecule A [p1(j,1)] and 
c         molecule B1 [p1(j,2)] to form complex A:B1 [p1(j,4)]
         

           if (kf(1) .gt. 0)  then


            del_complx = (t_inc/float(valence))
     x                *( kf(1)*float(valence)*p1(j,1)
     x                         *p1(j,2) - kr(1)*p1(j,4) )

               p1(j,4) = p1(j,4) + del_complx
               p1(j,1) = p1(j,1) - del_complx
               p1(j,2) = p1(j,2) - del_complx



c        Calculate bivalent complex B1:A:B1 [p1(j,0)]

         if(valence .eq. 2) then

            del_complx = (t_inc/float(valence))*
     x                   ( kf(1)*p1(j,2)*p1(j,4) 
     x                     - 2.*kr(1)*p1(j,0))


              p1(j,0) = p1(j,0) + del_complx
              p1(j,2) = p1(j,2) - del_complx
              p1(j,4) = p1(j,4) - del_complx

          end if

 300      continue

           end if



c
c         Calculate interaction between  molecule A [p1(j,1)] and 
c         molecule B2 [p1(j,3)] to form complex A:B2 [p1(j,5)]
c
            
          if (kf(2) .gt. 0)  then
            del_complx = t_inc*( kf(2)*p1(j,1)*p1(j,3) 
     x                     - kr(2)*p1(j,5) )

               p1(j,5) = p1(j,5) + del_complx
               p1(j,1) = p1(j,1) - del_complx
               p1(j,3) = p1(j,3) - del_complx

          end if


            
 105            continue

c
c         Calculate interaction between  dimer A:A [p1(j,0)] with
c         itself to form higher order aggregates A:A [p1(j,6)] ...
c         
 

           if (kf(3) .gt. 0)  then

            del_complx = t_inc*( kf(3)*p1(j,0)*p1(j,0) 
     x                     - kr(3)*p1(j,6) )

               p1(j,6) = p1(j,6) + del_complx
               p1(j,0) = p1(j,0) - 2.*del_complx

c
c       Calculate change in concentration of 3-mer: D-D-D where D represents an
c       A-A dimer, D-D-D == p1(j,7);  note that p(7) can break down
c       to p(0) at either of two sites, hence the factor of two in
c       the dissociation term below.
c

            del_complx = t_inc*( kf(3)*p1(j,0)*p1(j,6) 
     x                     - 2.*kr(3)*p1(j,7) )

               p1(j,7) = p1(j,7) + del_complx
               p1(j,0) = p1(j,0) - del_complx
               p1(j,6) = p1(j,6) - del_complx


c
c
c       Calculate change in concentration of p(8): D-D-D-D from association
c       of p(0) and p(7).  Breakdown of p(8) can occur with equal likelihood
c       at any of three sites resulting in formation of equimolar units of
c       p(0), p(6), and p(7).   Three moles of p(8) required to form
c       two moles of break-down products.
c
               del_formed = t_inc*kf(3)*p1(j,0)*p1(j,7)
               del_break  = 3.*t_inc*kr(3)*p1(j,8)

               p1(j,8) = p1(j,8) + del_formed - del_break
               p1(j,0) = p1(j,0) - del_formed + 2.*del_break/3.
               p1(j,6) = p1(j,6) + 2*del_break/3.
               p1(j,7) = p1(j,7) - del_formed + 2.*del_break/3.




c
c       Calculate change in concentration of p(9): D^5 from association
c       of p(0) and p(8).  Breakdown of p(9) can occur with equal likelihood
c       at any of four sites resulting in formation of equimolar units of
c       p(0), p(6), p(7), and p(8).   Four moles of p(9) required to form
c       two moles of break-down products.
c
               del_formed = t_inc*kf(3)*p1(j,0)*p1(j,8)
               del_break  = 4.*t_inc*kr(3)*p1(j,9)

               p1(j,9) = p1(j,9) + del_formed - del_break
               p1(j,0) = p1(j,0) - del_formed + 2.*del_break/4.
               p1(j,6) = p1(j,6) + 2.*del_break/4.
               p1(j,7) = p1(j,7) + 2.*del_break/4.
               p1(j,8) = p1(j,8) - del_formed + 2.*del_break/4.




c
c       Calculate change in concentration of p(10): D^6 from association
c       of p(0) and p(9).  Breakdown of p(10) can occur with equal likelihood
c       at any of five sites resulting in formation of equimolar units of
c       p(0), p(6) -- p(9).   Five moles of p(10) are required to form
c       two moles of break-down products.
c
               del_formed = t_inc*kf(3)*p1(j,0)*p1(j,9)
               del_break  = 5.*t_inc*kr(3)*p1(j,10)

               p1(j,10) = p1(j,10) + del_formed - del_break
               p1(j,0) = p1(j,0) - del_formed + 2.*del_break/5.
               p1(j,6) = p1(j,6) + 2.*del_break/5.
               p1(j,7) = p1(j,7) + 2.*del_break/5.
               p1(j,8) = p1(j,8) + 2.*del_break/5.
               p1(j,9) = p1(j,9) - del_formed + 2.*del_break/5.






c
c       Calculate change in concentration of p(11): D^7 from association
c       of p(0) and p(10).  Breakdown of p(11) can occur with equal likelihood
c       at any of six sites resulting in formation of equimolar units of
c       p(0), p(6) -- p(10).   Six moles of p(11) are required to form
c       two moles of break-down products.
c
               del_formed = t_inc*kf(3)*p1(j,0)*p1(j,10)
               del_break  = 6.*t_inc*kr(3)*p1(j,11)

               p1(j,11) = p1(j,11) + del_formed - del_break
               p1(j,0) = p1(j,0) - del_formed + 2.*del_break/6.
               p1(j,6) = p1(j,6) + 2.*del_break/6.
               p1(j,7) = p1(j,7) + 2.*del_break/6.
               p1(j,8) = p1(j,8) + 2.*del_break/6.
               p1(j,9) = p1(j,9) + 2.*del_break/6.
               p1(j,10) = p1(j,10) - del_formed + 2.*del_break/6.






c
c       Calculate change in concentration of p(12): D^8 from association
c       of p(0) and p(11).  Breakdown of p(12) can occur with equal likelihood
c       at any of seven sites resulting in formation of equimolar units of
c       p(0), p(6) -- p(11).   Seven moles of p(12) are required to form
c       two moles of break-down products.
c
               del_formed = t_inc*kf(3)*p1(j,0)*p1(j,11)
               del_break  = 7.*t_inc*kr(3)*p1(j,12)

               p1(j,12) = p1(j,12) + del_formed - del_break
               p1(j,0) = p1(j,0) - del_formed + 2.*del_break/7.
               p1(j,6) = p1(j,6) + 2.*del_break/7.
               p1(j,7) = p1(j,7) + 2.*del_break/7.
               p1(j,8) = p1(j,8) + 2.*del_break/7.
               p1(j,9) = p1(j,9) + 2.*del_break/7.
               p1(j,10) = p1(j,10) + 2.*del_break/7.
               p1(j,11) = p1(j,11) - del_formed + 2.*del_break/7.






c
c       Calculate change in concentration of p(13): D^9 from association
c       of p(0) and p(12).  Breakdown of p(13) can occur with equal likelihood
c       at any of eight sites resulting in formation of equimolar units of
c       p(0), p(6) -- p(12).   Eight moles of p(13) are required to form
c       two moles of break-down products.
c
               del_formed = t_inc*kf(3)*p1(j,0)*p1(j,12)
               del_break  = 7.*t_inc*kr(3)*p1(j,13)

               p1(j,13) = p1(j,13) + del_formed - del_break
               p1(j,0) = p1(j,0) - del_formed + 2.*del_break/8.
               p1(j,6) = p1(j,6) + 2.*del_break/8.
               p1(j,7) = p1(j,7) + 2.*del_break/8.
               p1(j,8) = p1(j,8) + 2.*del_break/8.
               p1(j,9) = p1(j,9) + 2.*del_break/8.
               p1(j,10) = p1(j,10) + 2.*del_break/8.
               p1(j,11) = p1(j,11) + 2.*del_break/8.
               p1(j,12) = p1(j,12) - del_formed + 2.*del_break/8.






           end if



 100            continue
 150            continue



             if (n .eq. n_cycles-1)  then
                do j = 0,13
                   p_old(j) = p1(1,j)
                end do
             end if


             if ((n .eq. n_cycles) .and.
     x          (n_cycles .gt. 1)) then

                   p_temp = 0
                   p_temp_old = 0
                do j = 6, 13
                   p_temp = p_temp + mw(j)*p1(1,j)
                   p_temp_old = p_temp_old + mw(j)*p_old(j)
                end do


                 time = t_inc*float(n)*float(n_reps)
                 if (valence .eq. 1) then
                   write(6,810)
                   write(7,810)
                 else
                   write(6,812)
                   write(7,812)
                 end if



                 write(6,800) time, (mw(m)*p1(1,m),m=0,5), p_temp
                 write(7,800) time, (mw(m)*p1(1,m),m=0,5), p_temp
                 write(6,815)      (mw(m)*p_old(m),m=0,5), p_temp_old
                 write(7,815)      (mw(m)*p_old(m),m=0,5), p_temp_old
                 write(6,820) (mw(m)*(p1(1,m)-p_old(m)),m=0,5),
     x               (p_temp - p_temp_old)
                 write(7,820) (mw(m)*(p1(1,m)-p_old(m)),m=0,5),
     x               (p_temp - p_temp_old)


c
c          Compute complexes to check proper results from iteration
c
c      Case 1:  Mono- or bivalent interaction between A and B1
c
             if ((kf(1) .ne. 0) .and. (kf(0) .eq. 0)) then
                 term_1 = float(valence)*pc(1)/mw(1) 
     x              + pc(2)/mw(2) + kr(1)/kf(1)
                 term_2 = term_1/2.
                 term_3 = float(valence)*(pc(1)/mw(1))*(pc(2)/mw(2))
                 term_4 = term_2*term_2 - term_3
                 complex = term_2 - sqrt(term_4)


             if (valence .eq. 2) then
               fr_bound = complex/(float(valence)*(pc(1)/mw(1)))
               b_a_b = fr_bound*fr_bound*(pc(1)/mw(1))
               b_a_b = b_a_b*mw(0)

               a_ = (1-fr_bound)*(1-fr_bound)*(pc(1)/mw(1))
               a_ = a_*mw(1)

               a_b = 2*fr_bound*(1-fr_bound)*(pc(1)/mw(1))
               a_b = a_b*mw(4)

               b_ = (pc(2)/mw(2)) - complex
               b_ = b_ * mw(2)
 
               write(6,830) b_a_b,a_,b_,a_b
               write(7,830) b_a_b,a_,b_,a_b
             
              else
               
               a_ = (pc(1)/mw(1)) - complex
               a_ = a_ * mw(1)

               b_ = (pc(2)/mw(2)) - complex
               b_ = b_ * mw(2)

               a_b = complex * mw(4)
               b_a_b = 0

               write(6,830) b_a_b,a_,b_,a_b
               write(7,830) b_a_b,a_,b_,a_b

              endif

             endif



c
c    Case 2:  Monovalent interaction between A and B2
c
            if ((kf(0) .eq. 0) .and. (kf(1) .eq. 0) .and.
     x          (kf(2) .ne. 0) .and. (kf(3) .eq. 0)) then


                 term_1 = pc(1)/mw(1) 
     x              + pc(3)/mw(3) + kr(2)/kf(2)
                 term_2 = term_1/2.
                 term_3 = (pc(1)/mw(1))*(pc(3)/mw(3))
                 term_4 = term_2*term_2 - term_3
                 complex = term_2 - sqrt(term_4)

               
               a_ = (pc(1)/mw(1)) - complex
               a_ = a_ * mw(1)

               b_ = (pc(3)/mw(3)) - complex
               b_ = b_ * mw(3)

               a_b = complex * mw(5)

               b_a_b = 0

               write(6,835) b_a_b,a_,b_,a_b
               write(7,835) b_a_b,a_,b_,a_b
             end if

c
c    Case 3  Competitive interactions involving A, B1, and B2
c
c
          if ((kf(0) .eq. 0)  .and.  (kf(1) .ne. 0) .and.
     x        (kf(2) .ne. 0)  .and.  (kf(3) .eq. 0))  then

        write (6,*)  'Expected:      ! Not yet implemented  !'
        write (7,*)  'Expected:      ! Not yet implemented  !'


          end if
c 
c
c



c
c    Case 4  Self-association of A
c

        if ( (kf(0) .ne. 0) .and. (kf(1) .eq. 0) .and.
     x       (kf(2) .eq. 0) .and. (kf(3) .eq. 0) )  then

                 term_1 = pc(1)/mw(1) + (kr(0)/kf(0))/4.
                 term_2 = term_1/2.
                 term_3 = (pc(1)/mw(1))*(pc(1)/mw(1))
                 term_4 = term_2*term_2 - term_3/4.
                 complex = term_2 - sqrt(term_4)

               a_ = (pc(1)/mw(1)) - 2.*complex
               a_ = a_ * mw(1)

               a_a = complex * mw(0)

               write(6,835) a_a,a_
               write(7,835) a_a,a_

            end if

c
c
c       Case 5  Self-association and oligomerization of A.
c
c

            if ((kf(0) .ne. 0) .and. (kf(1) .eq. 0) .and.
     x          (kf(2) .eq. 0) .and. (kf(3) .ne. 0))  then

c
        write (6,*)  'Expected:      ! Not yet implemented  !'
        write (7,*)  'Expected:      ! Not yet implemented  !'
c

              end if


              endif   ! terminates if loop for post-equilibration output


 200        continue



 800    format(2x,f7.1,' s',3x,7f9.4,3x,f9.4)
 810    format(///,3x,'After pre-equilibration:',//,5x,
     x          'Time',9x,'A:A',8x,'A',8x,'B1',7x,'B2',5x,'A:B1',
     x          5x,'A:B2',4x,'(A:A)n')
 812    format(///,3x,'After pre-equilibration:',//,5x,
     x          'Time',8x,'B:A:B',7x,'A',8x,'B1',7x,'B2',5x,'A:B1',
     x          5x,'A:B2',4x,'(A:A)n')
 815    format(4x,'Last:    ',1x,7f9.4,3x,f9.4)
 820    format(4x,'Delta:   ',1x,7f9.4,3x,f9.4,/)
 830    format(4x,'Expected:'1x,3f9.4,9x,f9.4)
 835    format(4x,'Expected:'1x,2f9.4,9x,f9.4,9x,f9.4)
 850    format(//,4x,'Note:  The value of the parameter < t_inc > was',
     x   ' changed to',/,'  [', f12.8, ' ] seconds ',
     x   'to satisfy the criterion set by < delta_max >.')


                            end





c-----------------------------------------------------------------------------




        subroutine LOADER(ss,cl,mq,mqt,pc,saturate,mw,p1,valence)

         !
         !  routine to calculate molar concentration and
         !           to load sample
         !

         ! k1:      component #
         ! pc:      concentration of component k1 (MG/Ml)
         ! mw:      molecular weight of component k1
         ! prload:  molar concentration
         ! ss:      sample size
         ! p1:      concentration of sample k1 in cell k2
         ! k2:      cell #
         !
         !         06/11/93
         !

       integer k2,ss,cl,saturate,valence
       double precision p1(-10:6000,0:15)
       real pc(0:5),mw(0:15),prload,mq(0:200,0:15)




       do 100  k=1,3
        prload=pc(k)/mw(k)


       do  k2=1,ss
          p1(k2,k)=prload
       end do

100   continue
      

      write(6,900)
 900  format(/,/,'  Initiating Pre-Equilibration  ',/,/)

                end



c-----------------------------------------------------------------------------




       subroutine LOG_DAT(infile,sprd,prn,mqt,cl,cv,rl,ss,ndif,
     x   vl,dis, pc,mw,ec,dif,fil1,n_rec,kf,kr,t_tot,t_inc,
     x   valence,prl,delta_max,flo_rt,saturate,recfl)

              !
              ! Updates chronicle (log) file for cross-
              ! referencing of run paramaters with output
              ! files.
              !
              !           08/28/93
              !
              !
        parameter (max_records = 200)
        parameter (lists = 4)

        character*1 scratch
        character*3 satB2(max_records)
        character*4 fil1
        character*8 i_time, times(max_records)
        character*9 n_day, days(max_records)
        character*10 label(0:74)
        character*12 recfl
        character*14 infiles(max_records),outfiles(max_records)
        character*25 infile


        integer ss,cl,cls,mqt,ncycl,ndif,prn,rl,sprd,valence
        integer records, lists, prl, saturate

        real kf(0:3),kr(0:3),pc(0:15),vl(0:15),dis(0:15),ec(0:15),
     x   mw(0:15),dif(0:15)

         real pcA(max_records), mwA(max_records), vlA(max_records),
     x     ecA(max_records), dfA(max_records), dsA(max_records),
     x     pcB1(max_records), mwB1(max_records), vlB1(max_records),
     x     ecB1(max_records), dfB1(max_records), dsB1(max_records),
     x     pcB2(max_records), mwB2(max_records), vlB2(max_records),
     x     ecB2(max_records), dfB2(max_records), dsB2(max_records),
     x     vlAA(max_records), dfAA(max_records), dsAA(max_records),
     x     vlAB1(max_records),dfAB1(max_records),dsAB1(max_records),
     x     vlABB(max_records),dfABB(max_records),dsABB(max_records),
     x     vlAB2(max_records),dfAB2(max_records),dsAB2(max_records),
     x     vlC2(max_records), dfC2(max_records), dsC2(max_records),
     x     vlC3(max_records), dfC3(max_records), dsC3(max_records),
     x     vlC4(max_records), dfC4(max_records), dsC4(max_records),
     x     vlC5(max_records), dfC5(max_records), dsC5(max_records),
     x     vlC6(max_records), dfC6(max_records), dsC6(max_records),
     x     vlC7(max_records), dfC7(max_records), dsC7(max_records),
     x     kf_AA(max_records),kr_AA(max_records),kf_AB1(max_records),
     x     kr_AB1(max_records),kf_AB2(max_records),kr_AB2(max_records),
     x     kf_CC(max_records),kr_CC(max_records),tim_cyc(max_records),
     x     tim_inc(max_records),del_max(max_records),
     x     colvols(max_records),florts(max_records)

       integer clen(max_records),smpsz(max_records),monq(max_records),
     x     rnln(max_records),spread(max_records),difcyc(max_records),
     x     valnc(max_records),prls(max_records)

          
        double precision t_inc


        data label/'          ','Date:     ','Time:     ','Infile:   ',
     x             'A conc:   ','A mw:     ','A veloc:  ','A extcf:  ',
     x             'A dif:    ','A disp:   ','B1 conc:  ','B1 mw:    ',
     x             'B1 veloc: ','B1 extcf: ','B1 dif:   ','B1 disp:  ',
     x             'B2 conc:  ','B2 mw:    ','B2 veloc: ','B2 extcf: ',
     x             'B2 dif:   ','B2 disp:  ','AA veloc: ','AA dif:   ',
     x             'AA disp:  ','AB1 velc: ','AB1 dif:  ','AB1 disp: ',
     x             'AB2 velc: ','AB2 dif:  ','AB2 disp: ','ABB velc: ',
     x             'ABB dif:  ','ABB disp: ','C2 veloc: ','C2 dif:   ',
     x             'C2 disp:  ','C3 veloc: ','C3 dif:   ','C3 disp:  ',
     x             'C4 veloc: ','C4 dif:   ','C4 disp:  ','C5 veloc: ',
     x             'C5 dif:   ','C5 disp:  ','C6 veloc: ','C6 dif:   ',
     x             'C6 disp:  ','C7 veloc: ','C7 dif:   ','C7 disp:  ',
     x             'kf_AA:    ','kr_AA:    ','kf_AB1:   ','kr_AB1:   ',
     x             'kf_AB2:   ','kr_AB2:   ','kf_CC:    ','kr_CC:    ',
     x             'Col_len:  ','SampSz:   ','MonQ:     ','Run_len:  ',
     x             'Time_cyc: ','Time_inc: ','MaxDelta: ','GausSprd: ',
     x             'Dif_cycs: ','ValenceA: ','PreEquil: ','Saturate: ',
     x             'ColVl(ml):','FloRate:  ','Outfile:  '/

c
c     Read existing log file
c
        open (unit=09, file = 'scm10_log.dat',status = 'old')
          read(09,*) records

           if (records .gt. 0) then

          do 100 j = 1, 1 + records/lists
             
             n1 = 1 + (j-1)*lists
             n2 = j*lists
             if (n2 .gt. records) n2 = records
           read(09,810) scratch
           read(09,810) scratch
           read(09,*) label(0),(days(k), k = n1,n2)  
           read(09,*) label(0),(times(k), k = n1,n2) 
           read(09,*) label(0),(infiles(k), k = n1,n2) 
           read(09,800) label(0),(outfiles(k),k = n1,n2)

           read(09,820) label(0),(pcA(k), k = n1,n2)
           read(09,825) label(0),(mwA(k), k = n1,n2)
           read(09,820) label(0),(vlA(k), k = n1,n2)
           read(09,820) label(0),(ecA(k), k = n1,n2)
           read(09,820) label(0),(dfA(k), k = n1,n2)
           read(09,820) label(0),(dsA(k), k = n1,n2)
           read(09,840) label(0),(valnc(k), k = n1,n2)

           read(09,820) label(0),(pcB1(k), k = n1,n2)
           read(09,825) label(0),(mwB1(k), k = n1,n2)
           read(09,820) label(0),(vlB1(k), k = n1,n2)
           read(09,820) label(0),(ecB1(k), k = n1,n2)
           read(09,820) label(0),(dfB1(k), k = n1,n2)
           read(09,820) label(0),(dsB1(k), k = n1,n2)

           read(09,820) label(0),(pcB2(k), k = n1,n2)
           read(09,825) label(0),(mwB2(k), k = n1,n2)
           read(09,820) label(0),(vlB2(k), k = n1,n2)
           read(09,820) label(0),(ecB2(k), k = n1,n2)
           read(09,820) label(0),(dfB2(k), k = n1,n2)
           read(09,820) label(0),(dsB2(k), k = n1,n2)
           read(09,800) label(0),(satB2(k), k = n1,n2)

           read(09,820) label(0),(vlAA(k), k = n1,n2)
           read(09,820) label(0),(dfAA(k), k = n1,n2)
           read(09,820) label(0),(dsAA(k), k = n1,n2)

           read(09,820) label(0),(vlAB1(k), k = n1,n2)
           read(09,820) label(0),(dfAB1(k), k = n1,n2)
           read(09,820) label(0),(dsAB1(k), k = n1,n2)

           read(09,820) label(0),(vlAB2(k), k = n1,n2)
           read(09,820) label(0),(dfAB2(k), k = n1,n2)
           read(09,820) label(0),(dsAB2(k), k = n1,n2)

           read(09,820) label(0),(vlABB(k), k = n1,n2)
           read(09,820) label(0),(dfABB(k), k = n1,n2)
           read(09,820) label(0),(dsABB(k), k = n1,n2)

           read(09,820) label(0),(vlC2(k), k = n1,n2)
           read(09,820) label(0),(dfC2(k), k = n1,n2)
           read(09,820) label(0),(dsC2(k), k = n1,n2)

           read(09,820) label(0),(vlC3(k), k = n1,n2)
           read(09,820) label(0),(dfC3(k), k = n1,n2)
           read(09,820) label(0),(dsC3(k), k = n1,n2)

           read(09,820) label(0),(vlC4(k), k = n1,n2)
           read(09,820) label(0),(dfC4(k), k = n1,n2)
           read(09,820) label(0),(dsC4(k), k = n1,n2)

           read(09,820) label(0),(vlC5(k), k = n1,n2)
           read(09,820) label(0),(dfC5(k), k = n1,n2)
           read(09,820) label(0),(dsC5(k), k = n1,n2)

           read(09,820) label(0),(vlC6(k), k = n1,n2)
           read(09,820) label(0),(dfC6(k), k = n1,n2)
           read(09,820) label(0),(dsC6(k), k = n1,n2)

           read(09,820) label(0),(vlC7(k), k = n1,n2)
           read(09,820) label(0),(dfC7(k), k = n1,n2)
           read(09,820) label(0),(dsC7(k), k = n1,n2)
           read(09,810) scratch

           read(09,830) label(0),(kf_AA(k), k = n1,n2)
           read(09,830) label(0),(kr_AA(k), k = n1,n2)

           read(09,830) label(0),(kf_AB1(k), k = n1,n2)
           read(09,830) label(0),(kr_AB1(k), k = n1,n2)

           read(09,830) label(0),(kf_AB2(k), k = n1,n2)
           read(09,830) label(0),(kr_AB2(k), k = n1,n2)

           read(09,830) label(0),(kf_CC(k), k = n1,n2)
           read(09,830) label(0),(kr_CC(k), k = n1,n2)

           read(09,810) scratch
           read(09,840) label(0),(clen(k), k = n1,n2)
           read(09,840) label(0),(smpsz(k), k = n1,n2)
           read(09,840) label(0),(monq(k), k = n1,n2)
           read(09,840) label(0),(rnln(k), k = n1,n2)
           read(09,840) label(0),(prls(k), k = n1,n2)           
           read(09,820) label(0),(tim_cyc(k),k = n1,n2)
           read(09,820) label(0),(tim_inc(k),k = n1,n2)
           read(09,820) label(0),(del_max(k),k = n1,n2)
           read(09,840) label(0),(spread(k),k = n1,n2)
           read(09,840) label(0),(difcyc(k),k = n1,n2)
           read(09,820) label(0),(colvols(k),k = n1,n2)
           read(09,820) label(0),(florts(k),k = n1,n2)


 100       continue

           end if
           close(unit=09)

c
c    Add current values to log file array
c


        call time(i_time)
        call date(n_day)
        records = records + 1

        days(records) = n_day
        times(records) = i_time
        infiles(records) = infile
        outfiles(records) = recfl
        pcA(records) = pc(1)
        mwA(records) = mw(1)
        vlA(records) = vl(1)
        ecA(records) = ec(1)
        dfA(records) = dif(1)
        dsA(records) = dis(1)
        valnc(records) = valence

        pcB1(records) = pc(2)
        mwB1(records) = mw(2)
        vlB1(records) = vl(2)
        ecB1(records) = ec(2)
        dfB1(records) = dif(2)
        dsB1(records) = dis(2)

        pcB2(records) = pc(3)
        mwB2(records) = mw(3)
        vlB2(records) = vl(3)
        ecB2(records) = ec(3)
        dfB2(records) = dif(3)
        dsB2(records) = dis(3)
        if(saturate .eq. 1) then
           satB2(records) = 'Yes'
          else
           satB2(records) = 'No '
        end if

        if (valence .ne. 2) then
          vlAA(records) = vl(0)
          dfAA(records) = dif(0)
          dsAA(records) = dis(0)
        end if

        if (valence .eq. 2) then
          vlAA(records) = 0
          dfAA(records) = 0
          dsAA(records) = 0
        end if

        vlAB1(records) = vl(4)
        dfAB1(records) = dif(4)
        dsAB1(records) = dis(4)

        vlAB2(records) = vl(5)
        dfAB2(records) = dif(5)
        dsAB2(records) = dis(5)

        if (valence .ne. 2) then
          vlABB(records) = 0
          dfABB(records) = 0
          dsABB(records) = 0
        end if

        if (valence .eq. 2) then
          vlABB(records) = vl(0)
          dfABB(records) = dif(0)
          dsABB(records) = dis(0)
        end if


        if (kf(3) .gt. 0) then

          vlC2(records) = vl(6)
          dfC2(records) = dif(6)
          dsC2(records) = dis(6)

          vlC3(records) = vl(7)
          dfC3(records) = dif(7)
          dsC3(records) = dis(7)

          vlC4(records) = vl(8)
          dfC4(records) = dif(8)
          dsC4(records) = dis(8)

          vlC5(records) = vl(9)
          dfC5(records) = dif(9)
          dsC5(records) = dis(9)

          vlC6(records) = vl(10)
          dfC6(records) = dif(10)
          dsC6(records) = dis(10)

          vlC7(records) = vl(11)
          dfC7(records) = dif(11)
          dsC7(records) = dis(11)
          
        end if

        kf_AA(records) = kf(0)
        kr_AA(records) = kr(0)

        kf_AB1(records) = kf(1)
        kr_AB1(records) = kr(1)

        kf_AB2(records) = kf(2)
        kr_AB2(records) = kr(2)

        kf_CC(records) = kf(3)
        kr_CC(records) = kr(3)

        clen(records)    = cl-ss
        smpsz(records)   = ss
        monq(records)    = mqt
        prls(records)    = prl
        rnln(records)    = rl
        tim_cyc(records) = t_tot
        tim_inc(records) = t_inc
        del_max(records) = delta_max
        spread(records)  = sprd
        difcyc(records)  = ndif
        colvols(records) = cv
        florts(records)  = flo_rt

c
c     Write log file array to new log file
c

        open (unit=09, file = 'scm10_log.dat',status = 'new')

        write(9,*) records

          do 200 j = 1, 1 + records/lists
             
             n1 = 1 + (j-1)*lists
             n2 = j*lists
             if (n2 .gt. records) n2 = records

           write(09,810) ' '
           write(09,810) ' '
           write(09,800) label(1),(days(k), k = n1,n2) 
           write(09,800) label(2),(times(k), k = n1,n2) 
           write(09,800) label(3),(infiles(k), k = n1,n2)
           write(09,800) label(74),(outfiles(k),k = n1,n2)

           write(09,820) label(4),(pcA(k), k = n1,n2)
           write(09,825) label(5),(mwA(k), k = n1,n2)
           write(09,820) label(6),(vlA(k), k = n1,n2)
           write(09,820) label(7),(ecA(k), k = n1,n2)
           write(09,820) label(8),(dfA(k), k = n1,n2)
           write(09,820) label(9),(dsA(k), k = n1,n2)
           write(09,840) label(69),(valnc(k), k = n1,n2)
           
           write(09,820) label(10),(pcB1(k), k = n1,n2)
           write(09,825) label(11),(mwB1(k), k = n1,n2)
           write(09,820) label(12),(vlB1(k), k = n1,n2)
           write(09,820) label(13),(ecB1(k), k = n1,n2)
           write(09,820) label(14),(dfB1(k), k = n1,n2)
           write(09,820) label(15),(dsB1(k), k = n1,n2)

           write(09,820) label(16),(pcB2(k), k = n1,n2)
           write(09,825) label(17),(mwB2(k), k = n1,n2)
           write(09,820) label(18),(vlB2(k), k = n1,n2)
           write(09,820) label(19),(ecB2(k), k = n1,n2)
           write(09,820) label(20),(dfB2(k), k = n1,n2)
           write(09,820) label(21),(dsB2(k), k = n1,n2)
           write(09,800) label(71),(satB2(k), k = n1,n2)
           
           write(09,820) label(22),(vlAA(k), k = n1,n2)
           write(09,820) label(23),(dfAA(k), k = n1,n2)
           write(09,820) label(24),(dsAA(k), k = n1,n2)

           write(09,820) label(25),(vlAB1(k), k = n1,n2)
           write(09,820) label(26),(dfAB1(k), k = n1,n2)
           write(09,820) label(27),(dsAB1(k), k = n1,n2)

           write(09,820) label(28),(vlAB2(k), k = n1,n2)
           write(09,820) label(29),(dfAB2(k), k = n1,n2)
           write(09,820) label(30),(dsAB2(k), k = n1,n2)

           write(09,820) label(31),(vlABB(k), k = n1,n2)
           write(09,820) label(32),(dfABB(k), k = n1,n2)
           write(09,820) label(33),(dsABB(k), k = n1,n2)

           write(09,820) label(34),(vlC2(k), k = n1,n2)
           write(09,820) label(35),(dfC2(k), k = n1,n2)
           write(09,820) label(36),(dsC2(k), k = n1,n2)

           write(09,820) label(37),(vlC3(k), k = n1,n2)
           write(09,820) label(38),(dfC3(k), k = n1,n2)
           write(09,820) label(39),(dsC3(k), k = n1,n2)

           write(09,820) label(40),(vlC4(k), k = n1,n2)
           write(09,820) label(41),(dfC4(k), k = n1,n2)
           write(09,820) label(42),(dsC4(k), k = n1,n2)

           write(09,820) label(43),(vlC5(k), k = n1,n2)
           write(09,820) label(44),(dfC5(k), k = n1,n2)
           write(09,820) label(45),(dsC5(k), k = n1,n2)

           write(09,820) label(46),(vlC6(k), k = n1,n2)
           write(09,820) label(47),(dfC6(k), k = n1,n2)
           write(09,820) label(48),(dsC6(k), k = n1,n2)

           write(09,820) label(49),(vlC7(k), k = n1,n2)
           write(09,820) label(50),(dfC7(k), k = n1,n2)
           write(09,820) label(51),(dsC7(k), k = n1,n2)
           write(09,810) ' '

           write(09,830) label(52),(kf_AA(k), k = n1,n2)
           write(09,830) label(53),(kr_AA(k), k = n1,n2)

           write(09,830) label(54),(kf_AB1(k), k = n1,n2)
           write(09,830) label(55),(kr_AB1(k), k = n1,n2)

           write(09,830) label(56),(kf_AB2(k), k = n1,n2)
           write(09,830) label(57),(kf_AB2(k), k = n1,n2)

           write(09,830) label(58),(kf_CC(k), k = n1,n2)
           write(09,830) label(59),(kf_CC(k), k = n1,n2)

           write(09,810) ' '
           write(09,840) label(60),(clen(k), k = n1,n2)
           write(09,840) label(61),(smpsz(k), k = n1,n2)
           write(09,840) label(62),(monq(k), k = n1,n2)
           write(09,840) label(63),(rnln(k), k = n1,n2)
           write(09,840) label(70),(prls(k), k = n1,n2)           
           write(09,820) label(64),(tim_cyc(k),k = n1,n2)
           write(09,820) label(65),(tim_inc(k),k = n1,n2)
           write(09,820) label(66),(del_max(k),k = n1,n2)
           write(09,840) label(67),(spread(k),k = n1,n2)
           write(09,840) label(68),(difcyc(k),k = n1,n2)
           write(09,820) label(72),(colvols(k),k = n1,n2)
           write(09,820) label(73),(florts(k),k = n1,n2)


 200       continue
           close (unit = 09)
           

 800       format(2x,a10,9a15)
 810       format(1a)
 820       format(2x,a10,9f15.6)
 825       format(2x,a10,9f15.0)
 830       format(2x,a10,9e15.2)
 840       format(2x,a10,9i15)


 900     format(//)

             return
             end


c-----------------------------------------------------------------------------



        subroutine MIGRATE(prn,p1,p2,sprd,cl,ss,
     x                     vl,dis,gs,saturate,pc,mw)

        !
        !     routine to translate protein along column with 
        !               gaussian dispersion
        !

        !     k:         component i.d.
        !     sprd:      spread (cells) for gaussian dispersion
        !     prfrs:     proteins zone -- leading edge supplemented
        !     vl:        constituent velocities
        !     dis:       dispersion factors
        !      gs:        gaussian table values
        !      fa:
        !      fb:
        !      fgz:
        !      p11:
        !      p3:
        !      p4:
        !      vx:
        !      dx:
        !      z1:
        !      k1:
        !      z2:
        !      vl1:
        !      ndx
        !      ndq:
        !      kv:
        !      p2inc:
        !      p1
        !      p2
        !
        !             09/11/91

        real fa,fb,fgz(2),vx,dx,gs(705),mw(0:15),pc(0:5),
     x     vl(0:15),dis(0:15)
        integer j,z1,k1,k,z2,vl1,ndx(2,0:1),ndq(2,0:1),kv
        integer sprd,prfrs,cl,ss,prn,saturate
        double precision p2inc,p1(-10:6000,0:15),p2(0:6000,0:15),
     x                 p11,p3,p4

          prfrs = cl + ss


        do 100  k = 0, prn


        do 300 j=prfrs,1-10*saturate,-1

         vx=vl(k)
           if (j .le. ss) vx = 1
          vl1=vx
          do 290 k1 = -sprd, sprd
           do 260 z1 = 1,2
            do 250 z2 = 0,1
             ndx(z1,z2) = z2+100*((1.5-z1)+k1-vl1+vx)/DIS(K)
              if (ndx(z1,z2) .gt. 350) ndx(z1,z2)=350
              if (ndx(z1,z2) .lt. -350) ndx(z1,z2)=-350
             ndx(z1,z2)=ndx(z1,z2)+351
250         continue
260        continue

          do 270 z1=1,2
           ndq(z1,0)=ndx(z1,0)
           ndq(z1,1)=ndx(z1,1)
            fr=ndx(z1,0)-ndq(z1,0)
            fa=gs(ndq(z1,0))
            fb=gs(ndq(z1,1))
            fgz(z1)=fa+fr*(fb-fa)
270       continue

          fgz1=fgz(1)-fgz(2)
          if (abs(fgz1) .lt. .0001) fgz1=0.
          idum=j+vl1-k1
           if (idum .gt. 0) then
             kv=k
             if (p1(j,k) .lt. .1e-20) p1(j,k)=0.
             p2(idum,kv)=(fgz1*p1(j,k))+p2(idum,kv)
           end if
290       continue
300      continue


         do 5 j=1,cl+SS
           p1(j,k)=p2(j,k)
           p2(j,k)=0.
 5       continue

 100     continue

         if(saturate .eq. 1) then

            do j =0,10
            p1(-j,3) = pc(3)/mw(3)
            end do

         end if


                  end
    


c-----------------------------------------------------------------------------


         subroutine MONITOR(prn,p1,p2,mqt,cl,mq)

          !
          !  routine to emulate sample elution
          !

          !   k:     component i.d.
          !   mqt:   length of monitor queue
          !   cl:    column length
          !   mq:    concentration (molar) of component
          !           k at cell j in the monitor queue
          !   p1:
          !   p2:
          !
          !         09/16/91
          !

       integer k,mqt,cl,prn
       real mq(0:200,0:15)
        double precision p1(-10:6000,0:15),p2(0:6000,0:15)


        do 100 k =0, prn

          do 2 j=mqt,2,-1
          mq(j,k)=mq(j-1,k)
2         continue

         mq(1,k)=0
        do 5 j=1,mqt
          mq(j,k)=mq(j,k)+p1(j+cl,k)
          p1(j+cl,k)=0.
5        continue

 100     continue
         end



c-----------------------------------------------------------------------------

       subroutine PRO_FILE(infile_pro,infiles,max_files)

       character*20 infile_pro
       character*25 infiles(100)
       integer max_files

       open (unit=03, file = infile_pro, status = 'old')

       n = 0
 199     n = n + 1
          read(03,900,end=99)  infiles(n)
           if (infiles(n)(1:1) .eq. ' ') goto 99

          go to 199
 99       max_files = n - 1


          close (unit = 03)



 900      format(a25)


              end

c-----------------------------------------------------------------------------


       subroutine PROTEINS(pro_conc,mol_wt,velocity,extinc,
     x                    diffus_coef,dispers,kf,kr,
     x                    protein_num,infile,valence,mess_t_inc)


              !
              ! Routine opens file (proteins_in.dat) and obtains
              ! values that represent two interacting proteins.
              ! kf/r(0) rate constants for self-association of molecule A
              ! kf/r(1) rate constants between molecule A and molecule B1.
              ! kf/r(2) rate constants between molecule A and molecule B2.
              ! kf/r(3) rate constants between Cs where C = A:A
              !
              !        05/24/93
              !

         character*20  infile
         integer protein_num, valence
         real pro_conc(0:15),mol_wt(0:15),velocity(0:15),
     x          extinc(0:15),diffus_coef(0:15),dispers(0:15),
     x          kf(0:3),kr(0:3)


         open(unit=03,file=infile,status='old')

c
c     Read concentrations, etc. for molecules A, B1, and B2.
c
         do j=1,3
            read(03,100) pro_conc(j),mol_wt(j),velocity(j),
     x        extinc(j),diffus_coef(j),dispers(j)
         end do

c
c     Read velocity, diffusion, and dispersion for dimer AA
c

         read (03,110) velocity(0), diffus_coef(0), dispers(0)

             read(03,110) (velocity(j), j=4,8)
             read(03,110) (velocity(j), j=9,13)

             read(03,120) (dispers(j), j=4,8)
             read(03,120) (dispers(j), j=9,13)

             read(03,120) (diffus_coef(j), j=4,8)
             read(03,120) (diffus_coef(j), j=9,13)

c
c     Read association and dissociation rates for:
c
c                         dimerization of molecule A: kf/r(0)
c                            interaction of A and B1: kf/r(1)
c                            interaction of A and B2: kf/r(2)
c                           aggregation of dimers AA: kr/r(3)
c

           do j=0,3   
             read(03,*) kf(j)
             read(03,*) kr(j)
           end do

           read(03,*) valence


          close(unit=03)

           protein_num = 5
          if (kf(3) .gt. 0.) protein_num = 13 
          if (valence .ne. 2) valence = 1
          if ((kf(0) .ne. 0) .or. (kf(2) .ne. 0) .or.
     x        (kf(1) .eq. 0)) valence = 1

          mess_t_inc = 0

100     format(f7.3,1x,f8.1,1x,4(f7.3,1x))
110     format(5(f6.2,1x))
120     format(5(f5.2,1x))
                   end



c-----------------------------------------------------------------------------

        subroutine PLOT_WRITE(n_rec,fil1,n_unit)
           !
           !
           !routine opens file (fil1//n_rec.out)
           !
           !

        character*4 fil1
        character*12 recfl


         write(unit=recfl, fmt = '(a4,i4,a4)'),fil1,n_rec,'.out'
         open (unit = n_unit, file=recfl, status='new')
         write(unit=recfl, fmt = '(a4,i4,a4)'),fil1,n_rec,'.mat'
         open (unit = n_unit+1, file=recfl, status='new')

                        end

c-----------------------------------------------------------------------------


       subroutine REPORT(ncycl,prel,prel_t,mqt,cl,ss,p1,pc,mw,
     x                count,rl,t_tot,prn,eluted,eluted_t,
     x                column_vol, flo_rate,
     x                position_label,saturate,valence)



c 
c      Writes elution values to terminal and to file
c
c      prel    === eluted protein
c      prel(0) === dimer of protein A
c      prel(1) === protein A
c      prel(2) === protein B1
c      prel(3) === protein B2    B1 and B2 compete for the same site on A
c      prel(4) === complex A:B1
c      prel(5) === complex A:B2
c      prel(6) === summed aggregates of A:A
c
       
        character*5 position_label
        integer ss,cl,rl,prn,saturate,valence
        real eluted(0:15),mw(0:15),pq(100),elu_time,t_tot,pc(0:15)
        double precision  prel(0:15), p1(-10:6000,0:15)

        if (ncycl .le. 1) then
           if (valence .eq. 1) then
             write(6,820) position_label
             write(7,820) position_label
           else
              write(6,825) position_label
              write(7,825) position_label
           end if
         end if



         count = 0
         p1_tot = 0




         if ((saturate .eq. 1) .and. (ncycl .lt. mqt)) then
           do j = mqt,mqt-ncycl,-1
            pq(j) = pc(3)
           end do
         end if



         do j = 0,prn
            p1_tot = p1_tot + prel(j)
           end do

         do j = mqt,2,-1
           pq(j) = pq(j-1)
         end do
           pq(1) = p1_tot



       if (count .gt. float(rl)/200.) count = 0.


      if (count .eq. 0) then

         if (position_label .eq. 'Mins.') then
          elu = t_tot*float(ncycl)/60.
          write(6,800) elu, (prel(m),m=0,6), pq(mqt)
         end if

         if (position_label .eq. 'Vol. ') then
          elu = t_tot*float(ncycl)*flo_rate/60.
          write(6,800) elu, (prel(m),m=0,6), pq(mqt)
         end if

         if (position_label .eq. 'Ve/Vt') then
          elu = t_tot*float(ncycl)*flo_rate/(60.*column_vol)
          write(6,800) elu, (prel(m),m=0,6), pq(mqt)
         end if


      end if


      count = count + 1

         write(07,800) elu, (prel(m),m=0,6), pq(mqt)
         write(08,805) elu, pq(mqt)
         write(09,810) pq(mqt)



 800   format(f7.2,7f9.4,2x,f9.4)
 805   format(f8.4,7f9.4,2x,f9.4)
 806   format(2x,i5,f8.4,7f9.4,2x,f9.4)
 810   format(1x,7f9.4,2x,f9.4)
 820   format(/,3X,a5,3X,'A:A',7X,'A',7x,'B1',7X,'B2',6X,'A:B1',
     x        5X,'A:B2', 4X,'(A:A)n',5X,'@MONTR')
 825   format(/,3X,a5,2X,'B:A:B',6X,'A',7x,'B1',7X,'B2',6X,'A:B1',
     x        5X,'A:B2', 4X,'(A:A)n',5X,'@MONTR')

                       end



c-----------------------------------------------------------------


         subroutine SATURAT(pc,mw,vl,p1,p2,saturate,prn,sprd,cl,ss,
     x                dis,gs,mqt,mq,ndif,dif)

         integer prn,sprd,cl,ss,saturate,mqt
         real mw(0:15),pc(0:5),vl(0:15),dis(0:15),gs(705),
     x        mq(0:200,0:15)
         double precision p1(-10:6000,0:15), p2(0:6000,0:15)
       

         if(saturate .eq. 1) then

            n_reps = 1.3*cl
            pr_load = pc(3)/mw(3)

            do j = 1,n_reps
               do k = -5,0
                 p1(0,3) = pr_load
               end do
               call MIGRATE(3,p1,p2,sprd,cl,ss,vl,dis,
     x                      gs,saturate,pc,mw)
               call MONITOR(prn,p1,p2,mqt,cl,mq)
               call DIFFUSION_1(ndif,p1,p2,j,prn,cl+ss,dif)

            end do


            end if
                     end


c-----------------------------------------------------------------


       subroutine SUMS (eluted,eluted_t,message_t_inc,t_inc)

         integer message_t_inc
         real eluted(0:15)
         double precision  t_inc


         write(6,100)  (eluted(j), j=0,6), eluted_t
         write(7,100)  (eluted(j), j=0,6), eluted_t

         write(6,850) t_inc
         write(7,850) t_inc

         close(unit = 07)
         close(unit = 08)

c
c     Clear < eluted > array to prepare for next run
c

         do j = 1,6
            eluted(j) = 0
         end do

         eluted_t = 0


 100     format(/,1x,'Sums',7f9.4,2x,f9.4)
 850     format(//,4x,'Note:  The value of the parameter < t_inc > was',
     x    ' changed to',/,' [', f12.8, ' ] seconds',
     x    ' to satisfy the criterion set by < delta_max >.')


                       end

c-----------------------------------------------------------------------------


       subroutine TERM_WRITE(cl,rl,mqt,ss,sprd,ndif,vl,
     x      dis,t_inc,delta_max,t_tot,pc,mw,ec,dif,kf,kr,n_rec,
     x      recfl,fil1,n_unit,valence)

       double precision  t_inc

           !
           !
           !  Routine writes documentation to terminal screen
           !  if n_unit = 06 or 05, or writes to file 
           !  (fil1//n_rec.chr) if n_unit = other value
           !
           !        09/14/91
           !
           !

        integer ss,cl,rl,sprd,valence
        real kf(0:3),kr(0:3),vl(0:15),dis(0:15),ec(0:15),mw(0:15),
     x       dif(0:15),pc(0:15),k12
        character*4 fil1
        character*9 daymonyr
        character*12 recfl

        if ((n_unit .ne. 5) .and. (n_unit .ne. 6)) then
         write(unit=recfl, fmt = '(a4,i4,a4)'),fil1,n_rec,'.chr'
         open (unit = n_unit, file=recfl, status='new')
        end if

        call date(daymonyr)


        write(n_unit,130) daymonyr


        run_time = t_tot*float(rl)/60.
        write(n_unit,100) cl-ss,rl,t_tot,run_time,t_inc,
     x   100*delta_max, mqt,ss,sprd,ndif,(kf(j),kr(j), j=0,3)


        write(n_unit,125) valence
        if (valence .eq. 1) then
          write(n_unit,110)
        else
          write(n_unit,112)
        end if

        WRITE (n_unit,120)
     x   VL(0), VL(1), VL(2), VL(3),
     X   DIS(0),DIS(1),DIS(2),DIS(3),
     X   PC(0), PC(1), PC(2), PC(3),
     X   MW(0), MW(1), MW(2), MW(3),
     X   EC(0), EC(1), EC(2), EC(3),
     X   DIF(0),DIF(1),DIF(2),DIF(3)

        write(n_unit,115)

        WRITE (n_unit,120)
     x   VL(4), VL(5),VL(6),VL(7),
     X   DIS(4),DIS(5),DIS(6),DIS(7),
     X   PC(4), PC(5), PC(6), PC(7),
     X   MW(4), MW(5), MW(6), MW(7),
     X   EC(4), EC(5), EC(6), EC(7),
     X   DIF(4),DIF(5),DIF(6),DIF(7)


100     format(/,5x,'COL LEN (cells): ',i4,
     x           7x,'RUN LEN(cycles): ',i4,
     x         /,5x,'T/CYCLE (secs):  ',f7.3,
     x          ,4x,'RUN TIME (mins)  ',f7.3,
     x         /,5x,'T_INC (secs):    ',f8.4,
     x           3x,'DELTA_MAX (%):   ',f5.1,
     x         /,5x,'MONITOR QUEUE:   ',i4,
     x           7x,'SAMPLE SIZE:     ',i4,
     x         /,5x,'SPRD +/-:        ',i4,
     x           7x,'DIF CYCLES:      ',i4,/,
     x         /,5x,'kf(0) =          ',e7.2,
     x           4x,'kr(0) =          ',e7.2,
     x         /,5x,'kf(1) =          ',e7.2,
     x           4x,'kr(1) =          ',e7.2,
     x         /,5x,'kf(2) =          ',e7.2,
     x           4x,'kr(2) =          ',e7.2,
     x         /,5x,'kf(3) =          ',e7.2,
     x           4x,'kr(3) =          ',e7.2,/)

 110    format(///,21x,'A:A',8x,'A',8x,'B1',8x,'B2')
 112    format(///,20x,'B:A:B',7x,'A',8x,'B1',8x,'B2')
 115    format(///,21x,'A:B1',6x,'A:B2',5x,'(A:A)2',5x,'(A:A)3')

 120    format(/'  VELOCITY    :',4F10.2,'  cells/cycle',
     X         /'  DISP.       :',4F10.2,'  cells/cycle',
     X /'  CONC.       :',4F10.2,'  mg/ml',
     X /'  MOLEC. WT.  :',4F10.0,'  Daltons',
     X /'  EXTN. COEF. :',4F10.2,'  a.u. / (mg/ml)',
     X /'  DIFF.       :',4F10.2,'  cells/cycle',/)

 125    format(/,2x,'Valence of A for B1 is:  ',i3)
 130    format(//,10x,a9,/)        

               end




c-----------------------------------------------------------------------------




















