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

Timestamp:

Mar 25, 2026, 11:19:02 AM (3 weeks ago)

Author:

jbclement

Message:

PEM:

• Fix outputs "diagevo.nc" for 3D data.
• Fix sign in computing exchanges due to adsorption/ice table and in balancing H2O flux from/into atmosphere.
• Few cleanings.

JBC

File:

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

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

@@ -378 +378 @@
 end if
 
-h2o_ice = h2o_ice + d_h2oice_new*dt
-zshift_surf = d_h2oice_new*dt/rho_h2oice
+h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice_new(:,:)*dt
+zshift_surf(:,:) = d_h2oice_new(:,:)*dt/rho_h2oice
 
 END SUBROUTINE evolve_h2oice
@@ -385 +385 @@
 
 !=======================================================================
-SUBROUTINE 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)
+SUBROUTINE balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_balanced)
 !-----------------------------------------------------------------------
 ! NAME
@@ -414 +414 @@
 real(dp), dimension(:,:), intent(in)    :: h2o_ice
 real(dp), dimension(:,:), intent(inout) :: d_h2oice
-real(dp), dimension(:,:), intent(out)   :: d_h2oice_new
+real(dp), dimension(:,:), intent(out)   :: d_balanced
 
 ! LOCAL VARIABLES
 ! ---------------
 integer(di) :: i, islope
-real(qp)    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice ! Balance variables
+real(qp)    :: S_corr_subl, S_corr_cond, S_target_subl, S_target_cond, S_ghostice ! Balance variables
 real(dp)    :: d_target
 
 ! CODE
 ! ----
-S_target = (S_atm_2_h2o + S_h2o_2_atm)/2._dp
-S_target_cond_h2oice = S_atm_2_h2oice + S_target - S_atm_2_h2o
-S_target_subl_h2oice = S_h2oice_2_atm + S_target - S_h2o_2_atm
-
-d_h2oice_new = 0._dp
+! Compute global targets
+S_corr_cond = (S_h2o_2_atm - S_atm_2_h2o)/2._qp ! Correction = deviation to the mean
+S_corr_subl = (S_atm_2_h2o - S_h2o_2_atm)/2._qp ! Correction = deviation to the mean
+S_target_cond = abs(S_atm_2_h2oice + S_corr_cond)
+S_target_subl = abs(S_h2oice_2_atm + S_corr_subl)
+
+! Compute balanced tendencies with positivity limiter
+d_balanced(:,:) = 0._dp
 S_ghostice = 0._qp
 do i = 1,ngrid
     do islope = 1,nslope
-        if (d_h2oice(i,islope) > 0._dp) then ! Condensing
-            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) = d_target
-            else ! Not enough ice to sublimate everything
-                ! We sublimate what we can
-                d_h2oice_new(i,islope) = h2o_ice(i,islope)/dt
-                ! It means the tendency is zero next time
+        if (d_h2oice(i,islope) > 0._dp) then ! Condensation
+            d_balanced(i,islope) = d_h2oice(i,islope)*real(S_target_cond/S_atm_2_h2oice,dp)
+        else if (d_h2oice(i,islope) < 0._dp) then ! Sublimation
+            d_target = d_h2oice(i,islope)*real(S_target_subl/S_h2oice_2_atm,dp)
+            if (abs(d_target*dt) <= h2o_ice(i,islope)) then ! Enough ice to sublimate
+                d_balanced(i,islope) = d_target
+            else ! Not enough ice to sublimate
+                ! Sublimate all the available ice
+                d_balanced(i,islope) = -h2o_ice(i,islope)/dt
+                ! If fully depleted, zero tendency for next time
                 d_h2oice(i,islope) = 0._dp
-                ! We compute the amount of H2O ice that we could not make sublimate
-                S_ghostice = S_ghostice + abs(d_target*dt) - h2o_ice(i,islope)
+                ! Compute the amount of H2O ice unable to sublimate
+                S_ghostice = S_ghostice + real(abs(d_target*dt),qp) - real(h2o_ice(i,islope),qp)
             end if
         end if
@@ -450 +453 @@
 end do
 
-! 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)
+! Enforce conservation removing the ghost ice from places where ice condensed
+where (d_balanced(:,:) > 0._dp) d_balanced(:,:) = d_balanced(:,:)*real((S_target_cond - S_ghostice)/S_target_cond,dp)
 
 END SUBROUTINE balance_h2oice_reservoirs