Changeset 3424 for trunk/LMDZ.COMMON/libf/evolution

Timestamp:

Sep 2, 2024, 5:32:52 PM (3 months ago)

Author:

jbclement

Message:

PEM:
Modification to the layering algorithm to make dust lag prevent further sublimation.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• layering_mod.F90 (modified) (8 diffs)

trunk/LMDZ.COMMON/libf/evolution/changelog.txt

@@ -413 +413 @@
 == 26/08/2024 == JBC
 Correction for the launching script due to r3403: the job scheduler detection was missing in the case of a new cycle.
+
+== 02/09/2024 == JBC
+Modification to the layering algorithm to make dust lag prevent further sublimation.

trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90

@@ -27 +27 @@
 real, parameter :: rho_dust = 2500.             ! Density of dust [kg.m-3]
 real, parameter :: rho_rock = 3200.             ! Density of rock [kg.m-3]
+
+! Lag layer parameters -> see Levrard et al. 2007
+real, parameter :: hmin_lag = 0.5  ! Minimal height of the lag deposit to reduce the sblimation rate
+real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate
 
 ! Stratum = node of the linked list
@@ -465 +469 @@
 
 !---- Local variables
-real :: htend_co2ice, htend_h2oice, htend_dust
-real :: h_co2ice_old, h_h2oice_old, h_dust_old
-real :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
+real                   :: htend_co2ice, htend_h2oice, htend_dust
+real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
+real                   :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
+real                   :: h_toplag
+type(stratum), pointer :: currentloc
 
 !---- Code
@@ -545 +551 @@
             h_h2oice_old = current1%h2oice_volfrac*current1%thickness
             h_dust_old = current1%dust_volfrac*current1%thickness
+
+            ! How much is CO2 ice sublimation inhibited by the top dust lag?
+            currentloc => current1%up
+            h_toplag = 0.
+            do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                h_toplag = h_toplag + currentloc%thickness
+                currentloc => currentloc%up
+            enddo
+            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
+
             ! How much CO2 ice can sublimate in the considered stratum?
             if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
@@ -550 +566 @@
                 h2subl_tot = 0.
                 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
-            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we stop here
-                !current1 => current1%down ! Move to the underlying stratum
-                !new_lag1 = .true.
-                exit
+            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
+                ! Is the considered stratum a dust lag?
+                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
+                    current1 => current1%down
+                    new_lag1 = .true.
+                else ! No, so we stop here
+                    exit
+                endif
             else ! Only a fraction can sublimate and so we move to the underlying stratum
                 h2subl = h_co2ice_old
@@ -587 +607 @@
             h_h2oice_old = current1%h2oice_volfrac*current1%thickness
             h_dust_old = current1%dust_volfrac*current1%thickness
+
+            ! How much is CO2 ice sublimation inhibited by the top dust lag?
+            currentloc => current1%up
+            h_toplag = 0.
+            do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                h_toplag = h_toplag + currentloc%thickness
+                currentloc => currentloc%up
+            enddo
+            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
+
             ! How much CO2 ice can sublimate in the considered stratum?
             if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
@@ -592 +622 @@
                 h2subl_tot = 0.
                 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
-            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we stop here
-                !current1 => current1%down ! Move to the underlying stratum
-                !new_lag1 = .true.
-                exit
+            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
+                ! Is the considered stratum a dust lag?
+                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
+                    current1 => current1%down
+                    new_lag1 = .true.
+                else ! No, so we stop here
+                    exit
+                endif
             else ! Only a fraction can sublimate and so we move to the underlying stratum
                 h2subl = h_co2ice_old
@@ -614 +648 @@
             h_h2oice_old = current2%h2oice_volfrac*current2%thickness
             h_dust_old = current2%dust_volfrac*current2%thickness
+
+            ! How much is CO2 ice sublimation inhibited by the top dust lag?
+            currentloc => current2%up
+            h_toplag = 0.
+            do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                h_toplag = h_toplag + currentloc%thickness
+                currentloc => currentloc%up
+            enddo
+            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
+
             ! How much H2O ice can sublimate in the considered stratum?
             if (h2subl_tot - h_h2oice_old <= eps) then ! There is enough to sublimate
@@ -619 +663 @@
                 h2subl_tot = 0.
                 call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
-            else if (h_h2oice_old < eps) then ! There is nothing to sublimate so we stop here
-                !current2 => current1%down ! Move to the underlying stratum
-                !new_lag2 = .true.
-                exit
+            else if (h_h2oice_old < eps) then ! There is nothing to sublimate
+                ! Is the considered stratum a dust lag?
+                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
+                    current1 => current1%down
+                    new_lag1 = .true.
+                else ! No, so we stop here
+                    exit
+                endif
             else ! Only a fraction can sublimate and so we move to the underlying stratum
                 h2subl = h_h2oice_old