subroutine cosat(nsize,t,p,qsat,qsurf_n2,qsurf_ch4)

  IMPLICIT NONE

!=======================================================================
!
!  CO mass mixing ratio at saturation (kg/kg) for a given pressure (Pa)
!   and Temperature (K) array
!
!=======================================================================


!   INPUT
      integer nsize 
      real t(nsize) , p(nsize)
!   OUTPUT
      real qsat(nsize),qsurf_n2(nsize),qsurf_ch4(nsize)
      INTEGER i

          
      do i=1,nsize

!          from Fray and schmitt fit by formulation L = 2.74e5 J/kg
           qsat(i)=0.1537*exp((28*274./8.314)*(1/68.1-1/t(i)))*100000*28.0/(28.0*p(i))

!          Raoult law if mixte in N2 ice and CH4 ice : assume 0.3% CO in N2
!          (Merlin) and 0.3% in CH4
!          if ((qsurf_n2(i).gt.0.001).or.(qsurf_ch4(i).gt.0.001)) then

!       Raoult law if mixte in N2 ice : assume 0.3% CO in N2
        if ((qsurf_n2(i).gt.0.001)) then
           qsat(i)=qsat(i)*0.3/100.
        endif

!       security :
        qsat(i)=min(qsat(i),0.99)
      enddo

  return
end subroutine cosat