Changeset 4074 for trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90

Timestamp:

Feb 17, 2026, 2:45:53 PM (10 days ago)

Author:

jbclement

Message:

PEM:

• Correct management of H2O ice tendency in 1D when there is not enough ice anymore.
• Clean initialization of allocatable module arrays (especially needed when no slope)
• One more renaming for consistency + few small updates thoughout the code.

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90 (modified) (5 diffs)

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

@@ -62 +62 @@
 if (.not. allocated(is_h2o_perice_PCM)) allocate(is_h2o_perice_PCM(ngrid))
 if (.not. allocated(co2_perice_PCM)) allocate(co2_perice_PCM(ngrid,nslope))
+is_h2o_perice_PCM(:) = .false.
+co2_perice_PCM(:,:) = 0._dp
 
 END SUBROUTINE ini_surf_ice
@@ -359 +361 @@
 
 if (ngrid == 1) then ! In 1D
-    h2o_ice = h2o_ice + d_h2oice*dt
-    zshift_surf = d_h2oice*dt/rho_h2oice
+    where (d_h2oice(:,:) < 0._dp .and. abs(d_h2oice(:,:)*dt) > h2o_ice(:,:)) ! Not enough ice to sublimate everything
+        ! We sublimate what we can
+        d_h2oice_new(:,:) = h2o_ice(:,:)/dt
+        ! It means the tendency is zero next time
+        d_h2oice(:,:) = 0._dp
+    else where
+        d_h2oice_new(:,:) = d_h2oice(:,:)
+    end where
 else ! In 3D
     call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit)
     if (stopcrit%stop_code() > 0) return
-
     call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
-    h2o_ice = h2o_ice + d_h2oice_new*dt
-    zshift_surf = d_h2oice_new*dt/rho_h2oice
 end if
+
+h2o_ice = h2o_ice + d_h2oice_new*dt
+zshift_surf = d_h2oice_new*dt/rho_h2oice
 
 END SUBROUTINE evolve_h2oice
@@ -408 +416 @@
 ! ---------------
 integer(di) :: i, islope
-real(qp)    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice, d_target ! Balance variables
+real(qp)    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice ! Balance variables
+real(dp)    :: d_target
 
 ! CODE
@@ -421 +430 @@
     do islope = 1,nslope
         if (d_h2oice(i,islope) > 0._dp) then ! Condensing
-            d_h2oice_new(i,islope) = d_h2oice_new(i,islope)*real(S_target_cond_h2oice/S_atm_2_h2oice,dp)
+            d_h2oice_new(i,islope) = d_h2oice(i,islope)*real(S_target_cond_h2oice/S_atm_2_h2oice,dp)
         else if (d_h2oice(i,islope) < 0._dp) then ! Sublimating
             d_target = d_h2oice(i,islope)*real(S_target_subl_h2oice/S_h2oice_2_atm,dp)
-            if (abs(d_target*dt) < h2o_ice(i,islope)) then ! Enough ice to sublimate everything
-                d_h2oice_new(i,islope) = real(d_target,dp)
+            if (abs(d_target*dt) <= h2o_ice(i,islope)) then ! Enough ice to sublimate everything
+                d_h2oice_new(i,islope) = d_target
             else ! Not enough ice to sublimate everything
                 ! We sublimate what we can
@@ -438 +447 @@
 end do
 
-! We need to remove this ice unable to sublimate from places where ice condenseds in order to keep balance
+! We need to remove this ice unable to sublimate from places where ice condensed in order to keep balance
 where (d_h2oice_new > 0._dp) d_h2oice_new = d_h2oice_new*real((S_target_cond_h2oice - S_ghostice)/S_target_cond_h2oice,dp)