Changeset 3842 for trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90

Timestamp:

Jul 9, 2025, 11:52:15 AM (6 days ago)

Author:

jbclement

Message:

PEM:

• Handling the situation in the layering algorithm when CO2 ice sublimation and H2O ice condensation happen simultaneously.
• Correction of the incorporation of the sub-surface sublimation tendency into the overall tendency to be more robust.
• Revision of the layering initialization in the case where there is no "startpem.nc" file.

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (modified) (4 diffs)

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

@@ -598 +598 @@
 else if (current%h_co2ice > 0. .and. current%h_co2ice > current%h_h2oice) then ! No, there is more CO2 ice than H2O ice
     h_toplag = 0. ! Because it matters only for H2O ice
-    if (h_toplag < h_patchy_ice) then ! But the dust lag is too thin to considered ice as subsurface ice
-        co2_ice = current%h_co2ice
-    endif
+    if (h_toplag < h_patchy_ice) co2_ice = current%h_co2ice ! But the dust lag is too thin to considered ice as subsurface ice
 endif
 
@@ -713 +711 @@
 
 !---- Local variables
-real                   :: dh_co2ice, dh_h2oice, dh_dust
+real                   :: dh_co2ice, dh_h2oice, dh_dust, dh_dust_co2, dh_dust_h2o
 real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
 real                   :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
@@ -724 +722 @@
 dh_h2oice = d_h2oice/rho_h2oice
 dh_dust = 0.
-if (dh_h2oice >= 0. .and. dh_co2ice >= 0.) then ! To put a dust sedimentation tendency only when ice is condensing
+if (dh_h2oice >= 0.) then ! To put a dust sedimentation tendency only when ice is condensing
     if (impose_dust_ratio) then
-        dh_dust = dust2ice_ratio*(dh_h2oice + dh_co2ice)
-    else 
+        if (dh_co2ice >= 0.) then
+            dh_dust = dust2ice_ratio*(dh_h2oice + dh_co2ice)
+        else
+            dh_dust = dust2ice_ratio*dh_h2oice
+        endif
+    else
         dh_dust = d_dust/rho_dust
     endif
@@ -838 +840 @@
 !-----------------------------------------------------------------------
 ! 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 ! 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
-        unexpected = .true.
-!~     endif
+else if (dh_co2ice <= 0. .and. dh_h2oice >= 0.) then
+    if (dh_dust >= 0.) then
+        dh_dust_co2 = 0.
+        dh_dust_h2o = dh_dust
+    else
+        dh_dust_co2 = dh_dust
+        dh_dust_h2o = 0.
+    endif
+    if (abs(dh_co2ice) > dh_h2oice) then ! CO2 ice sublimation is dominant
+        ! CO2 ICE SUBLIMATION: retreating a stratum
+        h2subl_co2ice_tot = abs(dh_co2ice)
+        h2lift_tot = abs(dh_dust_co2)
+        do while ((h2subl_co2ice_tot > 0. .or. h2lift_tot > 0.) .and. associated(current))
+            h_co2ice_old = current%h_co2ice
+
+            ! 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_sublco2**(h_toplag/hmin_lag_co2)
+                else
+                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
+                endif
+            endif
+           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
+
+            ! 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
+            ! Ice sublimation
+            if (h2subl_co2ice > 0.) call sublimate_ice(this,current,h2subl_co2ice,0.,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. h2lift_tot > 0.) .and. current%h_co2ice < 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 (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
+
+        ! H2O ICE CONDENSATION: building a stratum
+        h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity)
+        h_tot = dh_h2oice + dh_dust_h2o + h_pore
+        ! New stratum at the top of the layering
+        if (new_str) then
+            call add_stratum(this,this%top%top_elevation + h_tot,0.,dh_h2oice,dh_dust_h2o,h_pore,0.)
+            new_str = .false.
+        else
+            this%top%top_elevation = this%top%top_elevation + h_tot
+            this%top%h_h2oice      = this%top%h_h2oice      + dh_h2oice
+            this%top%h_dust        = this%top%h_dust        + dh_dust_h2o
+            this%top%h_pore        = this%top%h_pore        + h_pore
+        endif
+    else ! H2O ice condensation is dominant
+        ! H2O ICE CONDENSATION: building a stratum
+        h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity)
+        h_tot = dh_h2oice + dh_dust_h2o + h_pore
+        ! New stratum at the top of the layering
+        if (new_str) then
+            call add_stratum(this,this%top%top_elevation + h_tot,0.,dh_h2oice,dh_dust_h2o,h_pore,0.)
+            new_str = .false.
+        else
+            this%top%top_elevation = this%top%top_elevation + h_tot
+            this%top%h_h2oice      = this%top%h_h2oice      + dh_h2oice
+            this%top%h_dust        = this%top%h_dust        + dh_dust_h2o
+            this%top%h_pore        = this%top%h_pore        + h_pore
+        endif
+
+        ! CO2 ICE SUBLIMATION: retreating a stratum
+        h2subl_co2ice_tot = abs(dh_co2ice)
+        h2lift_tot = abs(dh_dust_co2)
+        do while ((h2subl_co2ice_tot > 0. .or. h2lift_tot > 0.) .and. associated(current))
+            h_co2ice_old = current%h_co2ice
+
+            ! 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_sublco2**(h_toplag/hmin_lag_co2)
+                else
+                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
+                endif
+            endif
+           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
+
+            ! 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
+            ! Ice sublimation
+            if (h2subl_co2ice > 0.) call sublimate_ice(this,current,h2subl_co2ice,0.,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. h2lift_tot > 0.) .and. current%h_co2ice < 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 (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
+    endif
 !-----------------------------------------------------------------------
 ! NOT EXPECTED SITUATION