
   
      implicit real*8(a-h,p,r-y),integer (z),integer*2 (i)          
c                                                                               
      common/ cst11 /fh2o,fco2    

      dimension ac(2000), ah(2000), xs(2000)

      delx = dlog(0.49999999d0)/150.
      x = 0.5 -   0.4999999d0

      do i  = 1, 150
         xs(i) = dexp(x)
         x = x + delx
      end do 

      do i = 151, 302
         xs(i) = 1. - dexp(x)
         x = x - delx
      end do 

      
      loop = 302 

      open (77,file='activities')
      write (77,*) 1
      write (77,*) 
      write (77,*) 2
      write (77,*) 10.,0.,10.,0.
      write (77,2000) loop*2,1,1,1,-1.,-1.
2000  format (/,/,/,/,4(I4,1x),/,/,2(f5.0,1x))


              
10    write (6,*) 'enter p(bar), t(k)'
      read (5,*) p,tk
c                      get stuff for activity calculation
         xc=0.0
         call mrk (p,tk,xc)
         fhp=dexp(fh2o)

         xc=1.0
         call mrk (p,tk,xc)
         fcp=dexp(fco2)    

         do i = 1, loop 
            xc= xs(i)
            call mrk (p,tk,xc)
               fhm=dexp(fh2o)
               fcm=dexp(fco2)
            
               ah(i) =fhm/fhp
               ac(i) =fcm/fcp 
           
         end do    
         write (77,1100) (ac(i),ah(i),xs(i)*p,i=1,loop)
       write (77,*) 1, 1,1, 1,1, 1,1
      write (77,*) 0,0,0,0,0
                      
1100  format (g14.5,2x,g14.5,2x,g14.5)
  
99    stop                          
      end                                                     
c                 
      subroutine mrk (p,t,xco2) 
c----------------------------------------------------------------------         
c                                                                               
c----------------------------------------------------------------------         
      implicit real*8 (a-g,o-y),integer (z),integer*2 (h-m)                     
c                                                                               
      common/ cst11 /f(2)         
c                                                                               
      dimension b(2),ad(2),x(2)                                           
c                                                                               
      data b/14.6,29.7/                                                         
c----------------------------------------------------------------------         
      x(2)=xco2                                                                 
      x(1)=1.d0-xco2                                                            
      tc=t-273.15d0                                                             
c                             compute dispersion term  for h2o                  
      if (tc.lt.465.0d0) ad(1)=7.4174d7+2.1325d5*tc-4.421d2*tc*tc               
     ++0.28623*tc**3                                                            
      if (tc.gt.465.d0) ad(1)=4.17278127d9-3.0825d7*tc+8.6359d4*tc*tc           
     +-1.05914d2*tc**3+4.79854d-2*tc**4                                         
      if (tc.gt.600.d0) ad(1)=1.382d8-4.936d4*tc                                
c                             compute dispersion term  for co2                  
      ad(2)=7.209d7-7.05d4*tc+21.29*tc**2                                       
      xk=dexp(-11.218+(6032.d0/t)-(2.782d6/t**2)+(4.708d8/t**3))                
c                             compute dispersion term  for h2o-co2              
c      tt=t/100.d0                                                              
c      xk=-16.92931d-4*tt+0.902564d-4*tt**2+82.31842d-4                         
      rt=83.14d0*t                                                              
      dsqrtt=dsqrt(t)                                                           
      ch=xk*6912.2595d0*t**2.5+79.267647d6                                      
c                             compute covolume bx and dispersion terms          
c                             for mixtures                                      
c                             (r**2=6912.2595)                                  
c                             (2.0*sqrt(ah2o*aco2)=79.26747d6)                  
c                             (ah2o=34.6d6 aco2=45.4d0)                         
      ch=xk*6912.2595d0*t**2.5+79.267647d6                                      
      bx=x(1)*14.6d0+xco2*29.7d0                                                
      aij=x(1)*x(1)*ad(1)+xco2*xco2*ad(2)+xco2*x(1)*ch                          
c                             solve for the molar volume of the mixture         
      c1=rt/p                                                                   
      c3=aij/p/dsqrtt                                                           
      c2=-bx*bx-c1*bx+c3                                                        
      c4=c2/3.d0-c1**2/9.d0                                                     
      c5=(4.5d0*c1*c2-c1**3)/27.d0-c3*bx/2.0                                    
      e1=c5*c5                                                                  
      e2=c4**3                                                                  
      if ((e1+e2).lt.0.d0) goto 10                                              
c                                          segment to solve for v               
c                                          if only one root.                    
      c9=dsqrt(e1+e2)                                                           
      c6=c9-c5                                                                  
      c7=-c5-c9                                                                 
      if (c6.lt.0.d0) goto 30                                                   
      c6=c6**.33333333d0                                                        
      goto 40                                                                   
30    c6=-((-c6)**.33333333d0)                                                  
40    if (c7.lt.0.d0) goto 50                                                   
      c7=c7**.33333333d0                                                        
      goto 60                                                                   
50    c7=-((-c7)**.33333333d0)                                                  
60    v=c6+c7+c1/3.d0                                                           
      goto 20                                                                   
c                                  segment to solve for v if                    
c                                  there are 3 roots.                           
10    if (c4.gt.0.d0) goto 100  
c
      write (6,*) '**warning** ptx conditions are in or near a two ' 
     *           ,'phase region predicted by the rk',xco2  
c
      phi=dsqrt(-e1/e2)                                                         
      goto 110                                                                  
100   phi=-dsqrt(-e1/e2)                                                        
c                                       get the arc cosine of phi in radians
c                                       for cdc and dec compilers use dacos,
c                                       for ibm and watfiv use darcos.
110   phi=dacos(phi)                                                           
      yzk=0.d0                                                                  
      if (xco2.lt.0.1.and.p.lt.300.) then 
         v = 1.d10
      else 
         v= -1.d10
      end if      
      do 90 h=1,3   
                                                            
      vv=2.d0*(dsqrt(-c4))*dcos(phi/3.d0+yzk)+c1/3.d0 
    
      if (xco2.lt.0.1.and.p.lt.300.) then 
         v = max(vv,v)
      else 
         v = min(vv,v)
      end if                       
                                                         
90    yzk=yzk+2.0943951d0                                                       
c                                       segment to compute the                  
c                                       h2o and co2 fugacities.                 
20    d2=dlog((v+bx)/v)                                                         
      d1=rt*dsqrtt*bx                                                           
      d3=d2-bx/(bx+v)                                                           
      d4=v-bx                                                                   
      d5=aij*d3/d1/bx                                                           
      d6=dlog(rt/d4)                                                            
c                                                                               
      do 70 l=1,2                                                               
      if (x(l).gt.0.d0) goto 80                                                 
c                                                                               
      f(l)=0.d0                                                                 
      goto 70                                                                   
c                                                                               
80    aj2=2.d0*x(l)*ad(l)+(1.d0-x(l))*ch                                        
      f(l)=dlog(x(l))+b(l)/d4-aj2/d1*d2+b(l)*d5+d6                              
70    continue                                                                  
c                                                                               
      return                                                                    
      end                                                                       
                     
      block data                                                                
c----------------------------------------------------                           
c                                                                               
c----------------------------------------------------                           
      implicit real*8 (a-g,o-y),integer (z),integer*2 (h-m)                     
c                                                                               
      common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ mnrvar /pp,tt,yy,rr              
     *      / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9                 
      data r,rr/8.3144126d00,83.14d0/                                           
c                              the following data statement assigns             
c                              logical unit numbers for i/o.                    
c                              for interact (batch modes):                      
c     data n1,n2,n3,n4,n5,n6,n7,n8,n9/21,22,6,24,25,26,5,6,29/                  
c                              for cms (interactive modes):                     
      data n1,n2,n3,n4,n5,n6,n7,n8,n9/21,22,23,24,25,26,27,28,29/               
c                                                                               
      end                                                                       
