Changeset 4205 for trunk/LMDZ.VENUS

Timestamp:

Apr 29, 2026, 4:03:35 PM (9 days ago)

Author:

mlefevre

Message:

Venus PCM. Update of the polysulfur condensation scheme. ML

File:

• trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/polysulfur_condens.F (modified) (1 diff)

trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/polysulfur_condens.F

@@ -20 +20 @@
 !----------------------------------------------------------------------------
       real :: rat
+      ! Saturation vapor pressure of Sx (Pa)
       real :: psats2,psats3,psats4,psats8
+      ! Saturation vapor pressure of Sx (Pa) from Lyons et al., 2008
       real :: psats2l,psats3l,psats4l,psats8l
+      ! Saturation vapor pressure of Sx (Pa) from Zahnle et al., 2016
       real :: psats2z,psats3z,psats4z,psats8z
+      ! VMR equivalent of P_sat of Sx
       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
+      ! Ratios of Sx molec mass over Venus atm mean molec mass
+      real, parameter :: epss2 = 64.13/43.44
+      real, parameter :: epss3 = 96.198/43.44
+      real, parameter :: epss4 = 128.26/43.44
+      real, parameter :: epss8 = 256.52/43.44
+      ! Local total VMR of Sx
+      real, dimension(nblon,nblev) :: s2t,s3t,s4t,s8t
 
 ! >>> 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
-        
+      ! Init of total Sx
       s2t(:,:) = s2(:,:) + s2con(:,:)
       s3t(:,:) = s3(:,:) + s3con(:,:)
       s4t(:,:) = s4(:,:) + s4con(:,:)
       s8t(:,:) = s8(:,:) + s8con(:,:)
-      s2con(:,:) = 0.0
-      s3con(:,:) = 0.0
-      s4con(:,:) = 0.0
-      s8con(:,:) = 0.0
+      ! Condensed parts then set up to minimal value
+      s2con(:,:) = 1e-30
+      s3con(:,:) = 1e-30
+      s4con(:,:) = 1e-30
+      s8con(:,:) = 1e-30
+      ! Gas part set up to total
+      s2(:,:) = s2t(:,:)
+      s3(:,:) = s3t(:,:)
+      s4(:,:) = s4t(:,:)
+      s8(:,:) = s8t(:,:)
 
-      DO i = 1,nblon
+
+      ! Loop all over Venus
       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
+        DO i = 1,nblon
+          !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
+          !Zahnle 2016 does not have S3, S4. 
+          !Using the average of the ratio between the two sets
+          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)
+          ! Check the atm P condition to Psat S2
+          ! If condition is met, condensation
+          ! Otherwise  Sx_gas and Sx_cond stays as initialized          
+          if (PP(i,k).gt.psats2) then
+            qsats2 = psats2/(PP(i,k))
+            if (s2t(i,k).gt.qsats2) then
+              s2con(i,k) = s2t(i,k) - qsats2
+              s2(i,k) = qsats2
+            endif
+          endif
 
-        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
+          ! Check the atm P condition to Psat S3
+          ! If condition is met, condensation
+          ! Otherwise  Sx_gas and Sx_cond stays as initialized          
+          if (PP(i,k).gt.psats3) then
+            qsats3 = psats3/(PP(i,k))
+            if (s3t(i,k).gt.qsats3) then
+              s3con(i,k) = s3t(i,k) - qsats3
+              s3(i,k) = qsats3
+            endif
+          endif
+
+          ! Check the atm P condition to Psat S4
+          ! If condition is met, condensation
+          ! Otherwise  Sx_gas and Sx_cond stays as initialized          
+          if (PP(i,k).gt.psats4) then
+            qsats4 = psats4/(PP(i,k))
+            if (s4t(i,k).gt.qsats4) then
+              s4con(i,k) =s4t(i,k) - qsats4
+              s4(i,k) = qsats4
+            endif
+          endif
+
+          ! Check the atm P condition to Psat S8
+          ! If condition is met, condensation
+          ! Otherwise  Sx_gas and Sx_cond stays as initialized          
+          if (PP(i,k).gt.psats8) then
+            qsats8 = psats8/(PP(i,k))
+            if (s8t(i,k).gt.qsats8) then
+              s8con(i,k) = s8t(i,k) - qsats8
+              s8(i,k) = qsats8
+            endif
+          endif
+        ENDDO
       ENDDO