SUBROUTINE polysulfur_condens(
     + nblev, nblon,
     + TT,PP,
     + s2,s3,s4,s8,
     + s2con,s3con,s4con,s8con)

      USE chemparam_mod
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: nblon  ! nombre de points horizontaux
      INTEGER, INTENT(IN) :: nblev  ! nombre de couches verticales

!----------------------------------------------------------------------------
!     Ambient air state variables:
      REAL, INTENT(INOUT), DIMENSION(nblon,nblev) :: s2,s3,s4,s8,TT,PP
      REAL, INTENT(INOUT), DIMENSION(nblon,nblev) :: s2con,s3con,s4con,
     +                                               s8con      
!----------------------------------------------------------------------------
      INTEGER :: i,k
!----------------------------------------------------------------------------
      real :: rat
      real :: psats2,psats3,psats4,psats8
      real :: psats2l,psats3l,psats4l,psats8l
      real :: psats2z,psats3z,psats4z,psats8z
      real :: qsats2,qsats3,qsats4,qsats8
      real :: mms2,mms3,mms4,mms8
      real*8, dimension(nblon,nblev) :: s2t,s3t,s4t,s8t
      real*8 :: epss2,epss3,epss4,epss8

! >>> Program starts here:

      mms2 = 64.13 ! S2 molecular mass
      mms3 = 96.198 ! S3 molecular mass 
      mms4 = 128.26 ! S4 molecular mass
      mms8 = 256.52 ! S8 molecular mass
      epss2 = mms2/43.44
      epss3 = mms3/43.44
      epss4 = mms4/43.44
      epss8 = mms8/43.44
        
      s2t(:,:) = s2(:,:) + s2con(:,:)
      s3t(:,:) = s3(:,:) + s3con(:,:)
      s4t(:,:) = s4(:,:) + s4con(:,:)
      s8t(:,:) = s8(:,:) + s8con(:,:)
      s2con(:,:) = 0.0
      s3con(:,:) = 0.0
      s4con(:,:) = 0.0
      s8con(:,:) = 0.0

      DO i = 1,nblon
      DO k = 1,nblev
        !Lyons 2008 ; Rimmer 2021
        psats2l = 1e5*10**(7.0240 - 6091.2/(TT(i,k))) ! S2
        psats3l = 1e5*10**(6.3428 - 6202.2/(TT(i,k))) ! S3
        psats4l = 1e5*10**(6.0028 - 6047.5/(TT(i,k))) ! S4
        psats8l = 1e5*10**(4.1879 - 3269.1/(TT(i,k))) ! S8
        !Zahnle 2016
        if (TT(i,k).gt.413) then
                psats2z = 1e5*exp (16.1 - 14000/TT(i,k) )
                psats8z = 1e5*exp (9.6 - 7510/TT(i,k) )
        else
                psats2z = 1e5*exp (27. - 18500/TT(i,k) )
                psats8z = 1e5*exp (20. - 11800/TT(i,k) )
        endif
        psats2 = psats2z
        psats8 = psats8z        
        rat = (psats2z/psats2l + psats8z/psats8l)/2.
        psats3 = rat*psats3l
        psats4 = rat*psats4l

        qsats2 = epss2*psats2/(PP(i,k) - psats2)    
        qsats3 = epss3*psats3/(PP(i,k) - psats3)
        qsats4 = epss4*psats4/(PP(i,k) - psats4)
        qsats8 = epss8*psats8/(PP(i,k) - psats8)

        if (s2t(i,k).gt.qsats2) then
          s2con(i,k) = s2t(i,k) - qsats2
          s2(i,k) = qsats2
        endif
        if (s3t(i,k).gt.qsats3) then
          s3con(i,k) = s3t(i,k) - qsats3
          s3(i,k) = qsats3
        endif
        if (s4t(i,k).gt.qsats4) then
          s4con(i,k) = s4t(i,k) - qsats4
          s4(i,k) = qsats4
        endif  
        if (s8t(i,k).gt.qsats8) then 
          s8con(i,k) = s8t(i,k) - qsats8
          s8(i,k) = qsats8
        endif                   
       
      ENDDO
      ENDDO

      END SUBROUTINE polysulfur_condens