Changeset 3842

Timestamp:

Jul 9, 2025, 11:52:15 AM (38 hours 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

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• deftank/run_PEM.def (modified) (3 diffs)
• layering_mod.F90 (modified) (4 diffs)
• pem.F90 (modified) (6 diffs)
• pemetat0.F90 (modified) (4 diffs)

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

@@ -728 +728 @@
 - Update of "visu_evol_layering.py", in particular to show value at cursor for 2D heatmaps.
 - Few cleanings.
+
+== 09/07/2025 == JBC
+- 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.

trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

@@ -78 +78 @@
 #---------- Ice management parameters ----------#
 # Amount of H2O ice to initialize the huge reservoir if the variable is not present in "startfi_PEM.nc". Default = 9200. kg.m-2 (= 10 m)
-# ini_huge_h2oice=1.e4
+# ini_huge_h2oice=9200.
 
 # Threshold to consider the amount of H2O ice as an infinite reservoir. Default = 460. kg.m-2 (= 0.5 m)
-# inf_h2oice_threshold=2.e3
+# inf_h2oice_threshold=460.
 
 # Do you want H2O frost to transform into perennial H2O ice? Default = .false.
@@ -87 +87 @@
 
 # Threshold to consider frost is becoming perennial H2O ice. Default = 460. kg.m-2 (= 0.5 m)
-# metam_h2oice_threshold=5.e-2
+# metam_h2oice_threshold=460.
 
 # Do you want the H2O ice to flow along subslope inside a cell? Default = .true.
@@ -96 +96 @@
 
 # Threshold to consider frost is becoming perennial CO2 ice. Default = 16500. kg.m-2 (= 10 m)
-# metam_co2ice_threshold=16.e3
+# metam_co2ice_threshold=16500.
 
 # Do you want the CO2 ice to flow along subslope inside a cell? Default = .true.

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

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

@@ -70 +70 @@
 use co2condens_mod,             only: CO2cond_ps
 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
+                                      ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str, print_layering
 use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
 use parse_args_mod,             only: parse_args
@@ -294 +294 @@
 
 ! CODE
+! Elapsed time with system clock
+call system_clock(count_rate = cr)
+call system_clock(c1)
+
+! Parse command-line options
+call parse_args()
+
+! Header
 write(*,*) '  *    .          .   +     .    *        .  +      .    .       .      '
 write(*,*) '           +         _______  ________  ____    ____      *           + '
@@ -302 +310 @@
 write(*,*) '            +       |_____|  |________||_____||_____|   +     .         '
 write(*,*) '  .      *          .   *       .   +       *          .        +      .'
-
-! Elapsed time with system clock
-call system_clock(count_rate = cr)
-call system_clock(c1)
-call parse_args()
 
 ! Some user info
@@ -618 +621 @@
     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
+    where (h2o_ice_depth > 0. .and. zdqsdif_ssi_tot < -1.e-10) d_h2oice = zdqsdif_ssi_tot
 endif
 do i = 1,ngrid
@@ -781 +782 @@
             do ig = 1,ngrid
                 call make_layering(layerings_map(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p)
+                !call print_layering(layerings_map(ig,islope))
                 co2_ice(ig,islope) = 0.
                 h2o_ice(ig,islope) = 0.
@@ -999 +1001 @@
         do islope = 1,nslope
             totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)*cell_area(ig)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-        enddo 
+        enddo
     enddo
     write(*,'(a,f8.3,a)') " > Relative total CO2 mass balance = ", 100.*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini), ' %'

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

@@ -23 +23 @@
 use layering_mod,               only: layering, array2stratif, nb_str_max, layering_algo, add_stratum, ini_layering
 use callkeys_mod,               only: startphy_file
+use glaciers_mod,               only: rho_co2ice, rho_h2oice
 
 #ifndef CPP_STD
@@ -54 +55 @@
 real, dimension(ngrid,nslope),           intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
 real, dimension(ngrid,nslope),           intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
-type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map             ! Layerings
+type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map       ! Layerings
 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
@@ -184 +185 @@
                     do islope = 1,nslope
                         call ini_layering(layerings_map(ig,islope))
-                        call add_stratum(layerings_map(ig,islope),ini_huge_h2oice,0.,ini_huge_h2oice,0.,0.,0.)
+                        call add_stratum(layerings_map(ig,islope),1.05*ini_huge_h2oice/rho_h2oice,0.,ini_huge_h2oice/rho_h2oice,0.05*ini_huge_h2oice/rho_h2oice,0.,0.)
                     enddo
                 else
                     do islope = 1,nslope
                         call ini_layering(layerings_map(ig,islope))
-                        if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),perennial_co2ice(ig,islope),perennial_co2ice(ig,islope),0.,0.,0.,0.)
+                        if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),1.05*perennial_co2ice(ig,islope)/rho_co2ice,perennial_co2ice(ig,islope)/rho_co2ice,0.,0.05*perennial_co2ice(ig,islope)/rho_co2ice,0.,0.)
                     enddo
                 endif
@@ -392 +393 @@
                 do islope = 1,nslope
                     call ini_layering(layerings_map(ig,islope))
-                    call add_stratum(layerings_map(ig,islope),ini_huge_h2oice,0.,ini_huge_h2oice,0.,0.,0.)
+                    print*, 'coucou', 1.05*ini_huge_h2oice/rho_h2oice, layerings_map(ig,islope)%top%top_elevation
+                    call add_stratum(layerings_map(ig,islope),1.05*ini_huge_h2oice/rho_h2oice,0.,ini_huge_h2oice/rho_h2oice,0.05*ini_huge_h2oice/rho_h2oice,0.,0.)
                 enddo
             else
                 do islope = 1,nslope
                     call ini_layering(layerings_map(ig,islope))
-                    if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),perennial_co2ice(ig,islope),perennial_co2ice(ig,islope),0.,0.,0.,0.)
+                    if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),1.05*perennial_co2ice(ig,islope)/rho_co2ice,perennial_co2ice(ig,islope)/rho_co2ice,0.,0.05*perennial_co2ice(ig,islope)/rho_co2ice,0.,0.)
                 enddo
             endif