Changeset 3782 for trunk/LMDZ.COMMON/libf

Timestamp:

May 28, 2025, 1:48:58 PM (3 weeks ago)

Author:

jbclement

Message:

PEM:

• Correction to identify sublimating ice (surface vs subsurface) at initialization
• Gathering the processes of "ice condensation" and "dust sedimentation" in the same framework of the layering algorithm
• Reordering the workflow to manage the different situations of the layering algorithm

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• layering_mod.F90 (modified) (10 diffs)
• pem.F90 (modified) (3 diffs)

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

@@ -673 +673 @@
 == 26/05/2025 == JBC
 Information about checking of CO2 mass balance.
+
+== 28/05/2025 == JBC
+- Correction to identify sublimating ice (surface vs subsurface) at initialization
+- Gathering the processes of "ice condensation" and "dust sedimentation" in the same framework of the layering algorithm
+- Reordering the workflow to manage the different situations of the layering algorithm

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

@@ -29 +29 @@
 real, parameter :: h_patchy_dust         = 0.05    ! Height under which a dust statum is considered patchy [m]
 
-! Lag layer parameters -> see Levrard et al. 2007
-real, parameter :: hmin_lag  = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate
+! Parameters for CO2 flux correction (inspired by Levrard et al. 2007)
+real, parameter :: hmin_lag  = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate [m]
 real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate
 
@@ -66 +66 @@
 !     > add_stratum
 !     > insert_stratum
-!     > rm_stratum
+!     > del_stratum
 !     > get_stratum
 !     > ini_layering
@@ -206 +206 @@
 !=======================================================================
 ! To remove a stratum in a layering
-SUBROUTINE rm_stratum(this,str)
+SUBROUTINE del_stratum(this,str)
 
 implicit none
@@ -218 +218 @@
 
 !---- Code
-! Protect the sub-surface strata from removing
+! Protect the sub-surface strata from deletion
 if (str%top_elevation <= 0.) return
 
@@ -237 +237 @@
 str => new_str
 
-END SUBROUTINE rm_stratum
+END SUBROUTINE del_stratum
 
 !=======================================================================
@@ -626 +626 @@
 
 !---- 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, h_toplag
-type(stratum), pointer :: currentloc
+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
 
 !---- Code
@@ -638 +638 @@
 zshift_surf = this%top%top_elevation
 zlag = 0.
-
+unexpected = .false.
+
+!-----------------------------------------------------------------------
 ! NOTHING HAPPENS
 if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then
     ! So we do nothing
 !-----------------------------------------------------------------------
-! DUST SEDIMENTATION with possibly ice condensation
-else if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice .and. dh_co2ice >= 0. .and. dh_h2oice >= 0.) then
-    write(*,*) '> Layering: dust sedimentation'
-    h_pore = dh_dust*dry_regolith_porosity/(1. - dry_regolith_porosity)
-    h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore
-    ! New stratum at the top of the layering
-    if (new_str) then
-        call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0.)
-        new_str = .false.
+! BUILDING A STRATUM (ice condensation and/or dust desimentation)
+else if (dh_co2ice >= 0. .and. dh_h2oice >= 0. .and. dh_dust >= 0.) then
+    ! Which porosity is considered?
+    if (dh_co2ice > dh_h2oice .and. dh_co2ice > dh_dust) then ! CO2 ice is dominant
+        h_pore = dh_co2ice*co2ice_porosity/(1. - co2ice_porosity)
+    else if (dh_h2oice > dh_co2ice .and. dh_h2oice > dh_dust) then ! H2O ice is dominant
+        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)
     else
-        this%top%top_elevation = this%top%top_elevation + h_tot
-        this%top%h_co2ice      = this%top%h_co2ice      + dh_co2ice
-        this%top%h_h2oice      = this%top%h_h2oice      + dh_h2oice
-        this%top%h_dust        = this%top%h_dust        + dh_dust
-        this%top%h_pore        = this%top%h_pore        + h_pore
-    endif
-!-----------------------------------------------------------------------
-! DUST LIFTING with possibly ice sublimation
-!else if (dh_dust < dh_co2ice .and. dh_dust < dh_h2oice .and. dh_co2ice <= 0. .and. dh_h2oice <= 0.) then
-else if (dh_dust < 0. .and. abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) 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 > eps .and. abs(dh_co2ice) < tol) .or. (current%h_h2oice > eps .and. abs(dh_h2oice) < tol)) then ! No, there is ice cementing the dust
-        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 CONDENSATION
-else if ((dh_co2ice > dh_dust .or. dh_h2oice > dh_dust) .and. dh_dust >= 0.) then
-    write(*,*) '> Layering: ice condensation'
-    ! Which porosity is considered?
-    if (dh_co2ice > dh_h2oice) then ! CO2 ice is dominant
-        h_pore = dh_co2ice*co2ice_porosity/(1. - co2ice_porosity)
-    else ! H2O ice is dominant
-        h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity)
+        unexpected = .true.
     endif
     h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore
@@ -702 +670 @@
     endif
 !-----------------------------------------------------------------------
-! ICE SUBLIMATION
-else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. abs(dh_dust) < tol) then
-    write(*,*) '> Layering: ice sublimation'
-    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))
-        h_co2ice_old = current%h_co2ice
-        h_h2oice_old = current%h_h2oice
-        h_dust_old = current%h_dust
-
-!~         ! How much is CO2 ice sublimation inhibited by the top dust lag?
-!~         h_toplag = 0.
-!~         if (associated(current%up)) then ! If there is an upper stratum
-!~             currentloc => current%up
-!~             ! Is the considered stratum a dust lag?
-!~             do while (.not. (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) .and. associated(currentloc%up))
-!~                 h_toplag = h_toplag + thickness(currentloc)
-!~                 currentloc => currentloc%up
-!~             enddo
-!~         endif
-!~         h2subl_co2ice_tot = h2subl_co2ice_tot*fred_subl**(h_toplag/hmin_lag)
-
-        ! 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
-        ! Is there H2O ice to sublimate?
-        h2subl_h2oice = 0.
-        if (h_h2oice_old > 0. .and. h2subl_h2oice_tot > 0.) then
-            h2subl_h2oice = min(h2subl_h2oice_tot,h_h2oice_old)
-            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 
-            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
-    if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0
+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
 !-----------------------------------------------------------------------
-!~ ! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION
-!~ else if (dh_co2ice < 0. .and. dh_h2oice > 0. .and. (abs(dh_dust) < abs(dh_co2ice) .or. abs(dh_dust) < abs(dh_h2oice))) then
-
-! TO DO????
-
+    ! ICE SUBLIMATION
+    else if (abs(dh_dust) < tol) then
+        write(*,*) '> Layering: ice sublimation'
+        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))
+            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
+            ! Is there H2O ice to sublimate?
+            h2subl_h2oice = 0.
+            if (h_h2oice_old > 0. .and. h2subl_h2oice_tot > 0.) then
+                h2subl_h2oice = min(h2subl_h2oice_tot,h_h2oice_old)
+                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 
+                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
+        if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0
+    else
+        unexpected = .true.
+    endif
 !-----------------------------------------------------------------------
-! UNUSUAL (WEIRD) SITUATION
+! 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
+
+!~     else
+        unexpected = .true.
+!~     endif
+!-----------------------------------------------------------------------
+! NOT EXPECTED SITUATION
 else
+    unexpected = .true.
+endif
+
+if (unexpected) then
     write(*,'(a)') 'Error: situation for the layering construction not managed!'
     write(*,'(a,es12.3)') '    > dh_co2ice [m/y] = ', dh_co2ice
     write(*,'(a,es12.3)') '    > dh_h2oice [m/y] = ', dh_h2oice
     write(*,'(a,es12.3)') '    > dh_dust   [m/y] = ', dh_dust, tol
-!~     error stop
+    error stop
 endif
 
@@ -786 +792 @@
 
 ! Remove the eroding dust lag if there is no dust anymore
-if (str%h_dust < tol) call rm_stratum(this,str)
+if (str%h_dust < tol) call del_stratum(this,str)
 
 END SUBROUTINE lift_dust_lag
@@ -829 +835 @@
 
 ! Remove the sublimating stratum if there is no ice anymore
-if (str%h_co2ice < tol .and. str%h_h2oice < tol) call rm_stratum(this,str)
+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

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

@@ -674 +674 @@
         enddo
     enddo
+    ! We put the sublimating tendency coming from subsurface ice into the overall tendency
+    where (zdqsdif_ssi_tot < 0.)
+        d_h2oice = zdqsdif_ssi_tot
+    end where
 endif
 do i = 1,ngrid
@@ -684 +688 @@
             co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
         endif
-        if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
-            is_h2oice_sublim_ini(i,islope) = .true.
-            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
+        if (d_h2oice(i,islope) < 0.) then
+            if (h2o_ice(i,islope) > 0.) then
+                is_h2oice_sublim_ini(i,islope) = .true.
+                h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
+            else if (h2o_ice_depth(i,islope) > 0.) then
+                is_h2oice_sublim_ini(i,islope) = .true.
+            endif
         endif
     enddo
@@ -692 +700 @@
 write(*,*) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini
 write(*,*) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf
-
-! We put the sublimating tendency coming from subsurface ice
-where (zdqsdif_ssi_tot > 0.)
-    d_h2oice = -zdqsdif_ssi_tot
-end where
 
 totmass_adsco2_ini = 0.