Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3784)
+++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3785)
@@ -681,2 +681,10 @@
 == 30/05/2025 == JBC
 Correction of 'h2o_ice_depth' update in the main loop.
+
+== 30/05/2025 == JBC
+- Rework of ice sublimation in the layering algorithm to sublimate ice and to lift dust at the same time
+- Addition of the possibility to make H2O ice or mixed dust/H2O ice lag when CO2 ice is sublimating
+- Addition of a lag resistance for CO2 ice sublimation
+- Consideration of a 'wet_lag_porosity' for H2O ice lag layer
+- Consideration of a threshold under which H2O ice lag layer is considered patchy
+- Deletion of unused function 'thickness_toplag'
Index: /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3784)
+++ /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3785)
@@ -19,13 +19,15 @@
 
 ! Physical parameters
-real, parameter :: d_dust                = 5.78e-2 ! Tendency of dust [kg.m-2.y-1]
-real, parameter :: dry_lag_porosity      = 0.4     ! Porosity of dust lag
-real, parameter :: dry_regolith_porosity = 0.4     ! Porosity of regolith
-real, parameter :: breccia_porosity      = 0.1     ! Porosity of breccia
-real, parameter :: co2ice_porosity       = 0.      ! Porosity of CO2 ice
-real, parameter :: h2oice_porosity       = 0.      ! Porosity of H2O ice
-real, parameter :: rho_dust              = 2500.   ! Density of dust [kg.m-3]
-real, parameter :: rho_rock              = 3200.   ! Density of rock [kg.m-3]
-real, parameter :: h_patchy_dust         = 0.05    ! Height under which a dust statum is considered patchy [m]
+real, parameter :: d_dust            = 5.78e-2 ! Tendency of dust [kg.m-2.y-1]
+real, parameter :: dry_lag_porosity  = 0.4     ! Porosity of dust lag
+real, parameter :: wet_lag_porosity  = 0.1     ! Porosity of water ice lag
+real, parameter :: regolith_porosity = 0.4     ! Porosity of regolith
+real, parameter :: breccia_porosity  = 0.1     ! Porosity of breccia
+real, parameter :: co2ice_porosity   = 0.      ! Porosity of CO2 ice
+real, parameter :: h2oice_porosity   = 0.      ! Porosity of H2O ice
+real, parameter :: rho_dust          = 2500.   ! Density of dust [kg.m-3]
+real, parameter :: rho_rock          = 3200.   ! Density of rock [kg.m-3]
+real, parameter :: h_patchy_dust     = 0.05    ! Height under which a dust lag is considered patchy [m]
+real, parameter :: h_patchy_ice      = 0.05    ! Height under which a H2O ice lag is considered patchy [m]
 
 ! Parameters for CO2 flux correction (inspired by Levrard et al. 2007)
@@ -77,5 +79,4 @@
 ! Procedures to get information about a stratum:
 !     > thickness
-!     > thickness_toplag
 !     > is_co2ice_str
 !     > is_h2oice_str
@@ -85,5 +86,5 @@
 ! Procedures for the algorithm to evolve the layering:
 !     > make_layering
-!     > lift_dust_lag
+!     > lift_dust
 !     > subl_ice_str
 !=======================================================================
@@ -294,9 +295,9 @@
     do i = index_breccia - 1,2,-1
         h_tot = layer_PEM(i) - layer_PEM(i - 1)
-        h_pore = h_tot*dry_regolith_porosity
+        h_pore = h_tot*regolith_porosity
         h_soil = h_tot - h_pore
         call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.)
     enddo
-    h_pore = layer_PEM(1)*dry_regolith_porosity
+    h_pore = layer_PEM(1)*regolith_porosity
     h_soil = layer_PEM(1) - h_pore
     call add_stratum(this,0.,0.,0.,h_soil,h_pore,0.)
@@ -482,28 +483,4 @@
 
 END FUNCTION thickness
-
-!=======================================================================
-! To get the top lag height
-FUNCTION thickness_toplag(this) RESULT(h_toplag)
-
-implicit none
-
-!---- Arguments
-type(layering), intent(in) :: this
-real :: h_toplag
-
-!---- Local variables
-type(stratum), pointer :: current
-
-!---- Code
-h_toplag = 0.
-current => this%top
-! Is the considered stratum a dust lag?
-do while (is_dust_lag(current) .and. associated(current))
-    h_toplag = h_toplag + thickness(current)
-    current => current%down
-enddo
-
-END FUNCTION thickness_toplag
 
 !=======================================================================
@@ -610,4 +587,93 @@
 !=======================================================================
 !=======================================================================
+! To lift dust in a stratum
+SUBROUTINE lift_dust(this,str,h2lift)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+real, intent(in) :: h2lift
+
+!---- Code
+! Update of properties in the eroding dust lag
+str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust)
+str%h_pore        = str%h_pore        - h2lift*str%h_pore/str%h_dust
+str%h_dust        = str%h_dust        - h2lift
+
+! Remove the eroding dust lag if there is no dust anymore
+if (str%h_dust < tol) call del_stratum(this,str)
+
+END SUBROUTINE lift_dust
+
+!=======================================================================
+! To sublimate ice in a stratum
+SUBROUTINE sublimate_ice(this,str,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+logical,                intent(inout) :: new_lag
+real,                   intent(inout) :: zlag
+real, intent(in) :: h2subl_co2ice, h2subl_h2oice
+
+!---- Local variables
+real                   :: h_ice, h2subl, h_pore, h_tot
+real                   :: hlag_dust, hlag_h2oice
+type(stratum), pointer :: current
+
+!---- Code
+h_ice = str%h_co2ice + str%h_h2oice
+h2subl = h2subl_co2ice + h2subl_h2oice
+
+! How much dust and H2O ice does ice sublimation involve to create the lag?
+hlag_dust = str%h_dust*h2subl/h_ice
+hlag_h2oice = 0.
+if (h2subl_co2ice > 0. .and. h2subl_h2oice < tol) hlag_h2oice = str%h_h2oice*h2subl/h_ice
+
+! Update of properties in the sublimating stratum
+str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hlag_dust - hlag_h2oice
+str%h_co2ice      = str%h_co2ice      - h2subl_co2ice
+str%h_h2oice      = str%h_h2oice      - h2subl_h2oice - hlag_h2oice
+str%h_dust        = str%h_dust        - hlag_dust
+str%h_pore        = str%h_pore        - h2subl*h_pore/h_ice
+
+! Correct the value of top_elevation for the upper strata
+current => str%up
+do while (associated(current))
+    current%top_elevation = current%down%top_elevation + thickness(current)
+    current => current%up
+enddo
+
+! Remove the sublimating stratum if there is no ice anymore
+if (str%h_co2ice < tol .and. str%h_h2oice < tol) call del_stratum(this,str)
+
+! Which porosity is considered?
+if (hlag_dust >= hlag_h2oice) then
+    h_pore = hlag_dust*dry_lag_porosity/(1. - dry_lag_porosity)
+else
+    h_pore = hlag_h2oice*wet_lag_porosity/(1. - wet_lag_porosity)
+endif
+h_tot = hlag_dust + hlag_h2oice + h_pore
+
+! New stratum above the current stratum as a lag layer
+if (new_lag) then
+    call insert_stratum(this,str,str%top_elevation + h_tot,0.,hlag_h2oice,hlag_dust,h_pore,0.)
+    new_lag = .false.
+else
+    str%up%top_elevation = str%up%top_elevation + h_tot
+    str%up%h_h2oice      = str%up%h_h2oice      + hlag_h2oice
+    str%up%h_dust        = str%up%h_dust        + hlag_dust
+    str%up%h_pore        = str%up%h_pore        + h_pore
+endif
+
+zlag = zlag + h_tot
+
+END SUBROUTINE sublimate_ice
+
+!=======================================================================
 ! Layering algorithm
 SUBROUTINE make_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current)
@@ -626,9 +692,10 @@
 
 !---- Local variables
-real    :: dh_co2ice, dh_h2oice, dh_dust
-real    :: h_co2ice_old, h_h2oice_old, h_dust_old
-real    :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
-real    :: h2lift, h2lift_tot, h_pore, h_tot
-logical :: unexpected
+real                   :: dh_co2ice, dh_h2oice, dh_dust
+real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
+real                   :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
+real                   :: h2lift, h2lift_tot, h_pore, h_tot, h_toplag, R_subl
+logical                :: unexpected
+type(stratum), pointer :: currentloc
 
 !---- Code
@@ -653,5 +720,5 @@
         h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity)
     else if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice) then ! Dust is dominant
-        h_pore = dh_dust*dry_regolith_porosity/(1. - dry_regolith_porosity)
+        h_pore = dh_dust*regolith_porosity/(1. - regolith_porosity)
     else
         unexpected = .true.
@@ -670,32 +737,33 @@
     endif
 !-----------------------------------------------------------------------
-else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. dh_dust <= 0.) then ! A stratum is retreating
-    ! DUST LIFTING
-    if (abs(dh_co2ice) < tol .and. abs(dh_co2ice) < tol) then
-        write(*,*) '> Layering: dust lifting'
-        h2lift_tot = abs(dh_dust)
-        do while (h2lift_tot > 0. .and. associated(current))
-            ! Is the considered stratum a dust lag?
-            if (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) then ! No, there is ice cementing the dust
-                dh_dust = 0. ! The tendency is set to 0
-                exit
-            else ! Yes, we can lift dust
-                h2lift = min(h2lift_tot,current%h_dust) ! How much dust can be lifted?
-                h2lift_tot = h2lift_tot - h2lift
-                call lift_dust_lag(this,current,h2lift)
-                if (h2lift_tot > 0.) current => current%down ! Still dust to be lifted so we move to the underlying stratum
-            endif
-        enddo
-        if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
-!-----------------------------------------------------------------------
-    ! ICE SUBLIMATION
-    else if (abs(dh_dust) < tol) then
-        write(*,*) '> Layering: ice sublimation'
+ ! RETREATING A STRATUM (ice sublimation and/or dust lifting)
+else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. dh_dust <= 0.) then
         h2subl_co2ice_tot = abs(dh_co2ice)
         h2subl_h2oice_tot = abs(dh_h2oice)
-        do while (h2subl_co2ice_tot > 0. .and. h2subl_h2oice_tot > 0. .and. associated(current))
+        h2lift_tot = abs(dh_dust)
+        do while (h2subl_co2ice_tot > 0. .and. h2subl_h2oice_tot > 0. .and. h2lift_tot > 0. .and. associated(current))
             h_co2ice_old = current%h_co2ice
             h_h2oice_old = current%h_h2oice
-    
+
+            ! How much is CO2 ice sublimation inhibited by the top dust lag?
+            R_subl = 1.
+            if (associated(current%up)) then ! If there is an upper stratum
+                h_toplag = 0.
+                currentloc => current%up
+                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
+                    h_toplag = h_toplag + thickness(currentloc)
+                    currentloc => currentloc%up
+                enddo
+                if (currentloc%h_h2oice > h_patchy_ice) then
+                    R_subl = 0.
+                else if (0. < currentloc%h_h2oice .and. currentloc%h_h2oice <= h_patchy_ice) then
+                    h_toplag = h_toplag + thickness(currentloc)
+                    R_subl = fred_subl**(h_toplag/hmin_lag)
+                else
+                    R_subl = fred_subl**(h_toplag/hmin_lag)
+                endif
+            endif
+           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
+
             ! Is there CO2 ice to sublimate?
             h2subl_co2ice = 0.
@@ -710,8 +778,25 @@
                 h2subl_h2oice_tot = h2subl_h2oice_tot - h2subl_h2oice
             endif
-            ! Sublimation
-            if (h2subl_co2ice > 0. .or. h2subl_h2oice > 0.) call subl_ice_str(this,current,h2subl_co2ice,h2subl_h2oice,zlag,new_lag)
-            ! Still ice to be sublimated and no more ice in the current stratum so we move to the underlying stratum
-            if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol) then 
+            ! Ice sublimation
+            if (h2subl_co2ice > 0. .or. h2subl_h2oice > 0.) call sublimate_ice(this,current,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)
+
+            ! Is there dust to lift?
+            h2lift = 0.
+            if (is_dust_lag(current) .and. h2lift_tot > 0.) then
+                h2lift = min(h2lift_tot,current%h_dust)
+                h2lift_tot = h2lift_tot - h2lift
+                ! Dust lifting
+                if (h2lift > 0.) call lift_dust(this,current,h2lift)
+            else if (associated(current%up)) then
+                if (is_dust_lag(current%up) .and. h2lift_tot > 0.) then
+                    h2lift = min(h2lift_tot,current%up%h_dust)
+                    h2lift_tot = h2lift_tot - h2lift
+                    ! Dust lifting
+                    if (h2lift > 0.) call lift_dust(this,current%up,h2lift)
+                endif
+            endif
+
+            ! Still ice to be sublimated or dust to be lifted and no more ice in the current stratum so we move to the underlying stratum
+            if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0. .or. h2lift_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol) then
                 current => current%down
                 new_lag = .true.
@@ -722,34 +807,11 @@
         if (h2subl_co2ice_tot > 0.) dh_co2ice = 0. ! No CO2 ice available anymore so the tendency is set to 0
         if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0
-    else
-        unexpected = .true.
-    endif
+        if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
 !-----------------------------------------------------------------------
 ! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION
 else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then
-!~     if (abs(dh_co2ice) > dh_h2oice .and. dh_co2ice < dh_dust .and. dh_dust <= 0.) then
-!~         h2subl_co2ice_tot = abs(dh_co2ice)
-!~         do while (h2subl_co2ice_tot > 0. associated(current))
-!~             h_co2ice_old = current%h_co2ice
-!~             h_h2oice_old = current%h_h2oice
-    
-!~             ! Is there CO2 ice to sublimate?
-!~             h2subl_co2ice = 0.
-!~             if (h_co2ice_old > 0. .and. h2subl_co2ice_tot > 0.) then
-!~                 h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old)
-!~                 h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice
-!~             endif
-!~             ! Sublimation
-!~             if (h2subl_co2ice > 0.) call subl_ice_str(this,current,h2subl_co2ice,0.,zlag,new_lag)
-!~             ! Still ice to be sublimated and no more ice in the current stratum so we move to the underlying stratum
-!~             if (h2subl_co2ice_tot > 0. .and. current%h_co2ice < tol .and. current%h_h2oice < tol) then 
-!~                 current => current%down
-!~                 new_lag = .true.
-!~             else
-!~                 exit
-!~             endif
-!~         enddo
-!~         if (h2subl_co2ice_tot > 0.) dh_co2ice = 0. ! No CO2 ice available anymore so the tendency is set to 0
-!~     else if (dh_h2oice > abs(dh_co2ice) .and. 0 <= dh_dust .and. dh_dust < dh_h2oice) then
+!~     if (abs(dh_co2ice) > dh_h2oice .and. dh_co2ice < dh_dust .and. dh_dust <= 0.) then ! CO2 sublimation is dominant
+
+!~     else if (dh_h2oice > abs(dh_co2ice) .and. 0 <= dh_dust .and. dh_dust < dh_h2oice) then ! H2O condensation is dominant
 
 !~     else
@@ -774,81 +836,3 @@
 END SUBROUTINE make_layering
 
-!=======================================================================
-! To lift dust in a stratum
-SUBROUTINE lift_dust_lag(this,str,h2lift)
-
-implicit none
-
-!---- Arguments
-type(layering),         intent(inout) :: this
-type(stratum), pointer, intent(inout) :: str
-real, intent(in) :: h2lift
-
-!---- Code
-! Update of properties in the eroding dust lag
-str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust)
-str%h_pore        = str%h_pore        - h2lift*str%h_pore/str%h_dust
-str%h_dust        = str%h_dust        - h2lift
-
-! Remove the eroding dust lag if there is no dust anymore
-if (str%h_dust < tol) call del_stratum(this,str)
-
-END SUBROUTINE lift_dust_lag
-
-!=======================================================================
-! To sublimate ice in a stratum
-SUBROUTINE subl_ice_str(this,str,h2subl_co2ice,h2subl_h2oice,zlag,new_lag)
-
-implicit none
-
-!---- Arguments
-type(layering),         intent(inout) :: this
-type(stratum), pointer, intent(inout) :: str
-logical,                intent(inout) :: new_lag
-real,                   intent(inout) :: zlag
-real, intent(in) :: h2subl_co2ice, h2subl_h2oice
-
-!---- Local variables
-real                   :: hsubl_dust, h_pore, h_ice, h_tot, h2subl
-type(stratum), pointer :: current
-
-!---- Code
-h_ice = str%h_co2ice + str%h_h2oice
-h2subl = h2subl_co2ice + h2subl_h2oice
-
-! How much dust does ice sublimation involve to create the lag?
-hsubl_dust = str%h_dust*h2subl/h_ice
-
-! Update of properties in the sublimating stratum
-str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hsubl_dust
-str%h_co2ice      = str%h_co2ice      - h2subl_co2ice
-str%h_h2oice      = str%h_h2oice      - h2subl_h2oice
-str%h_dust        = str%h_dust        - hsubl_dust
-str%h_pore        = str%h_pore        - h2subl*h_pore/h_ice
-
-! Correct the value of top_elevation for the upper strata
-current => str%up
-do while (associated(current))
-    current%top_elevation = current%down%top_elevation + thickness(current)
-    current => current%up
-enddo
-
-! Remove the sublimating stratum if there is no ice anymore
-if (str%h_co2ice < tol .and. str%h_h2oice < tol) call del_stratum(this,str)
-
-! New stratum above the current stratum as a dust lag
-if (hsubl_dust > 0.) then ! Only if the dust lag is building up
-    h_pore = hsubl_dust*dry_lag_porosity/(1. - dry_lag_porosity)
-    h_tot = hsubl_dust + h_pore
-    if (new_lag) then
-        call insert_stratum(this,str,str%top_elevation + h_tot,0.,0.,hsubl_dust,h_pore)
-        new_lag = .false.
-    else
-        str%up%top_elevation = str%up%top_elevation + h_tot
-    endif
-endif
-zlag = zlag + h_tot
-
-END SUBROUTINE subl_ice_str
-
 END MODULE layering_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3784)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3785)
@@ -69,7 +69,6 @@
 use writediagpem_mod,           only: writediagpem, writediagsoilpem
 use co2condens_mod,             only: CO2cond_ps
-use layering_mod,               only: ptrarray, stratum, layering, del_layering, make_layering, layering_algo, &
-                                      get_nb_str_max, nb_str_max, thickness_toplag, subsurface_ice_layering,   &
-                                      is_dust_lag, is_co2ice_str, is_h2oice_str
+use layering_mod,               only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, &
+                                      ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str
 use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
 use version_info_mod,           only: print_version_info
